On interfaces between cell populations with different mobilities

. Partial diﬀerential equations describing the dynamics of cell population densities from a ﬂuid mechanical perspective can model the growth of avascular tumours. In this framework, we consider a system of equations that describes the interaction between a population of dividing cells and a population of non-dividing cells. The two cell populations are characterised by diﬀerent mobilities. We present the results of numerical simulations displaying two-dimensional spherical waves with sharp interfaces between dividing and non-dividing cells. Furthermore, we numerically observe how diﬀerent ratios between the mobilities change the morphology of the interfaces, and lead to the emergence of ﬁnger-like patterns of invasion above a threshold. Motivated by these simulations, we study the existence of one-dimensional travelling wave solutions.

1. Introduction.Partial differential equations (PDE) describing the dynamics of cell population densities from a fluid mechanical perspective can model the growth of avascular tumours [2,20,23].These models rely on the observation that, when in mechanical contact with other cells, proliferating cells exert a pressure on their neighbours.This pressure results in a mechanical form of cell motion (i.e.neighbouring cells are pushed and forced to move) and, above a certain threshold, it can also induce quiescence due to competition for space [5,12,22].
Such PDEs can be used to investigate possible mechanisms that underlie cancer invasion, such as: the effects of local pressure gradients which are generated by proliferation and competition for space and resources; conditions for preserving the initial radial symmetry of the tumour, and conditions for developing irregular shapes and invading the surrounding tissue [1,3,8,9,6].Moreover, it is now known that models of this type are equivalent, in the incompressible limit, to models of tumour growth which are formulated as free boundary problems [18,19].
In this framework, we consider a PDE model describing the time dynamics of a population of dividing cells and a population of non-dividing cells, coupled by pressure forces.The two populations are characterised, respectively, by the local density of dividing cells m and the local density of non-dividing cells n, governed by the following system of equations: The second terms on the left hand sides model the tendency of cells to move down pressure gradients, and rely on the definition of the cell velocity fields through Darcy's law [7,14].The parameters µ > 0 and ν > 0 stand for the mobility (i.e. the quotient of permeability and viscosity) of dividing cells and non-dividing cells, respectively.In the compressible description, the pressure p is given by a constitutive relation whose desirable properties are discussed, for instance, in [10].For simplicity, we follow the lines of [8] and [19], and we use the constitutive relation where the parameter γ ≥ 1 controls the stiffness of the pressure law.Finally, the net growth rate G satisfies Assumptions (3) mean that competition for space decreases the cell division rate according to the local pressure.The parameter P M > 0 models the threshold pressure above which dividing cells are entering a quiescent state (i.e. the so-called homeostatic pressure) [8,22].
It is worth noting that, if we use m for the local density of cancer cells within a solid tumour and n for the concentration of cells in homeostatic equilibrium that surround the solid tumour, the case µ > ν is the most biologically relevant one.In fact, the loss of junctional contacts, such as the loss of E-cadherin mediated adhesion induced by the nuclear up-regulation of soluble β-catenin [21], endows cancer cells with a much higher mobility compared to the cells of the surrounding tissues.
Also, some biological scenarios are not captured by the system of equations under study.For instance, there are other PDE models for the dynamics of dividing and non-dividing cells which allow cells to pass from one state to the other, and include a death term for non-dividing cells as well.We refer the interested reader to [13,25], and references therein.
Our work follows earlier papers that are devoted to study the existence of travelling wave solutions with composite shapes and discontinuities for cell-density models of avascular tumour growth [18,19,26].These papers consider one single nonlinear PDE, while we deal with a nonlinear strongly-coupled system of PDEs.This represents a novel fundamental difficulty, and the current literature still lacks a systematic characterisation of such systems of equations, with the exception of semilinear cooperative systems.In this respect, it has been speculated that the system (1) could exhibit various behaviours depending on the values of the parameters µ and ν [12].A first observation is that segregation can occur: dividing and non-dividing cells remain well separated through a sharp interface.Phenomena of this type have already been analyzed for related parabolic systems [11,17].Furthermore, twodimensional spherical waves could be stable in the case when µ < ν, while they could become unstable for µ > ν.The rationale behind this intuition is that, since mobility is inversely proportional to viscosity, the case µ > ν corresponds to the situation where one "less viscous fluid" (dividing cells) expand into a "more viscous fluid" (non-dividing cells).This might trigger the emergence of patterns that would be reminiscent of the viscous fingering resulting from Saffman-Taylor-like instabilities [24], i.e. those instabilities that can be observed when a less viscous fluid is injected into a more viscous fluid in a radial Hele-Shaw cell.
The aforementioned intuitions are supported by our numerical simulations which we report in Section 2. They display sharp interfaces and show that two-dimensional spherical waves are stable for µ < ν and can be unstable for µ > ν.Motivated by these results, our goal here is to investigate the profile of one-dimensional travelling waves for the system (1) by using numerical simulations and qualitative analyses.
One-dimensional travelling waves may provide a tractable way to explain both segregation and the profiles of two-dimensional spherical waves through, for instance, transversal instability [15].They are solutions of the form where the constant σ > 0 represents the travelling wave velocity.Substituting these travelling wave ansatz into equations (1) gives An observation we will use throughout our analysis is that, since p vanishes at ±∞, equation ( 5) gives In Section 3, we construct travelling wave solutions such that, in (6), we have n ∞ = 0 for all values of µ and ν.These solutions have compact support in n, and they display a remarkable qualitative difference between the cases µ < ν (when they exhibit a segregation effect) and µ > ν (when stable mixing occurs).Notice that any two functions m and n such that K γ (n + m) γ = P M are solutions in the case where σ = 0.In Section 4, we extend the construction to the case where n ∞ > 0. In this setting, travelling wave solutions still exist but they are numerically unstable objects, and a numerical constraint has to be imposed to capture them.
2. Numerical observations.Based on numerical simulations, we first investigate the behaviour of two-dimensional spherical waves that, as mentioned before, can be unstable for µ > ν.To this end, we numerically solve the mathematical problem defined by completing (1) with zero Neumann boundary conditions and the following initial conditions, which mimic a biological scenario where dividing cells expand in an embedding medium made of non-dividing cells: m(x, y, t = 0) := a m e −bm (x 2 +y 2 ) and n(x, y, t = 0) := a n e −bn (x 2 +y 2 ) , with Moreover, in view of the considerations drawn in [12,26], and to satisfy assumptions (3), we use the following definition for the net growth rate G: Finally, in line with the ideas presented in [18,19], we choose γ = 30 (i.e.we focus on the incompressible limit γ → ∞).
Numerical computations are performed in Matlab.We select a uniform discretisation consisting of 450 2 points on the square [−L, L] × [−L, L], with L = 45, as the spatial domain.The method for calculating numerical solutions is based on a time splitting scheme between the conservative parts and the reaction term.A finite volume method is used for solving the conservative parts.Convection terms are approximated through an upwind scheme, and the cell edge states are calculated by means of a high order extrapolation procedure [16].
The numerical results are summarised in Fig. 1 and Fig. 2.These results highlight how variations in the parameters µ and ν can lead to major changes in the morphology of the interface between dividing cells and non-dividing cells.Both figures show that the additional mass generated by cell divisions induces a pressure gradient that moves the dividing-cell population against non-dividing cells.On the one hand, when µ < ν, we observe the emergence of a spherical wave of dividing cells pushing the surrounding non-dividing cells (see left panel in Fig. 2), and the creation of an invasive front made of non-dividing cells that are induced to move by the expansion of dividing cells (see right panel in Fig. 1).On the other hand, when µ > ν, the non-dividing cells cannot move away with the same speed as the growth front.This triggers the formation of dendrites made of dividing cells that protrude through and dislocate the surrounding non-dividing cells.As a result, finger-like patterns of invasion emerge (see Fig. 2).Theorem 3.1.For all values of µ, ν and for M > 0 given, there are a speed σ and a width r > 0 such that a solution of (4)-( 5), ( 9) exists, with m non-increasing, discontinuous at x = 0, and such that n satisfies r 0 n(x)dx = M > 0 (given mass).
The pressure has a kink at x = 0 with sgn(p (0)) = sgn([p ]) = sgn(ν − µ).The construction gives a limit as γ → ∞ (Hele-Shaw or incompressible limit), and the speed satisfies The numerical solution obtained for µ < ν is shown in Fig. 4. If µ > ν, this travelling wave solution is unstable, as it can be seen very intuitively by using the proof below.Take a small perturbation of n near x = 0, for which n > 0 in region II.Such perturbation will move with a velocity close to −νp (0 + ) < −νp (0 − ) (see condition (16)).This means that the perturbation will further separate from the main core rather than join it.The numerical solution obtained for µ > ν is presented in Fig. 5.This figure displays a transient regime after which n is left behind and m propagates alone (see also Supplementary Movie S1) with a profile analogous to that observed in the case where n = 0.Such a profile resembles that of travelling wave solutions to combustion models with non-linear diffusion which has been previously presented by Berestycki et al. [4].
Proof.As depicted in Fig. 3, we work on two regions: region I that coincides with [0, r] = Supp(n), and region II that is defined as (−∞, 0].We use a shooting argument for the parameter σ that we fix at the beginning. Step 1.Since we are considering the case n ∞ = 0, (6) implies (σ + νp )n = 0.In region I, as n > 0, this leads to and therefore, The mass conservation gives the value of r by Step 2. Adding equations ( 4) and ( 5), we obtain Moreover, in region II, since n ≡ 0, we can rewrite equation (4) as and the solution has to satisfy and boundary condition ( 11) at x = 0.By the maximum principle p cannot have a local minimum therefore p (x) < 0 for x < 0. So the solution p is larger than p(0) on (−∞, 0) and the equation ( 14) is non-degenerate.Therefore the solution p is continuous on (−∞, 0) and there is a point x 0 < 0 such that p (x 0 ) > −∞.
Multiplying equation ( 13) by p and integrating over a bounded interval (x 0 , r) gives since m and p are bounded.Therefore, p ∈ L 2 loc (x 0 , r) and so p is continuous; as a consequence, m + n is continuous as well.This also means that m(0 − ) = n(0 + ).
Step 3. Integrating equation ( 13), we obtain Using (10) and step 2, and by the continuity of p and m + n, we deduce the jump condition at 0 by Step 4.Here we want to show monotonicity of p in σ.To this end, let us define F (p) := pG(p) and q := ∂p ∂σ .Deriving equation ( 14) w.r.t.σ we obtain with q(−∞) = 0 and q (0) = − 1 µ .
Deriving equation ( 14 We rewrite equations ( 17) and ( 18) as where the coefficients a i are implicitly defined.So we view the equations ( 17) and (18) as linear in q and p , respectively.Now, the function 1 σ p satisfies the boundary condition of equation (17).Since there is the additional negative term p on the right-hand side of the equation on q, we have ∂p ∂σ = q ≤ 1 σ p ≤ 0. With this monotonicity it follows that the overdetermined problem at hand has a unique solution for a single value σ.More precisely, we choose a σ and solve equation (14) with the boundary conditions ( 15) and ( 16) to obtain p. Then we rewrite equation (12) as where we implicitly define the constant C γ .Plugging this expression for r in condition (11) we obtain In order to find (σ, p) such that also this condition is satisfied, we use that the right-hand side of ( 19) is increasing in σ and that p(0) is strictly decreasing in σ.
In this way we obtain a unique solution (σ, p).
Step 5.One can go further in the Hele-Shaw (incompressible) limit γ → ∞.In fact, equation ( 14) is then reduced to find a solution of which we have to match with the condition p σ (0) = − σ µ .Now we show that such a value of σ is unique.
The above semilinear elliptic problem is standard and there is a unique solution of (20).By the comparison principle (which holds true because G (•) < 0), we know that, for σ 1 < σ 2 , we have p σ1 < p σ2 .Additionally, as p σ < 0, we deduce that p σ < 0, and thus we conclude that the profile is decreasing.Furthermore, by integration From the monotonicity in σ of p σ and because G ≤ 0, we deduce that p σ (0) is an increasing function of σ, while σ → − σ µ is decreasing.This proves a unique possible match.
4. An extension with n ∞ = 0. We now provide a case where segregation does not hold but a sharp interface is still present.In particular, we consider the case where n does not vanish at infinity, that is, we impose the conditions with n ∞ > 0.
Proof.We use the same notations as in Section 3. In region I i.e. [0, ∞), we can rewrite equation (6) as Therefore, starting from a value n(0) = n 0 > n ∞ , n is decreasing to n ∞ with exponential decay.
Notice that we can write n(x) = N σ ν x with N the (parameter independent) solution of By an immediate monotonicity argument, we can fix uniquely n 0 through the mass condition In other words, there is a function N 0 such that As before, we complete in region II for x < 0. Since m(0 − ) = n 0 − n ∞ (no jump on p and on m + n), we can find again the jump condition on p (0) and thus achieve, see condition (13), Therefore, thanks to (22), we find Then, it remains to solve the equation for m, or equivalently p since n = n ∞ , that is (14).It comes with the overdetermined conditions p(−∞) = P M , p(0) = K γ N 0 M σ ν γ and (24).The monotonicity argument remains true and thus the solution is again unique and m is decreasing.
For µ > ν, we have The travelling wave solution of Theorem 4.1 is highly numerically unstable, and we have to impose a numerical constraint to capture it, that is, we solve instead of solving (1).The numerical solution obtained for µ < ν is shown in Fig. 6.If µ > ν, the numerical results summarised in Fig. 7 suggest that the travelling wave solution becomes unstable, and the sharp interface of Theorem 4.1 cannot propagate.

Figure 1 .
Figure 1.Numerical observations in the case µ < ν.Plots of the computed m (left panel) and n (right panel) at time t = 1 for ν = 2 and µ = 1.We observe the emergence of a spherical wave of dividing cells pushing the surrounding non-dividing cells (left panel), and an invasive front made of non-dividing cells that are induced to move by the expansion of dividing cells (right panel).

Figure 2 .
Figure 2. Numerical observations in the case µ > ν.Plots of the computed m (left panel) and n (right panel) at time t = 1 for ν = 1 and µ = 2.We observe the appearance of numerical instabilities which result in finger-like patterns of dividing cells (left panel) that protrude through and dislocate the surrounding nondividing cells (right panel).

3 . 9 )Figure 3 .
Figure 3.The profile of p for the travelling wave when n has a finite support that coincides with [0, r].

