## Abstract

In this paper, I review several free boundary problems that arise in the mathematical modelling of biological processes. The biological topics are quite diverse: cancer, wound healing, biofilms, granulomas and atherosclerosis. For each of these topics, I describe the biological background and the mathematical model, and then proceed to state mathematical results, including existence and uniqueness theorems, stability and asymptotic limits, and the behaviour of the free boundary. I also suggest, for each of the topics, open mathematical problems.

## 1. Introduction

Models in mathematical biology aim to explain biological processes and then serve as a tool to promote new research in biology. When the predictions of a mathematical model are in agreement with experimental data, the model can then be used to perform *in silico* experiments and to suggest new hypotheses. However, the biological processes are typically extremely complicated, and most mathematical models can only provide just a ‘proof of concept’; this is nevertheless useful as an initial step in the development of a more comprehensive model. Historically, applications in the physical sciences have been a driving force in the development of the mathematical sciences. The same is beginning to happen with the biological sciences. Biology-inspired theorems and theories are beginning to emerge in a number of diverse fields of mathematics. Quite often, even simple models give rise to challenging mathematical questions.

Free boundary problems have a long history, with most of them coming from physics and engineering. More recently, with the growing research in mathematical biology, free boundary problems are emerging in increasing number. Such problems arise when a live tissue begins to grow or shrink abnormally. In this paper, I consider mathematical models of free boundary problems that arise in tumour growth, wound healing, biofilms, granulomas and atherosclerosis. I give a brief biological introduction to each of the models; some of the model simulations actually agree quantitatively with experimental data. However, my emphasis here is on the mathematical questions that arise from the models. I report on mathematical results and suggest open problems.

## 2. Tumour growth

Denote a tumour region in by *Ω*(*t*). Within *Ω*(*t*), there are both tumour cells as well as other healthy normal cells, including, in particular, immune cells that attempt to slow down tumour growth, or even eliminate it. There are also many proteins within the tumour, secreted by the cells, which communicate signals of enhancement or inhibition of growth. A model that attempts to describe the interactions among these variables consists of a system of partial differential equations (PDEs) in *Ω*(*t*), which, in vector form, can be written as
Here *X*=(*X*_{1},…,*X*_{m}), ** u** is a velocity field,

*D*

_{X}is a diffusion matrix and

*f*(

*X*) is a nonlinear function which includes the biochemical communications among the variables; more precisely,

**may vary from one species to another.**

*u*We begin with an extremely simple model in which the tumour contains only cancer cells that are distributed uniformly in *Ω*(*t*), with density 1. The proliferation rate is proportional to the concentration of nutrients *c*, for example, oxygen. Then, after normalization, the function *c* satisfies a diffusion equation
2.1
where the right-hand side represents the consumption rate of the nutrients by the cells, with boundary condition
2.2
where *Γ*(*t*) is the boundary of *Ω*(*t*). Cells proliferate at a rate which depends on *c*; for instance,
where . Thus, cells grow if they receive nutrients at concentration above the critical number , but they shrink, or die, if the nutrients level is below . Cells that have died are removed from the tumour region, but the mathematical model does not account for this explicitly. Since the solution of (2.1) and (2.2) satisfies , it is natural to assume that , for otherwise the cells will die and the tumour will not materialize. Since the cell density is uniformly equal to 1, proliferation gives rise to a velocity field ** v** among the cells and, by conservation of mass,
Oxygen molecules are very small compared with cells, and their concentration is also very small. We therefore neglect their mass and, furthermore, assume that they move only by diffusion. We next assume that the tumour tissue is a porous medium so that, by Darcy's law,
where

*σ*is the pressure of the fluid-like cells in

*Ω*(

*t*). Hence 2.3 We are dealing here with a ‘solid’ tumour and so what keeps the tumour body together are the adhesion forces between cells at the boundary. This can be expressed by the boundary condition for the pressure: 2.4 where

*κ*is the mean curvature (

*κ*=1/

*R*(

*t*) if

*Ω*(

*t*) is a ball of radius

*R*(

*t*)). The boundary condition (2.4) represents surface tension, arising from cell-to-cell adhesion, which keeps the boundary of the tumour intact.

