Free boundary problems arising in biology

The present paper describes the general structure of free boundary problems for systems of PDEs modeling biological processes. It then proceeds to review two recent examples of the evolution of a plaque in the artery, and of a granuloma in the lung. Simplified versions of these models are formulated, and rigorous mathematical results and open questions are stated.

Introduction. Recent years have seen a dramatic increase in the number and variety of mathematical models describing biological processes. Such models aim to produce, by simulations, results that are in agreement with experimental data, and develop new hypotheses that could be tested experimentally. Some of the mathematical models are formulated in terms of systems of partial differential equations (PDEs) and give rise to new and interesting mathematical questions. These include not only existence, uniqueness and regularity of solutions, but also proving biologically-motivated properties of the solutions. This is particularly the case when the domain where the system of PDEs to be solved varies with time and its boundary is not known in advance. Such free boundary problems are well known in elasticity, in phase transition, and in fluid dynamics. In recent years new free boundary problems (FBPs) appeared in mathematical models of biological processes, for example in tumor growth and in healing of dermal wounds.
The present paper has two aims. The first one is to describe the structure of FBPs that arise in biology, a structure that encompasses many specific models. The second aim is to review several very recent models, explain how they address biological questions and, at the same time, how they inspired new rigorous mathematical results and new mathematical open problems.