Figure 4 .
Figure 4. Travelling waves of Theorem 3.1 for µ < ν.Profiles of p (left panel), and m (right panel, red curve) and n (right panel, blue curve) for the travelling wave in the case where n has a compact support and µ < ν.The dashed line in the left panel highlights the value of P M , while the dashed line in the right panel highlights the value of P M /K γ 1/γ .

Figure 5 .
Figure 5. Transient regime of Theorem 3.1 for µ > ν.Profiles of p (left panel), and m (right panel, red curve) and n (right panel, blue curve) in the case where n has a compact support and µ > ν.The dashed line in the left panel highlights the value of P M , while the dashed line in the right panel highlights the value of P M /K γ 1/γ .This figure shows a transient regime after which n is left behind and m propagates alone (see also Supplementary Movie S1).

Figure 6 .
Figure 6.Travelling waves of Theorem 4.1 for µ < ν.Profiles of p (left panel), and m (right panel, red curve) and n (right panel, blue curve) for the travelling wave in the case where n does not vanish at infinity and µ < ν.The dashed line in the left panel highlights the value of P M , while the dashed line in the right panel highlights the value of P M /K γ 1/γ .

Figure 7 .
Figure 7. Transient regime of Theorem 4.1 for µ > ν.Profiles of p (left panel), and m (right panel, red curve) and n (right panel, blue curve) in the case where n does not vanish at infinity and µ > ν.The dashed line in the left panel highlights the value of P M , while the dashed line in the right panel highlights the value of P M /K γ 1/γ .