We assume that the boundary *Γ*(*t*) moves in accordance with the velocity ** v**, so that its component in the direction of the outward normal

**is 2.5 where ∂**

*n**σ*/∂

*n*is the derivative of

*σ*in the direction

**. Note that if**

*n**μ*=0 (no proliferation) then the problem is reduced to the Hele-Shaw problem for the pressure

*σ*. Tumour models were introduced as early as the 1970s by Greenspan [1,2], but the description of tumour evolution as a free boundary problem (2.1)–(2.5) was first proposed by Byrne & Chaplain [3].

### Theorem 2.1 ([4])

*There exists a unique radially symmetric stationary solution to the system* (*2.1*)–(*2.5*), *given by*
*and*
*and the radius R of the free boundary is the unique solution of the equation*

Note that if *μ* increases then the proliferation rate will increase as long as and if *γ* decreases then the cohesive force that keeps the tumour spherical will decrease. Thus, the parameter *μ*/*γ* may be viewed as representing the aggressiveness of a ‘well-fed’ tumour.

Hence we expect that the spherical tumour will be stable if *μ*/*γ* is small, but unstable if *μ*/*γ* is large. More precisely we have the following result.

### Theorem 2.2 ([5,6])

*For any integer* ℓ≥2, *there exists a branch of stationary solutions, bifurcating from a point μ*/*γ*=*M*_{ℓ}*, with free boundary*
2.6
*for any small* |*μ*/*γ*−*M*_{ℓ}|, *where Y* _{ℓ,0} *are the spherical harmonics of mode* (ℓ,0),
*and*
*where I*_{m}(*r*) *is the modified Bessel function of order m.*

Each bifurcation branch can be viewed as ‘fingers’ emanating from the spherical tumour; a larger number of such fingers indicates an increased risk of cancer metastasis. The radially symmetric stationary solution (RSSS) cannot be asymptotically stable at *μ*/*γ*=*M*_{2}. Indeed, any stationary solution with free boundary (2.6) and |*μ*/*γ*−*M*_{2}| positive and arbitrarily small lies arbitrarily close to the spherical solution but does not converge to it as . The loss of asymptotic stability of the RSSS, however, may occur already for smaller values of *μ*/*γ* as asserted in the following theorem.

### Theorem 2.3 ([7–9])

*There exists a positive number R*_{*} *such that the following hold*:

(i)

*if R*>*R*_{*}*then the RSSS is asymptotically stable as long as μ/γ*<*M*_{2};(ii)

*if R*>*R*_{*}*then the stationary solution*(*2.6*)*with*ℓ=2*is linearly asymptotically stable if μ/γ*−*M*_{2}*is positive and small, but linearly asymptotically unstable if μ/γ*−*M*_{2}*is negative and small;*(iii)

*if R*<*R*_{*}*then the asymptotic stability of the RSSS holds if μ/γ*<*M*_{*}*for some number M*_{*}<*M*_{2}*, but ‘Hopf bifurcation’ occurs at μ/γ*=*M*_{*}.