Differential equations for cells and cytokines.
We denote different species of cells by X 1 , X 2 , . . . , X m , and different species of proteins and molecules (such as oxygen) by Y 1 , Y 2 , . . . , Y n . For simplicity we shall denote the densities of X i by the same symbol X i , and similarly by Y j the concentration of Y j . Some of the cells belong to the immune system; for example, macrophages which 'eat' (endocytose) foreign bacteria, T cells that kill parasites, and dendritic cells which recognize invading micro organisms and alert the T cells. Other cells such as muscle cells, endothelial cells and epithelial cells have specific tasks. Endothelial cells line up the 194 AVNER FRIEDMAN inner surface of blood vessels' walls, while epithelial cells cover the inner surface of cavities of all the body organs. Proteins produced by immune cells are called cytokine. They are released by immune cells in order to stimulate or inhibit the activity of other cells. Fibroblasts are cells that produce proteins, primarily collagens, in order to build and maintain the tissue structure where cells live; this tissue is called the extracellular matrix (ECM) Cells grow, divide and die, and they also move within the ECM. Cancer cells grow abnormally fast, and they move to occupy new territories. Their movement is enabled by metaloproteinase (MMP), a protein that degrades the ECM. The healing of a dermal wound involves a movement of macrophages from the healthy tissue into the wound in order to clear the infection. This movement is initiated by platelets derived growth factor (PDGF), a cytokine produced by the macrophages. These are just two examples of how cells and proteins interact.
In developing a mathematical model of a biological process, we first need to be clear on what is the biological question that we want to address. We then proceed to introduce the biological species that must be included in addressing this question, and develop a system of equations that represent the biological process. The next step is to use mathematical theory and simulations to derive predictions of the model. The final step is to verify that these predictions agree with known experimental results and, at the same time, can be used to address the biological question. Such an approach was adopted in two recent textbooks in mathematical biology, one for undergraduates, using ODEs [1], and another for graduate students using also PDEs [9].
The biological processes take place in a domain Ω(t), and the structure of the tissue of Ω(t) must be taken into consideration. This involves, in particular, the ECM density, ρ. We assume that the combined densities of the cells and ρ is approximately constant at each point x of Ω(t), and take (1.1) Because of varying birth and death rates of various cells, which also affect ρ, the assumption (1.1) implies that cells may push against each other or move away from each other, thus giving rise to a velocity field u by which all cells are advected. But some cells may be subjected to additional velocity, by chemotaxis. Specialized cytokines, called chemokines, attract specific cells in the direction of their concentration gradient. For example, when infection occurs in an organ, epithelial cells produce major chemoattractant protein (MCP-1) which attracts macrophages from the blood into the infected organ. The chemoattractive velocity of cells X i by a chemokine Y j (j = j(i)) is given by χ j ∇Y j where χ j represents the force of attraction.
We can now write the general structure of the PDE system for the cells in Ω(t) as follows: where χ j(i) > 0 if X i is attracted chemotactically to a chemokine Y j(i) , and χ j(i) = 0 otherwise. Here D Xi is the diffusion (or dispersion) coefficient for species X i , and represents the various interactions among cells and cytokines that affect the density X i . The equations for proteins and other molecules are simpler: Indeed, since proteins are much smaller than cells, their diffusion coefficients are much larger, and their advection ∇ · (Y j u) is negligible compared to their diffusion term, so it may be neglected. The equation for the ECM density, ρ, is given by the conservation law:  (1.2). We next have to derive an equation for u using the fact that div u = H. There are several ways of doing it.
(i) We assume that the tissue has the structure of a porous medium, so that, by Darcy's law, u = −∇p where p is the internal pressure; hence (1.6) (ii) We assume that the tissue is fluid-like, and use Stokes equation (iii) Another way of deriving an equation for u without using (1.5) is the following. We assume that the tissue is viscoelastic with negligible inertial force and in quasisteady state. Then [8,16] where P , the isotropic pressure, is a given function of ρ. In terms of the coordinates (u 1 , u 2 , u 3 ) of u, the last equation takes the form of a strongly elliptic system, the free boundary. The domain Ω(t) varies within a tissue G with fixed boundary ∂G. The boundary ∂Ω(t) of Ω(t) may possibly intersect ∂G, and this situation typically results in singularities at the boundary points of their common portions. We shall consider here only the case where ∂Ω(t) lies within G.
The standard kinematic condition then says that the free boundary moves with the internal velocity u. Thus, V n = u · n (1.9) where V n is the velocity of ∂Ω(t) in the direction of the outward normal n to ∂Ω(t).
1.4. Boundary conditions. The boundary conditions for X i , Y j depend on the specific biological process. We give two examples.
1.4.1. Cancer. Tumor cells produce a cytokine, vascular endothelial growth factor (VEGF), which induces endothelial cells to move into the tumor, bringing along with them new blood capillaries that provide oxygen to the fast proliferating cells. Setting V = VEGF, X 1 = endothelial cells (E), Y 1 = oxygen(w) we then have the boundary conditions: where K V , K E are constants and E 0 , w 0 are the concentrations of E, w in healthy state. We may take no-flux boundary conditions for the remaining X i , Y j . Note that no boundary condition is needed for ρ, since the free boundary is a characteristic surface for the hyperbolic equation (1.4).
With regard p and u, in case (1.5) a natural boundary condition is where κ is the mean curvature of the free boundary (κ = 1 R if the domain is a ball of radius R) and γ is a positive coefficient representing the cell-to-cell adhesion at the tumor boundary. In the case of (1.6) the boundary condition is S n = −γκ n where S is the stress tensor and I is the identity matrix. Free boundary problems of cancer models have been rigorously studied for the last 20 years; we refer to review papers [2,3] and an additional recent work [11].
1.4.2. Dermal wounds. Figure 1 is a schematic of a 3-dimensional dermal wound W (t) surrounded by a partially healed tissue Ω(t) which is assumed to be a viscoelastic material as in (1.8).
The plane x 3 = 0 is the thin epidermal layer exposed to air, and the boundary Γ (in {x 3 < 0}) borders the healthy tissue. We assume that u = 0 on Γ. On ∂Ω(t) ∩ {x 3 = 0} it is natural to assume that We can then extend u 1 , u 2 across x 3 = 0 by reflection, and u 3 by anti-reflection. Similarly we extend the PDE system into the domain We assume that ∂Ω(0) intersects {x 3 = 0} orthogonally, in order to ensure that the free boundary ∂ Ω(t) is initially smooth around x 3 = 0. In the extended FBP the free boundary does not meet the fixed boundary Γ and its reflection, and (1.8) is thus assumed to hold throughout ∂ Ω(t). The no-stress boundary condition associated with (1.8) is where n = (n 1 , n 2 , n 3 ). We may assume no-flux condition on the free boundary for all the variables except for the cytokine PDGF which appears in the initial stage of healing process. Existence and uniqueness of a smooth solution in 3-d was proved in [8] for a small time interval. More general results were proved in the 2-d case of radially symmetric flat wounds, lying in x 3 = 0, with In particular, if the flux of oxygen is seriously limited because of vascular damage at r = R 0 , then the wound will not heal, that is, there exists a time t * such that R(t) ↓ R * for t ↑ t * and R(t) = R(t * ) > 0 for t > t * [7].
2. Atherosclerosis: The risk of high cholesterol. Atherosclerosis 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 types of cells and proteins within the plaque as well as influx of cholesterols and immune cells from the arterial blood into the plaque. We distinguish between low density cholesterol (LDL) and high density cholesterol (HDL). The blood contains both species of cholesterols. The process of plaque development begins with a lesion in the endothelial layer of an arterial wall, which enables cholesterols to enter the intima layer of the wall and become oxidized by radicals. The oxidized LDL alerts macrophages from the blood to come in and 'eat' the oxidized LDL; by doing so the macrophages eventually become 'obese' cells, called foam cells. On the other hand, the oxidized HDL are able to ligate to receptors on the foam cells and initiate a process that transports the cholesterol inside the foam cells out of the plaque and into the blood for recycling or elimination. Thus, LDL are 'bad' cholesterols while HDL are 'good' cholesterols. The question is what is the risk of plaque growth for a person whose cholesterol blood levels are L 0 for LDL and H 0 for HDL.
A comprehensive model was developed in [4,13]. The model includes also smooth muscle cells which are attracted from the media layer of the arterial wall into the intima and T cells. The model simulations produced a risk map that predicts, for any pair (L 0 , H 0 ), the percentage of growth or shrinkage of a small plaque over a period of time. The simulations in [4,13] are based on 2-d geometries where the profile of the plaque is shown either in a cross section of the artery (Fig. 2(a)) or in a planar cross section along the artery (Fig. 2(b)). Biologically inspired mathematical questions are: (i) do there exist steady plaques, and are they stable; (ii) under what relations between L 0 and H 0 will a small plaque grow, or shrink.
These questions have been addressed in [5] for a simplified model based on the diagram in Fig. 3. For simplicity we identify the 'clean' macrophages M 2 , derived from foam cells, with the M 1 macrophages that came from the blood to 'eat' the oxidized LDL.
Then, based on Fig. 3 we have the following system of PDEs: More precisely, we have identified here the LDL and HDL with the oxidized LDL and oxidized HDL, respectively. We also assumed that H, by consuming some of the radical available to L, acts to reduce the growth of M (which arrive into the plaque when alerted by L), and this is phenominally represented by the term λM L/(δ+H).
We assume that all the variables in Eqs. (2.1)-(2.4) are radially symmetric, and denote the plaque by Ω(t) = {R(t) < r < R 0 }. We assume that on the free boundary r = R(t), Thus, the fluxes of L, H and M from the blood into the plaque depend on their concentrations in the blood, while foam cells cleared out into the blood.
The following results were proved in [5]: (i) If L 0 < aH 0 + b then any small plaque will disappear as t → ∞ (that is, R(t) → R 0 as t → ∞); (ii) If L 0 > aH 0 + b then any initial plaque will persist for all t > 0; (iii) There exist a unique -thin plaque for any sufficiently small , with The parameters a, b and γ are given explicitly in terms of the parameters of the system (2.1)-(2.6). It was also proved in [5] that there is a parameter S 1 , depending explicitly on the parameters of the system (2.1)-(2.6) such that: (iv) The steady solution in (iii) is linearly asymptotically stable if S 1 > 0 and linearly asymptotically unstable if S 1 < 0.
3. Dynamics of granulomas. A granuloma is a small mass of tissue (nodule) consisting of immune cells, including macrophages, that aggregate in response to infection and attempt to block the spread of the infection. The infection may result from an autoimmune disease, or from invading pathogens such as Micobacteria tuberculosis (Mtb). The boundary of the granuloma is a free boundary.
Mathematical models of granulomas in sarcoidosis and tuberculosis (TB) have recently been developed [14,15]. Here we focus on tuberculosis. Approximately 1/3 of the world's population have Mtb in their body, most commonly in the lung. However, the disease is only in latent form, and it remains inactive as long as the immune system is not severely compromised; the Mtb pathogens are confined within granulomas.
A comprehensive model of a TB granuloma in the lung was developed in [15]. The model was used to explore the efficacy of experimental drugs for individuals with varying degrees of strength of their immune response. The model includes several types of macrophages and T cells, dendritic cells, within host and extracellular Mtb, and various cytokines released by the cells. The system consists of 22 PDEs in a domain Ω(t) with free boundary. The boundary conditions include continuous flux of macrophages and T cells migrating from the lymph nodes into the granuloma.
For mathematical analysis, we will consider here a highly simplified system con- Since the granuloma Ω(t) consists only of M and B, we assume that and there is a velocity field u by which both M and B are advected. We model the dynamics of M and B by the following equations: The encounter between M and B results in a loss for both of them, at rates µ 1 and µ 2 , α is the death rate of M , and λ is the growth rate of B. We assume that macrophages migrate from the lymph nodes into the granuloma, while pathogens are removed across the boundary, so that the system (3.5)-(3.9) becomes a free boundary problem for B with free boundary moving with velocity which depends nonlinearly and nonlocally on B.
In [12] it was proved that this problem has a unique global solution for any initial condition B(r, 0) = B 0 (r), 0 ≤ B 0 (r) ≤ 1.
(3.10) Various monotone properties and asymptotic behavior of the solution were also established in [12].
A biologically important question is whether there exists steady state granulomas.
In [10] it was proved that if D B = D M then there exists steady state granulomas for any radius 0 < R < j 0 /λ 1/2 (3.11) where j 0 is the smallest zero of the zeroth-order Bessel function. We next extend the model to the non-radially symmetric case. We assume that Darcy's law holds, so that u = −∇p, and then in Ω(t), (3.12) and take p = γκ on ∂Ω(t).
(3.13) It is natural to ask whether with this extended model, the radially symmetric steady states asserted in (3.11) are asymptotically stable, or linearly stable. It was proved in [10] that if D B = D M then the steady state granulomas with sufficiently small radius R are linearly asymptotic stable with respect to any (3.14) perturbation with spherical harmonics of mode ≥ 2, but they are unstable for zero mode. Since bacteria are much smaller than cells, D B is definitely larger than D M . It would therefore be very interesting to extend the results (3.11), (3.14) to the case where D B > D M . 4. Conclusion. In this paper we considered biological processes which take place in a portion Ω(t) of a tissue of an organ. The location of Ω(t) is unknown in advance, and it depends on the biological processes taking place within Ω(t). Mathematical models of such processes can be formulated as free boundary problems for systems of PDEs. We described the general structure of the PDE systems and of the velocity of the free boundary, that is, the boundary of Ω(t). We also gave examples of boundary conditions prescribed at the free boundary.
The mathematical models aim to provide answers to biological questions. But at the same time the models may raise mathematical questions that inspire new mathematics. We gave two such recent examples. The first example deals with atherosclerosis, where a plaque develops in the artery, which may cause a heart attack, and the second example deals with the dynamics of a granuloma, an aggregation of immune cells that aim to block a pathogen. In both cases we described, for simplified models, new mathematical results motivated by biological questions. These examples illustrate the opportunities for mathematical research inspired by biology.