FREE BOUNDARY PROBLEMS FOR SYSTEMS OF STOKES EQUATIONS

. Recent years have seen a dramatic increase in the number and variety of new mathematical models describing biological processes. Some of these models are formulated as free boundary problems for systems of PDEs. Relevant biological questions give rise to interesting mathematical questions regarding properties of the solutions. In this review we focus on models whose formulation includes Stokes equations. They arise in describing the evolution of tumors, both at the macroscopic and molecular levels, in wound healing of cutaneous wounds, and in bioﬁlms. We state recent results and formulate some open problems.

1. Tumor growth at the macroscopic level. We denote a tumor domain in R 3 by Ω(t), and its boundary by Γ(t); Γ(t) is a "free boundary," i.e., it is not prescribed in advance, and it needs to be solved together with a system of PDEs that are to be satisfied in Ω(t). We assume that there are three types of cells in the tumor: proliferating cells with (mass) density r(x, t), quiescent cells with density q(x, t), and dead cells with density n(x, t). We assume that proliferating cells change into quiescent cells at a rate K R (c), which depends on the concentration c(x, t) of the nutrients within the tumor. Similarly, quiescent cells become proliferating cells at a rate K Q (c), and they become necrotic (dead) at a rate K N (c). Proliferating cells have a proliferation (or growth) rate K B (c). Necrotic cells are removed from the tumor at a constant rate K O . Naturally, in the tumor, and we denote the corresponding velocity by v. Then, by conservation of mass, We next assume that all the cells are of the same volume and mass, and that the total density of the cells is uniform throughout the tumor. Then, after normalization, we have r + q + n = 1.
Summing the three preceding conservation laws, we deduce that div v = K B (c)r− K O n. This equation can be used to replace the conservation law for n. We now substitute n = 1 − r − q into the expression for div v and, together with the conservation laws for r and q, we obtain the system div v = h(c, r, q) in Ωt), t > 0, (1.1) The nutrient concentration c is depleted as it is consumed by the live cells. We assume that it satisfies a quasi-stationary diffusion equation where λ is a positive constant. As in [2,18,19,22] we assume that the tumor tissue is fluid-like as it develops, for example, in brain cancer, or in breast cancer which is confined to the mammary glands. The fluid velocity v = (v 1 , v 2 , v 3 ) and the fluid pressure p then satisfy the Newtonian constitutive law where σ ij is the stress tensor, p = − 1 3 σ kk is the isotropic stress, is the strain tensor, ∆ = e kk = div v is the dilatation, and ν is the viscosity coefficient. If there are no body forces then ∂σ ij ∂x j = 0, and we can rewrite this equation as the Stokes equation We now turn to the boundary conditions at the free boundary Γ(t). We assume that the tumor is held together by the forces of cell-to-cell adhesion with constant intensity γ. Introducing the stress tensor we then have where n is the outward unit normal and κ is the mean curvature (κ > 0 if Γ(t) is the surface of a convex body). We also assume the kinematic condition where V n is the velocity of the free boundary Γ(t) in the direction n. Finally, we take We note that no boundary conditions are needed for r and q, since (by (1.8)) the characteristic curves initiating at Γ(0) will remain on Γ(t) for all t > 0.
The system (1.6), (1.7) has six-dimensional kernel V 0 consisting of rigid motions We must therefore add six scalar constraints. We can write these constraints in the form where A(t), B(t) are prescribed functions. Finally we prescribe initial conditions: 12) and assume that As in [20,21], it is convenient to introduce the Lagrangian variable ξ by where u(ξ, t) = v(X( ξ, t), t),p( ξ, t) = p(X( ξ, t), t).
Theorem 1.2. [9]. If (1.14), (1.15) and (1.16) with T 0 = ∞ hold, and all the data are radially symmetric, then there exists a unique radially symmetric solution of (1.1)-(1.13) for all t > 0, with free boundary r = R(t), satisfying the regularity properties asserted in Theorem 1.1; in particular, It is interesting to determine whether there exist stationary radially symmetric solutions, and whether they are asymptotically stable in the following sense: A radially symmetric stationary solution is asymptotically stable if, for any initial data which belong to a small neighborhood of the stationary solution, there exists a unique global solution and it converges to the radially symmetric solution, or to a nearby translate of it.
Another interesting question is to determine whether there are branches of nonradially symmetric stationary solutions bifurcating from radially symmetric solutions. These will correspond to tumors which develop "fingers" and are thus about to metastasize.
The above questions have been studied, so far, primarily only in the special case where there are only proliferating cells, i.e., q ≡ n ≡ 0, r ≡ 1 and div v = µ(c −c), 0 <c < c. (1.17) The condition 0 <c < c ensures that the tumor will not disappear as t → ∞. In the sequel we normalize c by taking λ = 1 in (1.5).
Theorem 1.3. [12]. In the special case (1.17) there exists a unique radially symmetric stationary solution with free boundary r = R, and it is given by and R is the solution of the equation The constant p is determined by the condition In order to state bifurcation results for the solution (1.18)-(1.22), we introduce the quotient where I m (R) is the modified Bessel function given by , . Then there exists a branch of symmetry breaking stationary solutions bifurcating from the point µ γ = M n and having the form for all small | |, where Y n,0 (θ) is the spherical harmonic of mode (n, 0).
The idea of the proof is as follows: Consider a family of domains with boundary r = R + S(θ) and let ( v, p, c) be the solution of the free boundary problem, except for the kinematic free boundary condition. Setting we then wish to determine the parameter µ/γ such that and apply the Crandall-Rabinowitz theorem [1]. But this requires some estimates.
To derive these estimates we reduce the dimension of the problem (from 3 space variables to one) by representing the solution as a series of vector spherical harmonics. The next question is the asymptotic stability of the radially symmetric stationary solution.
Note that Theorem 1.5 establishes only linear asymptotic stability; the proof of asymptotic stability for the full nonlinear system will require additional analysis; cf. [10,11].
2. Multiscale tumor model. The cell cycle is divided into four phases: S (for synthesis), M (for mitosis), and gap phases G 1 and G 2 . During the S phase the DNA is replicated, that is, each chromosome is duplicated. During the mitosis phase, M , the nuclear membrane breaks down, sister chromatids are separated, new nuclear membranes are formed, and the cell divides into two daughter cells. S and M are separated by the two gap phases, G 1 and G 2 , during which the cell continues to grow. The cell cycle is controlled at two check points, R 1 and R 2 ; R 1 and R 2 are also called restriction points.
We consider a tumor model with two time scales: the time t during which the tumor evolves, and the time s during which the cell grows and divides; within the cell cycle we have the time s i during which the tumor cells progress in their i-phase of the cell cycle. The decisions on transition from one phase to another are being made at the two checkpoints R 1 and R 2 which are located, respectively, near the end of the G 1 phase and the G 2 phase. These decisions depend on the state of the cell as well as on signals the cell receives from the microenvironment; for example, signals of overpopulation or hypoxia: if there are already too many cells in the microenvironment, or if there is not sufficient supply of oxygen -do not continue with the cell cycle, go into a quiescent mode. When as a result of mutations the cells ignore such cues, there will be abnormal replication of cells and a tumor may develop.
At checkpoint R 1 the cell decides on one of three options: (i) to commit suicide (apoptosis) if it senses that it has been damaged beyond repair during the growth phase G 1 ; (ii) to go into a quiescent phase G 0 and stay there for a while, if the microenvironment is hypoxic or overpopulated with other cells; or (iii) to proceed to the S phase. At checkpoint R 2 the cell decides either to go into apoptosis if irreparable damage has occurred during the DNA replication, or to continue toward the M phase. A cell remains in G 0 for a period of time, by the end of which it proceeds to the S phase. We introduce the following notation: The variable x varies in the tumor region Ω(t), in R 3 , with boundary Γ(t). We denote by w(x, t) the oxygen concentration and by Q(x, t) the combined density of live cells in phases G 1 , S, G 2 , M . Due to cells proliferation and death, there is a velocity field v(x, t), which is assumed to be common to all the cells. Then, by conservation of mass, 2) where λ i (w) are growth rates, which depend on the oxygen level w. We assume that Here w * is the necrotic level of oxygen, λ 0 is the death rate of cells in quiescent mode, λ 4 is the clearance rate of dead cells, and µ 1 , µ 2 are the rates at which cells are damaged and go into apoptosis at the check points R 1 , R 2 , respectively. From what was explained above with regard to the check points R 1 and R 2 we infer the following jump relations: The function K(w, Q) in (2.6) represents the effect of signals given by the microenvironment: under hypoxic or overpopulation conditions the cell, at R 1 , is induced to go into quiescent mode; hence and, of course, K(w, Q) + µ 1 ≤ 1. SMAD is a gene that blocks the cell (at R 1 ) from proceeding from G 1 to S under hypoxic conditions, and APC is a gene that blocks it under overpopulation conditions.
It is convenient to introduce the notation The total density of each population of cells that are in the same phase of the cell cycle, is given by and then the density of the total population of live cells is It will be convenient to introduce the vector i=0 . In order to determine the velocity v, we integrate each of the equations in (2.1)-(2.3) over s i , 0 < s i < A i and sum up the resulting equations. Using the relations (2.5)-(2.8) we then get We assume that the total density of all the cells, live and dead, is constant and, for simplicity, take this constant to be 1, that is, (2.10) Conversely, one can easily show that (2.11) together with (2.9), (2.12) imply (2.10) for all t > 0 provided (2.10) is satisfied at t = 0.
Substituting (2.11) into (2.1)-(2.3), we obtain The oxygen concentration satisfies a diffusion equation where λ is a positive constant, w is the average oxygen concentration in a healthy tissue, and b is a transmission function of oxygen from the vasculature into the tissue. We finally assume that the velocity v satisfies the Stokes equation ( In the radially symmetric case with free boundary r = R(t), there exists a number β * > 0 such that The interpretation of the first case is: If both genes are mutated then the tumor increases beyond bound, so cancer occurs. It is also interesting to see what happens if only one gene, say SMAD, is mutated, so that K = K(Q), and K(Q) can be arbitrary control function.
Numerical calculation shows that one can get ρ 1 and ρ 2 to be close to R(0). We conclude that one mutation is not enough to produce uncontrolled proliferation.
3. Wound healing. We consider a 3-dimensional dermal wound W (t) lying in the half space {z < 0} and surrounded by partially healed tissue Ω(t): and Ω(t) = D\W (t) where D is a cylindrical domain Following [23] we assume that the partially healed tissue is a viscoelastic material where the stress σ has the form σ = −P I + τ , P being the isotropic pressure and τ the deviatoric stress, and that where η is the shear viscosity. If the healing is slow and the inertial forces are negligible, then the momentum equation in Ω(t) reduces to ∇ · σ = 0, so that, setting x 1 = x, x 2 = y, x 3 = z, we get [17] We denote by Γ(t) the common boundary ∂Ω(t) ∩ ∂W (t); Γ(t) is a free boundary. If we represent Γ(t) in the form φ(t, x, y, z) = 0 then the condition V n = v · n, where n is the normal to Γ(t), can be written in the form We assume that there are no external forces on the wound surface, so that where n = (ν 1 , ν 2 , ν 3 ). If the free boundary can be written as z = Z(t, x, y), then φ = Z(t, x, y) − z and The boundary conditions for v at the fixed boundaries are: Consider first the case where the isotropic pressure P is a given function P (t, x) where x = (x, y, z) and assume that Γ(0) has the form x = X(λ) where λ varies in region ∧ and X ∈ C 3+α (∧) (3.5) for some 0 < α < 1, and Theorem 3.1. [15]. If (3.5), (3.6) hold and Γ(0) intersects z = 0 orthogonally, then there exists a unique solution to the free boundary problem (3.1)-(3.4) for a small time interval 0 ≤ t ≤ T with Γ(t) in C 2+α in λ, and its t-derivative in C 1+α in λ. The first step in the proof is to reflect v 1 and v 2 across x 3 = 0 and anti-reflect v 3 across x 3 . In this way we obtain a free boundary which does not intersect the fixed boundary.
We next introduce a family of surfaces where e n is the normal to Γ(0) at X(λ), and h(t, λ) is a continuously differentiable function in (t, λ) with three continuous derivatives in λ, and with h(0, λ) ≡ 0. For any given h, we find a unique solution v of (3.1), (3.3) in the domains Ω(t) bounded by the surfaces X(t, λ). We then introduce a new functionh(t, λ) by means of the free boundary condition (3.2). The functionh satisfies a hyperbolic equation with coefficients which depend on h and v. The last step in the proof of Theorem 3.1 is showing that the mapping h →h is a contraction, so that it has a unique fixed point. In Theorem 3.1 it was assumed that Γ(0) intersects the plane x 3 = 0 orthogonally. If Γ(0) does not intersect x 3 = 0 orthogonally, then we expect that singularities will occur at Γ(t) ∩ {x 3 = 0}, and this case has not been investigated. We would like to find conditions which will ensure that healing takes place and that the wound shrinks. Consider for example the case of axial symmetry, so that ϕ = ψ(t, r, x 3 ) where r = (x 2 1 + x 2 2 ) 1/2 . We can represent Γ(t) in the form x 3 = Z(t, r) and the free boundary condition by Then the wound will begin to shrink if Z t > 0 at t = 0, i.e., if Problem. Find conditions on Γ(0) and P (0, r, x 3 ) which ensure that (3.7) is satisfied.
Note that this open problem is purely an elliptic one, with no free boundary. So far we assumed that P is a given function. But, actually, P is a monotone increasing function of ρ, and ρ evolves according to the law of mass conservation, Here w is the concentration of oxygen, and f is the concentration of cells called fibroblasts. Furthermore, the variables ρ, w, f and v are involved in a larger system of reaction-diffusion equations which includes several cell types and several growth factors. The existence and uniqueness of a solution for the complete system was proved in [15] for a small time interval. Global existence remains an open problem already for the special case considered in Theorem 3.1.

Biofilm.
A biofilm is an aggregate of microorganisms in which cells stick to a surface and are typically embedded within 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 attacks 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, 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 biofilm; we refer to [18] [22] 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 [19] but their distribution within the EPS is of great interest [2]. The pure fluid phase is modeled 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 equation 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) a 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: B is a sphere of radius R 0 , The free boundary is the common boundary Γ(t) = ∂Ω(t) ∩ ∂D(t) between Ω(t) and D(t). The PDE system in D(t) is given by where v 1 and P 1 are the velocity and pressure of the fluid in the domain D(t).
In the Ω(t) phase, we have a more complicated system of PDEs, where θ 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, and P 2 is the pressure of the mixture in the domain Ω(t).
The boundary conditions are where V Γ(t) is the velocity of the boundary Γ(t) in the normal direction n. The initial condition for θ n is θ n|t=0 = θ 0 in Ω(0).
This compatibility condition, on the free boundary x = X(λ, t) at t = 0, is the condition that lim t→0 d dt ∂θ n ∂n (X(λ, t), t) = 0, where ∂θn ∂t is computed from the differential equation for θ n and all the derivatives D j x θ n are replaced, as t → 0, by D j x θ 0 . Similarly, one defines the compatibility condition on ∂G × {t = 0}.
Theorem 4.1. [16]. Under the conditions (4.5)-(4.8), there exists a unique solution to the free boundary problem (4.2)-(4.4) for a small time interval 0 ≤ t ≤ T such that the free boundary has the form x = X 0 (λ) + h(t, λ) e n (λ) where e n (λ) is the normal to the Γ(0) at X 0 (λ) and h t , D λ h t , D 2 λ h t , D 3 λ h are in L ∞ ([0, T ] × ∧). The proof of Theorem 4.1 follows the general procedure as the proof of Theorem 3.1. However, it is far more complicated due to two factors: the coupling between the two phases at the free boundary (i.e., the first three equations in the system (4.4)) and the diffusion equation for θ n in the system (4.3) which is coupled to the Stokes equations for v 2 and v 3 .
It would be interesting to find conditions on the data under which the biofilm begins initially to grow or to decrease, for example in the case of the simple geometry (4.1). It would also be interesting to explore the case when the free boundary initially intersect the fixed boundary of G: what is the local behavior of the solution near the intersection of Γ(0) and ∂G?
We note in conclusion that the full mathematical model of biofilm is somewhat more complicated than the one formulated above: the inhomogeneous terms are nonlinear functions of the velocities and θ n , and in (4.3) is often taken to be equal to zero.