The ‘Hopf bifurcation’ is understood in the following sense: If we linearized the free boundary problem about the stationary solution at *μ*/*γ*=*M*_{*}, then any solution of the linearized problem approximates a periodic solution as . This means that the tumour is oscillating nearly periodically as, in fact, has been observed in several tumours [10].

The numbers *R*_{*} and *M*_{*} are given as solutions of some transcendental equations.

If we restrict the initial conditions on *c* and *Ω*(0) to be radially symmetric, then global existence and uniqueness of a radially symmetric solution of (2.1)–(2.5) can be established [1], and the question arises:

Is the stationary solution globally asymptotically stable when *R*>*R*_{*} and *μ*/*γ*<*M*_{2}, or when *R*<*R*_{*} and *μ*/*γ*<*M*_{*}?

A partial result is proved in [4].

We next consider the case where the tumour includes three types of cells: proliferating (*p*), quiescent (*q*) and dead cells (*a*). We assume uniform distribution of cells, so that
here again we neglect the concentration of nutrients (e.g. oxygen) which is relatively very small. Hence, owing to proliferation of *p* cells and removal of dead cells, there is a velocity field ** v**. Cells can change from

*p*to

*q*or vice versa, or from

*p*,

*q*to

*a*, depending on the available supply of nutrient

*c*. By conservation of mass, we then have 2.7 and 2.8 with

*F*and

*G*which account for the transition among the cells, and 2.9 where

*K*(

*c*)

*p*is the proliferation rate of the

*p*cells and

*K*

_{0}is the removal rate of dead cells. The function

*c*satisfies a diffusion equation 2.10 As before we assume the Darcy law, so that

**=−∇**

*v**σ*, and then 2.11 The model (2.7)–(2.11) was first introduced and investigated by Pettet

*et al.*[11]

### Theorem 2.4 ([12])

*The system* (*2.7*)–(*2.11*) *with boundary conditions* (*2.2*), (*2.4*), (*2.5*) *and smooth initial conditions for p, q and c has a unique smooth solution for some time interval* 0<*t*<*T*.

This system also has a unique radially symmetric stationary solution [13] and it is linearly asymptotically stable [14]. However, this solution is not explicitly known, and the problem of developing bifurcation theory for this solution has not been addressed.

The mathematical theory of a generic tumour has been extended further in several directions, but open problems abound.

(i) Tumour tissue is fluid-like so that Darcy's law is to be replaced by the Stokes equation [15–17].

(ii) Tumour evolves in two time scales: the cell cycle time (of the order of days) and the tumour growth time (of the order of months and years) [18–20].

In the two-time-scale models, proliferating cells that are in cell cycle phase *G*_{i}, *p*_{i}(*x*,*t*,*s*_{i}), satisfy a mass conservation law
but jump relations hold at the end-time when cells have to decide whether to enter the next phase, *G*_{i+1}, or remain quiescent for a while, or commit suicide (apoptosis) when irreparable damage has occurred, for instance in DNA replication. In this context, DNA mutations come into play, and thus multidimensional considerations enter: How do one or several mutations in cells affect the growth of the whole tissue? Such questions have been addressed in [20,21].

## 3. Wound healing

Chronic wounds represent a major public health problem in the USA and worldwide. A major complicating factor in wound healing is ischaemia, that is, insufficient supply of oxygen to the wound area, which is a result of damage to the vascular system around the wound. There are several mathematical models of wound healing which incorporate the effect of oxygen insufficiency [22–25]. Here we introduce a model, proposed in [25,26], which includes the dynamics of the wound boundary as a free boundary. We consider a cutaneous wound, that is, a dermal wound lying in {*z*<0}, where {*z*=0} is the surface of the skin. We denote the wound region by *W*(*t*) and the partially healed region by *Ω*(*t*)=*D* ∖ *W*(*t*) where *D* is a fixed domain. We take
and
Then
the common boundary of *Ω*(*t*) and *W*(*t*), is the free boundary.

The wound healing process begins when blood platelets and immune cells, called macrophages, first arrive in the surroundings of the wound and secrete signalling proteins that attract fibroblasts and blood vessels into the wound area. The fibroblasts secrete proteins, called collagens, in order to build up the damaged tissue. The density of the healing tissue in *Ω*(*t*), *ρ*(*x*,*y*,*z*,*t*), can vary from 0 to a maximal value *ρ*_{m}, and it exerts internal isotropic pressure *P*, which we model by
3.1
where *ρ*_{0}∈(0,*ρ*_{m}), *β* is a positive constant and *H*_{ϵ} is a smooth approximation to the Heaviside function. Oxygen, *w*, is a critical nutrient for healing, and we assume that it is supplied from the boundary of ∂_{0}*D*=∂*D*∩{*z*<0} of the partially healed tissue:
3.2
where the derivative ∂/∂*n* is in the outward normal direction, and *w*_{0} is the concentration of blood in the healthy vascular system. The parameter *α* expresses the ischaemic level: *α*=1 means total ischaemia, oxygen does not flow into the region *D*, whereas *α*=0 means that the vascular system at ∂_{0}*D* is in a healthy state.

The function *ρ* satisfies the mass conservation law
3.3
where
3.4
where *k*, *K*, λ, *ρ*_{m} are positive constants, and *f* is the density of fibroblasts.

Here ** v**=(

*v*

_{1},

*v*

_{2},

*v*

_{3}) is the internal velocity within

*Ω*(

*t*) which is a result of tissue growth and movement. The variables

*f*and

*w*satisfy a system of PDEs together with several more variables which are involved in the complex dynamics that models the healing process. We note that if

*α*=0 then from the PDE for oxygen it follows that

*w*≡0, so

*G*=0,

*ρ*will not increase, and the wound will not heal.

The free boundary moves according to the same law as in (2.5) or, in another notation, if *Γ*(*t*) is defined by an equation *φ*(*x*,*y*,*z*,*t*)=0, then
3.5

As in [26] we assume that the healing tissue is undergoing stress *σ*, *σ*=−*PI*+*τ*, where *P* is the isotropic pressure introduced in (3.1), *I* is the identity matrix and *τ* is the deviatoric stress. We also assume that the healing process is quasi-steady with negligible inertia so that the momentum equation becomes ∇⋅*σ*=0. Assuming that the tissue is viscoelastic (more precisely, upper convected Maxwell fluid with negligible relaxation time), the equation *τ*=*η*(∇** v**+∇

*v*^{T}) was derived in [26], where ∇

*v*^{T}is the transpose of the matrix ∇

**and**

*v**η*is a material constant. Hence, after normalizing to make

*η*=1, we get 3.6 in the healing tissue. Since there is no physical pressure originating from the interior of the wound, we also have 3.7 where

**is the normal. Consider first the very special geometry of radially symmetric ‘flat’ wound, as in [26].**

*n*Assuming that all the variables are radically symmetric, it can be deduced from (3.6) and (3.7) that the radial velocity *v*(*r*,*t*) satisfies the equation
3.8
and the boundary conditions
3.9

Numerical simulations of the model displayed in [26] show quantitative agreement with experimental results reported in [27]. We also have the following mathematical result.

### Theorem 3.1 ([28])

*The ‘flat’ wound problem has a unique smooth solution with free boundary r*=*R*(*t*), *and*

(i)

(ii)

*if α is near 1 then there exist R*_{*}>0*and t*_{*}>0*such that*

The last statement means that ischaemic wounds do not heal.

In the case where *α* is near 0 we expect the wound to heal, in the sense that , but this has not been proved and it remains an open problem, even in the case where *α*=0.

We now turn to the three-dimensional case. In this case, equation (3.6) for ** v**=(

*v*

_{1},

*v*

_{2},

*v*

_{3}) takes the form of a strongly elliptic system [29,30] 3.10 the boundary condition (3.7) on the free boundary takes the form 3.11 where

**=(**

*n**n*

_{1},

*n*

_{2},

*n*

_{3}) is the outward normal, and 3.12 Furthermore, on the remaining boundary of

*Ω*(

*t*), which lies on

*z*=0, we have 3.13 These conditions mean that the healing tissue is quiescent at the skin: both vertical velocity and acceleration vanish at

*z*=0.

We can then extend the system across *z*=0:*v*_{1} and *v*_{2} are extended by reflection and *v*_{3} is extended by anti-reflection. This allows us to consider a problem in which the free boundary does not meet the fixed boundary, and there will be no singularities at *Γ*(*t*)∩{*z*=0} provided we assume that
3.14

### Theorem 3.2 ([30])

*Under the condition* (*3.14*), *the three-dimensional wound problem has a unique smooth solution for some time interval* 0<*t*<*T.*

The proof uses the following approach. Introduce a family of surfaces
3.15
where λ varies in the region
3.16
and (*e*_{r}(λ),*e*_{θ}(λ),*e*_{φ}(λ)) is an orthonormal frame in the direction of increasing *r*,*θ*,*φ*. Given a function *h*(*t*,λ), we solve the PDE system of the wound problem in the domain with fixed boundary (3.15) and then use the free boundary condition (3.5) to define a new function . It was proved in [30] that the mapping has a unique fixed point if 0<*t*≤*T*, which yields the solution asserted in theorem 3.2.

The extension of the solution to all times remains a difficult open problem, so here we formulate a much simpler one.

Consider the cylindrically symmetric case, so that
and, in fact, take
Then the free boundary moves according to the equation
where , and is the radial velocity. Consider the elliptic system (3.8)–(3.11) with *P* a given function of (*r*,*t*). Find conditions on the function *Z*(*r*,0) and the pressure *P*(*r*,*z*) such that the wound begins to heal, that is, *Z*_{t}>0 at *t*=0, i.e.
Note that this is a purely elliptic problem with fixed boundary!

## 4. Biofilms

A biofilm is an aggregate of microorganisms in which cells stick to a surface and are typically embedded within a self-produced matrix of extracellular polymeric substance (EPS). The polymeric substance, which is often referred to as ‘slime’, is composed of proteins and polysaccharides. The matrix shields the microbes from external attack by chemicals and antibiotics. Biofilms are ubiquitous. They are found on the surface of stagnate water, on rocks, in the bottom of rivers, in the food chain and on living organisms. Macrobial biofilms are implicated in 80% of all infections in humans. They are found in urinary tract infections, in middle ear infection, in cutaneous wounds, in dental plaque, in contact lenses and in catheters. They block the diffusion of drugs into infected or injured areas. It is important to understand how biofilms evolve in order to devise ways that will wash away or destroy their protective matrix. There are many mathematical models of biofilms; we refer to [31,32] for comprehensive reviews of the work done up to 2010. The more complete models include two regions. In one region we have just fluid, while in the other region we have a mixture of the EPS with its embedded bacteria and the solvent. The fluid is often pure or contaminated water, or biofluid. The fluid is typically in motion with velocity *v*_{s}, while the EPS (with its embedded bacteria) is moving with a different velocity, *v*_{n}. The bacterial cells account for a small portion of the biofilm [33] but their distribution within the EPS is of great interest [34]. The pure fluid phase is modelled by a single Stokes equation for *v*_{s}, while in the biofilm region we have a coupled system of two Stokes equations, for *v*_{s} and *v*_{n}; the coupling occurs because the fractional volume *θ*_{n} of the EPS varies with time, while the Stokes equations include *θ*_{n} as a coefficient. The biofilm is attached to a fixed surface; here we consider a special case where the interface between the two phases does not intersect the fixed surface. We also set *v*_{s}=*v*_{1} in fluid, *v*_{s}=*v*_{3} in biofilm and *v*_{n}=*v*_{2}.

We denote by *Ω*(*t*) the region occupied by the biofilm. This region lies over a solid body *B*, and we assume that *Ω*(*t*) is surrounded by a fluid region *D*(*t*) with outer boundary *S*. A simple initial geometry could be, for instance:
4.1
The free boundary is the common boundary *Γ*(*t*)=∂*Ω*(*t*)∩∂*D*(*t*) between *Ω*(*t*) and *D*(*t*). The PDE system in the fluid region *D*(*t*) is given by
4.2
where *v*_{1} and *P*_{1} are, respectively, the velocity and pressure of the fluid in the domain *D*(*t*).

In the biofilm region *Ω*(*t*), we have a more complicated system of PDEs,
4.3
where, as mentioned above, *θ*_{n} is the volume fraction of the polymer, (1−*θ*_{n}) is the volume fraction of the fluid, *v*_{2} and *v*_{3} are the velocities of the polymer and fluid, respectively, and *P*_{2} is the pressure of the mixture.

The boundary conditions are
4.4
where *V*_{n} is the velocity of the boundary *Γ*(*t*) in the outward normal direction ** n**.

The initial condition for *θ*_{n} is

The functions *f*_{i}, *g*_{j} and *G* include cross-linking components of the biofilm matrix, chemical potentials and free energy [31,32]. The first equation in (4.3) should actually be split into two equations, one for the bacteria and another for the EPS, but this does not affect the proof of the following theorem.

### Theorem 4.1 ([35])

*Under some regularity and consistency conditions on the initial and boundary conditions and on the functions f*_{i}*, g*_{j} *and G, the biofilm system* (*4.2*)–(*4.4*) *has a unique smooth solution for some time interval* 0<*t*<*T*.

More precisely, the free boundary has the form
where *e*_{n}(λ) is the normal to the surface *Γ*(0) at *X*^{0}(λ) and
where *Λ* was defined in (3.16).

The proof begins in the same way as the proof of theorem 3.2. However, it is far more complicated because of the mixed boundary conditions among the three phases (with their velocities *v*_{1},*v*_{2} and *v*_{3}) on the free boundary.

It would be interesting to find sufficient conditions on the data under which the biofilm region begins to grow or shrink, for example, under strong external jet flow aimed to wash out the biofilm.

Another open problem is to consider the case where the initial free boundary intersects the surface ∂*B* on which the biofilm is based, and study the behaviour of the solution near the intersection points.

We note, in conclusion, that in the full mathematical model of a biofilm the inhomogeneous functions *f*_{i}, *g*_{j} and *G* are nonlinear functions of the velocities and of *θ*_{n}, but the proof of theorem 4.1 easily extends to this case; however, in most biofilm models *ϵ* is taken to be zero, in which case the proof of theorem 4.1 runs into difficulties.

## 5. Granulomas

Granulomas are small nodules which develop when a collection of immune cells attempt to wall off foreign substances that they are unable to eliminate. Such substances include infectious bacteria and parasites. Granulomas develop, for instance, when *Mycobacterium tuberculosis* is inhaled and endocytosed by alveolar macrophages. The bacteria may multiply within these macrophages and, as the infected macrophages die, a swarm of new bacteria are released and are again endocytosed by alveolar macrophages. As the bacterial attack continues, proinflammatory immune cells, namely, classically activated macrophages and T cells, arrive at the scene and begin to kill, or at least wall off, the bacteria. A granuloma thus formed may either grow or disappear over time, or else become stationary.

In terms of a mathematical model, a granuloma consists of a region *Ω*(*t*) with free boundary *Γ*(*t*) and variables such as macrophages, T cells and bacterial densities, which satisfy a system of PDEs with nonlinear coupling terms in *Ω*(*t*). The bacteria are found both inside macrophages as well as in the extracellular space. All the cells, including the bacteria, are moving with common velocity ** v** and are subject to diffusion. The evolution of the granuloma depends on the flux of immune cells into the granuloma. The biologically important question is whether a granuloma will grow, shrink or become stationary. All these three options can occur biologically.

Here we consider an extremely simple model with two variables: macrophages with density *M* and extracellular bacteria with density *B*.

By conservation of mass,
5.1
and
5.2
The parameter *μ*_{1} is the rate at which macrophages are killed by bacteria and, correspondingly, *μ*_{2} is the rate by which external bacteria are endocytosed by macrophages. The parameter *α* is the natural death rate of macrophages, while λ is the growth rate of external bacteria. Since bacteria are far smaller than macrophages, they diffuse faster; hence the parameter *δ* is positive. We take all these parameters to be constant, but this is an oversimplification. Indeed, the growth rate of external bacteria, λ, should account for those internal bacteria *B*_{i} which are released from infected macrophages upon their death. Similarly, the term *μ*_{1}*MB* does not distinguish between macrophages that are not yet infected by bacteria and those already infected. However, already the simplified model (5.1) and (5.2) raises interesting mathematical questions.

In this model (5.1) and (5.2), the velocity ** v** is still unknown. We make the assumption of uniform cell density:
5.3

Then, adding equations (4.1) and (4.2) we get
5.4
and if we assume Darcy's law
where *p* is the pressure within *Ω*(*t*), then the system can be written in the following form:
5.5
where
5.6

We assume that macrophages enter *Ω*(*t*) from the boundary *Γ*(*t*) at rate *β*, that is,
5.7
where ** n** is the outward normal, so that
5.8

We also assume that the granuloma is held together by surface tension induced by adhesive forces between adjacent macrophages at the boundary, so that, as in (2.4), 5.9 Finally, we impose the usual free boundary condition 5.10 and an initial condition 5.11

### Theorem 5.1 ([36])

*Under radially symmetric initial conditions, the system* (*5.4*)–(*5.11*) *has a unique radially symmetric solution for all t*>0, *with free boundary r*=*R*(*t*) *satisfying, for some positive constants c*_{0}*, C*_{0} *and C, the inequalities*
*furthermore, if* ∂*B*_{0}/∂*r*≤0 *then* ∂*B*(*r,t*)/∂*r*≤0 *for all* 0<*r*<*R*(*t*),*t*>0.

The inequality ∂*B*(*r*,*t*)/∂*r*≤0 means that, as new macrophages keep entering the granuloma across the boundary, the bacteria recede into the centre of the granuloma.

Under different conditions on the parameters of the system and on the initial conditions, it was shown in [36] that the granuloma may disappear, or may persist. So the question whether stationary granulomas exist, as is biologically observed, is interesting. This question was addressed in [37] for the case *δ*=0.

We write the stationary system in the form:
5.12
and
5.13
where the function *f*(*σ*) was introduced in (5.6), and
5.14
where *r*=*R* is the stationary free boundary point; at this point, the function *B*(*r*) must satisfy the boundary condition
5.15
with a given *β*. But instead of imposing condition (5.15), we shall impose the initial condition
5.16
and then *β* will depend on the choice of *σ*.

If *μ*>λ then the function *f*(*r*) has a unique positive zero, *σ*_{0} and

### Theorem 5.2 ([37])

*If μ*>λ *then there exists a maximal interval* (*σ*_{0}*,σ*_{1}) *with σ*_{1}≤1 *such that for any σ*∈(*σ*_{0}*,σ*_{1}) *there exists a unique stationary solution B*(*r*) *of* (*5.12*)–(*5.14*) *and* (*5.16*) *with v*(*r*)>0 *and* ∂*B*(*r*)/∂*r*<0, *and the range of the radii R*=*R*(*σ*) *covers the interval* 0<*R*<*R*_{0} *where R*_{0}=*j*_{0}/λ^{1/2}*, j*_{0} *being the smallest zero of the zeroth-order modified Bessel function I*_{0}(*r*).

The proof shows that the mappings *σ*→*R*(*σ*) and *σ*→*β* (in (5.15)) are continuous.

### Theorem 5.3 ([37])

*If μ*>λ *and σ*−*σ*_{0} *is sufficiently small then the stationary solution is linearly asymptotically stable with respect to any perturbation of the initial data by spherical harmonics Y*_{nm}(*θ,φ*) *with mode* (*n,m*), *n*≥2, *but unstable if n*=*m*=0.

We note that perturbation of mode *n*=*m*=0 changes the volume of the stationary granuloma, and the instability, in this case, is also suggested by biological considerations, since the bacteria flux rate *β* across the boundary will change (under condition (5.16)).

The extension of theorems 5.2 and 5.3 to the general case where *δ*>0 remains open.

Sarcoidosis is a disease in which granulomas are formed in the lungs and other organs, but its primary cause remains a mystery. Sarcoidosis of the lungs contributes to disability and increased mortality in patients. A mathematical model of sarcoidosis in the lungs was developed by Hao *et al*. [38]. The model includes macrophages, various types of T cells, and cytokines which either are produced by, or activate, these cells. The disease is assumed to begin with inflammation affecting macrophages. Assuming spherical symmetry, the free boundary *r*=*R*(*t*) depends (indirectly) on the level of inflammation *f*. It would be interesting to estimate the growth or shrinkage of *R*(*t*) in terms of the level *f* of inflammation.

## 6. Atherosclerosis

Atherosclerosis, the leading cause of death in the USA, is a disease in which a plaque builds up inside the arteries. As the plaque continues to grow, the shear force of the blood flow through the decreasing cross section of the lumen increases. This force may eventually cause rupture of the plaque, resulting in the formation of thrombus, and possibly heart attack. The process of plaque formation involves several cells and proteins within the plaque as well as influx of cholesterols and immune cells from the arterial blood into the plaque; there is also reverse cholesterol transport from the plaque into the artery aided by the ‘good’ cholesterols. A mathematical model of plaque evolution was developed and a ‘risk map’ was displayed in [39,40]. The risk map predicts the rate of plaque growth or shrinkage for any given level of ‘good’ and ‘bad’ cholesterols in the blood, and it includes risk factors such as smoking and high blood pressure. Here, we present a simple version which includes the following variables: The HDL removes the bad cholesterols from foam cells and the foam cells then return to become macrophages.

The following system in the plaque region *Ω*(*t*) is a simplified version of the model developed in [39,40]:
6.1
6.2
6.3
6.4

The first term on the right-hand side of (6.1) represents LDL ingested by macrophages at a rate limited by *k*_{1} (as ). Similarly, in equation (6.2) the clearing rate of foam cells by HDL is limited by *k*_{2}. Equations (6.3) and (6.4) include exchanges between macrophages and foam cells: macrophages that have ingested LDL become foam cells, and foam cells that were cleared from their LDL by HDL become macrophages again. Additionally, we included in equation (6.3) the term λ*ML*/(*δ*+*H*) which represents a sequence of processes involving LDL, oxidized LDL, chemoattractant proteins, T cells and IFN-*γ*, that collectively enhance activation of macrophages [40]; this enhancement is triggered by the bad cholesterols but inhibited by the good cholesterols [41]. Equations (6.1)–(6.4) also include degradation of LDL and HDL, and death of macrophages and foam cells.

We assume that the cells are distributed in the plaque with uniform density and take
6.5
The part of the boundary of the plaque in contact with the blood, *Γ*(*t*), is a free boundary, while the remaining part, ∂_{0}*Ω*(*t*), consists of the original surface of the arterial wall. The boundary conditions on *Γ*(*t*) are
6.6
where *L*_{0} and *H*_{0} are the cholesterol concentrations in the blood of an individual, and
6.7
Note that we have assumed here that in the blood the density of macrophages is *M*_{0}, and that there are no foam cells in the blood.

Consider the special case where the plaque is an infinite cylinder, so that we can focus on radially symmetric geometry with free boundary *r*=*R*(*t*); the plaque occupies a region
We impose a no-flux condition at *r*=*L* on all the variables. Adding equations (6.3) and (6.4) and setting ** v**=

*v*(

*r*,

*t*)

*e*_{r}, where

*v*is the radial velocity, we get 6.8 and 6.9 It was proved in [41] that an initial ‘

*ϵ*-thin’ plaque will dissolve as if 6.10 but will persist if 6.11 provided

*ϵ*is sufficiently small. Furthermore, for any small

*ϵ*>0 there exists a unique

*ϵ*-thin stationary plaque with

*L*

_{0}which is uniquely determined in the form

*L*

_{0}=

*μ*

_{1}(

*δ*+

*H*

_{0})/λ+

*K*

_{*}

*ϵ*+

*O*(

*ϵ*

^{2}) for some constant

*K*

_{*}, and it is linearly asymptotically stable or unstable depending on whether

*S*

_{*}>0 or

*S*

_{*}<0; both constants

*K*

_{*}and

*S*

_{*}are expressed in terms of the model parameters.

The above results suggest, in particular, that a person with cholesterol levels satisfying (6.10) is risk-free, whereas a person for whom (6.11) holds is not risk-free of atherosclerosis.

Let us denote the free boundary *r*=*R*(*t*) corresponding to the parameters *L*_{0} and *H*_{0} by *r*=*R*(*t*;*L*_{0},*H*_{0}). We conjecture that
i.e. the plaque grows as the bad cholesterol count (in the blood) increases and the good cholesterol count decreases.

## 7. Conclusion

In mathematical biology, the primary aim is to develop models that can predict known experimental results, and then use the models to suggest new hypotheses that can be tested experimentally or clinically. But it is also expected that at least some of the models will inspire new mathematics. In this paper, I focused on models which involve a portion of tissue that varies in space and time, i.e. I focused on free boundary problems. Furthermore, I put emphasis on new mathematics that is inspired by the biology. In the case of wound healing [26] or sarcoidosis [38], the models actually make predictions that quantitatively agree with *in vivo* data, and in atherosclerosis [39,40] the risk map is in agreement with the American Heart Association guidelines, and may suggest therapeutic strategies which are personalized [39]. On the other hand, the tumour models in §2 are too simple to be useful to biologists, but they offer a rich source of interesting mathematical problems. In the case of biofilms, a rigorous mathematical basis was established, but the usefulness of the model to biology and to mathematics will be tested when studying specific experimentally driven or medically driven situations.

## Competing interests

I declare I have no competing interests.

## Funding

I received no funding for this study.

## Footnotes

One contribution of 15 to a theme issue ‘Free boundary problems and related topics’.

- Accepted November 23, 2014.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.