From individual-based mechanical models of multicellular systems to free-boundary problems

TL;DR

This paper develops a mechanical model for two interacting cell populations, deriving a free-boundary continuum model, proving its mathematical properties, and validating it through numerical simulations that align well with individual-based models.

Contribution

It introduces a novel link between discrete cell models and continuum free-boundary problems, including mathematical analysis and validation.

Findings

Continuum model accurately describes cell population dynamics.

Mathematical proof of existence and traveling-wave solutions.

Numerical results agree with individual-based simulations.

Abstract

In this paper we present an individual-based mechanical model that describes the dynamics of two contiguous cell populations with different proliferative and mechanical characteristics. An off-lattice modelling approach is considered whereby: (i) every cell is identified by the position of its centre; (ii) mechanical interactions between cells are described via generic nonlinear force laws; and (iii) cell proliferation is contact inhibited. We formally show that the continuum counterpart of this discrete model is given by a free-boundary problem for the cell densities. The results of the derivation demonstrate how the parameters of continuum mechanical models of multicellular systems can be related to biophysical cell properties. We prove an existence result for the free-boundary problem and construct travelling-wave solutions. Numerical simulations are performed in the case where the cellular interaction forces are described by the celebrated Johnson-Kendall-Roberts model of elastic contact, which has been previously used to model cell-cell interactions. The results obtained indicate excellent agreement between the simulation results for the individual-based model, the numerical solutions of the corresponding free-boundary problem and the travelling-wave analysis.

Figures

Equations

$$r_1(\tau )=s_0,\text{for all }\tau \ge 0.$$

$$(\Delta i{)}_i=\lceil \sum _{j=2}^{i}g(r_j-r_{j-1})\Delta \tau \rceil \text{ for }i=2,\dots ,m-1$$

$$r_i(\tau +\Delta \tau )=r_{i-(\Delta i{)}_i}(\tau )\text{ for }i=2,\dots ,m-1,$$

$$\frac{d{r}_i}{d\tau }$$

$$\frac{d{r}_1}{d\tau }$$

$$\frac{d{r}_m}{d\tau }$$

$$\frac{d{r}_n}{d\tau }$$

$$r_i(\tau )=r(\tau ,y_i)\text{with}{y}_i=i\delta ,$$

$$r_{i\pm 1}(\tau )=r(\tau ,y_{i\pm 1})=r(\tau ,y_i\pm \delta ),r_{i-(\Delta i{)}_i}(\tau )=r(\tau ,y_{i-(\Delta i{)}_i})=r(\tau ,y_i-(\Delta i{)}_i\delta ).$$

$$r(\tau ,y_1)=s_0,r(\tau ,y_m)=s_1(\tau ),r(\tau ,y_n)=s_2(\tau ),\text{for }\tau >0.$$

$$r(\tau +\Delta \tau ,y_i)=r(\tau ,y_i)+\frac{\partial r(\tau ,y_i)}{\partial \tau }\Delta \tau +o(\Delta \tau ),$$

$$r(\tau ,y_i-(\Delta i{)}_i\delta )=r(\tau ,y_i)-\frac{\partial r(\tau ,y_i)}{\partial y}(\Delta i{)}_i\delta +o(\delta ),$$

$$\lceil \sum _{j=2}^{i}g(r_j-r_{j-1})\Delta \tau \rceil \approx {\int }_{y_1}^{y_i}g\left(\frac{\partial r(\tau ,y^{\prime })}{\partial y^{\prime }}\delta \right)d{y}^{\prime }\frac{\Delta \tau }{\delta },$$

$$\frac{\partial r(\tau ,y_i)}{\partial \tau }\approx -\left({\int }_{y_1}^{y_i}g\left(\frac{\partial r(\tau ,y^{\prime })}{\partial y^{\prime }}\delta \right)d{y}^{\prime }\right)\frac{\partial r(\tau ,y_i)}{\partial y}\text{for }i=2,\dots ,m-1.$$

$$r(\tau ,y_i+\delta )=r(\tau ,y_i)+\frac{\partial r(\tau ,y_i)}{\partial y}\delta +\frac{1}{2}\frac{{\partial }^{2}r(\tau ,y_i)}{\partial y^2}{\delta }^2+o({\delta }^2),$$

$$r(\tau ,y_i-\delta )=r(\tau ,y_i)-\frac{\partial r(\tau ,y_i)}{\partial y}\delta +\frac{1}{2}\frac{{\partial }^{2}r(\tau ,y_i)}{\partial y^2}{\delta }^2+o({\delta }^2),$$

$$F_l(r_i-r_{i-1})-F_l(r_{i+1}-r_i)\approx -F_l^{\prime }\left(\frac{\partial r}{\partial y}\delta \right)\frac{{\partial }^{2}r}{\partial y^2}{\delta }^2,l=A,B.$$

$$\frac{\partial r}{\partial \tau }=-\frac{1}{{\eta }_A}{F}_A^{\prime }\left(\frac{\partial r}{\partial y}\delta \right)\frac{{\partial }^{2}r}{\partial y^2}{\delta }^2-\left({\int }_{y_1}^{y}g\left(\frac{\partial r}{\partial y^{\prime }}\delta \right)d{y}^{\prime }\right)\frac{\partial r}{\partial y}\text{ for }y\in (y_1,y_m),$$

$$\frac{\partial r}{\partial \tau }=-\frac{1}{{\eta }_B}{F}_B^{\prime }\left(\frac{\partial r}{\partial y}\delta \right)\frac{{\partial }^{2}r}{\partial y^2}{\delta }^2\text{ for }y\in (y_m,y_n).$$

$$\frac{\partial r}{\partial \tau }$$

$$\displaystyle\quad-\frac{1}{\eta_{B}}\left[F_{B}\Big{(}\frac{\partial r}{\partial y}\,\delta\Big{)}+\frac{1}{2}F^{\prime}_{B}\left(\frac{\partial r}{\partial y}\,\delta\right)\frac{\partial^{2}r}{\partial y^{2}}\,\delta^{2}\right]\;\;\quad\text{ at }\;y=y_{m},$$

$$\frac{\partial r}{\partial \tau }=\frac{1}{{\eta }_B}\left[F_B\left(\frac{\partial r}{\partial y}\delta \right)-\frac{1}{2}{F}_B^{\prime }\left(\frac{\partial r}{\partial y}\delta \right)\frac{{\partial }^{2}r}{\partial y^2}{\delta }^2\right]\text{ at }y=y_n.$$

$${\rho }_A(\tau ,y)={\left(\frac{\partial r}{\partial y}\delta \right)}^{-1}\text{mbox}fory\in [y_1,y_m],{\rho }_B(\tau ,y)={\left(\frac{\partial r}{\partial y}\delta \right)}^{-1}\text{mbox}fory\in [y_m,y_n].$$

$$\frac{j-i}{r_j(\tau )-r_i(\tau )}=\frac{j-i}{r(\tau ,y_j)-r(\tau ,y_i)}.$$

$${\rho }_l(\tau ,y_i)=\frac{1}{r_{i+1}(\tau )-r_i(\tau )}\approx \frac{1}{\delta }\frac{1}{\frac{\partial r(\tau ,y_i)}{\partial y}}.$$

$$\frac{\partial r}{\partial \tau }=-\frac{\partial r}{\partial y}\frac{\partial y}{\partial t}=-\frac{1}{\delta }{\frac{1}{\rho }}_l\frac{\partial y}{\partial t}.$$

$$\frac{\partial r}{\partial y}=\frac{1}{\delta }\frac{1}{{\rho }_l},\frac{{\partial }^{2}r}{\partial y^2}=\frac{1}{\delta }\frac{\partial }{\partial r}\left({\frac{1}{\rho }}_l\right)\frac{\partial r}{\partial y}=-\frac{1}{{\delta }^2}\frac{1}{{\rho }_l^3}\frac{\partial {\rho }_l}{\partial r},$$

$$\left({\int }_{y_1}^{y}g\left(\frac{\partial r}{\partial y^{\prime }}\delta \right)d{y}^{\prime }\right)\frac{\partial r}{\partial y}=\left({\int }_{s_0}^{r}g\left(1/{\rho }_A\right){\rho }_{A}d{r}^{\prime }\right)\frac{1}{{\rho }_A}$$

$$\frac{1}{\delta }\frac{1}{{\rho }_A}\frac{\partial y}{\partial t}=-\frac{F_A^{\prime }(1/{\rho }_A)}{{\eta }_A{\rho }_A^3}\frac{\partial {\rho }_A}{\partial r}+{\int }_{s_0}^{r}g\left(1/{\rho }_A\right){\rho }_{A}d{r}^{\prime }\frac{1}{{\rho }_A}$$

$$\frac{1}{\delta }\frac{1}{{\rho }_B}\frac{\partial y}{\partial t}=-\frac{F_B^{\prime }(1/{\rho }_B)}{{\eta }_B{\rho }_B^3}\frac{\partial {\rho }_B}{\partial r}+{\int }_{s_0}^{s_1}g\left(1/{\rho }_A\right){\rho }_{A}d{r}^{\prime }\frac{1}{{\rho }_B}.$$

$$\frac{\partial }{\partial r}\left(\frac{1}{\delta }\frac{\partial y}{\partial t}\right)=\frac{\partial }{\partial t}\left(\frac{1}{\delta }\frac{\partial y}{\partial r}\right)=\frac{\partial }{\partial t}{\left(\frac{\partial r}{\partial y}\delta \right)}^{-1}=\frac{\partial {\rho }_l}{\partial t},$$

$$\frac{\partial }{\partial r}{\int }_{s_0}^{r}g\left(1/{\rho }_A\right){\rho }_{A}d{r}^{\prime }=G\left({\rho }_A\right){\rho }_A\text{ and }\frac{\partial }{\partial r}{\int }_{s_0}^{s_1}g\left(1/{\rho }_A\right){\rho }_{A}d{r}^{\prime }=0,$$

Full text

1. Introduction

Continuum mechanical models of multicellular systems have been increasingly used to achieve a deeper understanding of the underpinnings of tissue development, wound healing and tumour growth [1, 2, 3, 6, 8, 10, 13, 15, 26, 32, 40, 44, 46, 52, 53]. These models are formulated in terms of nonlinear partial differential equations for cell densities (or cell volume fractions) and, as such, are amenable to both numerical and analytical approaches that enable insight to be gained into the biological system under study. From a mathematical perspective, over the past few years particular attention has been paid to the existence of travelling-wave solutions with composite shapes [5, 12, 31, 43, 48] and to the convergence to free-boundary problems in the asymptotic limit whereby cells are represented as an incompressible fluid [7, 30, 33, 41, 42].

Whilst continuum mechanical models of multicellular systems are usually defined on the basis of tissue-scale phenomenological considerations, off-lattice individual-based models enable representation of cell mechanics at the level of individual cells [18, 49]. However, as the numerical exploration of such individual-based models requires large computational times for biologically relevant cell numbers and the models are not analytically tractable, it is desirable to derive continuum models in an appropriate limit [4, 9, 17, 18, 25, 28, 29, 34, 35, 36, 37, 47]. Although mechanical interactions between interfacing cell populations with different characteristics arise in many biological contexts (e.g. tumour growth, development), relatively little prior work has explored the connection between off-lattice individual-based models and continuum models in such situations.

In this paper we propose an individual-based mechanical model for the dynamics of two contiguous cell populations with different proliferative and mechanical characteristics. In our model: (i) every cell is identified by the position of its centre; (ii) mechanical interactions between cells are described via generic nonlinear force laws; and (iii) cell proliferation is contact inhibited. Formally deriving a continuum counterpart of the discrete model, we obtain a free-boundary problem with nonstandard transmission conditions that governs the dynamics of the cell densities. Our derivation extends a previous method developed for the case of a single cell population [36, 37].

To prove an existence result for the free-boundary problem, a novel extension of methods previously developed for related free-boundary problems [21, 22, 23, 51] is required due to the specific structure of our boundary and transmission conditions. In particular, the jump in the density and in the flux across the moving interface between the two cell populations, along with the fact that there are no conditions of Dirichlet type, prevent us from using existing ideas which are partially based on continuity arguments [22, 23] and from applying the enthalpy method [16, 51].

Moreover, building on a recently presented method for a related system of nonlinear partial differential equations [12, 31], we also construct travelling-wave solutions for the free-boundary problem. In this respect, the novelty of our work lies in the fact that we consider a free-boundary problem subject to biomechanical transmission conditions which are different from those considered in [12, 31]. This requires a different approach when studying the properties of the solution at the interface between the two cell populations and introduces a significant difference in the qualitative properties of the travelling wave.

Numerical simulations are performed in the case where the cellular interaction forces are described by the celebrated Johnson-Kendall-Roberts (JKR) model of elastic contact [27], which has been shown to be experimentally accurate in some cases [14], and has been previously used to approximate mechanical interactions between cells [19, 20]. The results obtained support the findings of the travelling-wave analysis, and demonstrate excellent agreement between the individual-based model and the corresponding free-boundary problem.

The paper is organised as follows. In Section 2, we present our individual-based mechanical model and formally derive the corresponding free-boundary problem. In Section 3, we prove the existence result for the free-boundary problem. In Section 4, we develop the travelling-wave analysis. In Section 5, we compare simulation results for the individual-based model and numerical solutions of the free-boundary problem. Section 6 concludes the paper and provides a brief overview of possible research perspectives.

2. Formulation of the individual-based model and derivation of the corresponding free-boundary problem

2.1. Formulation of the individual-based model

We consider a one-dimensional multicellular system that consists of two populations of cells that are arranged along the real line $\mathbb{R}$ and characterised by different proliferative and mechanical properties. We label the two cell populations by the letters $A$ and $B$ and make the assumption that, during the considered time interval, the cells in population $A$ can proliferate, whereas the cells in population $B$ cannot. We denote the number of cells in population $B$ by $M>0$. Moreover, at time $\tau\geq 0$ we let the function $m(\tau)$ represent the number of cells in population $A$ and compute the total number of cells inside the system as $n(\tau)=m(\tau)+M$.

We adopt a discrete off-lattice modelling approach whereby every cell is identified by the position of its centre [49]. Building upon the ideas presented in [36, 37], we model the two cell populations as a chain of masses and springs with the masses corresponding to the cell centres, and assume the cell order to be fixed. We label each cell by an index $i=1,\ldots,n(\tau)$ and describe the position of the $i^{th}$ cell’s centre at time $\tau$ by means of the function $r_{i}(\tau)$. Without loss of generality, we let the cells of population $A$ be on the left of the cells of population $B$.

We assume that the centre of the first cell of population $A$ is pinned at a point $s_{0}\in\mathbb{R}$, i.e.

[TABLE]

We describe the effect of cell proliferation and mechanical interactions between cells on the dynamics of the multicellular system using the modelling strategies and the assumptions described hereafter.

Mathematical modelling of cell proliferation.

We assume that cell proliferation is contact dependent such that the proliferation rate $g$ of the $j^{th}$ cell in population $A$ depends on the position of neighbouring cells, i.e. $g\equiv g\left(r_{j}(\tau)-r_{j-1}(\tau)\right)$ with $j=2,\ldots,m(\tau)-1$. Hence in a time interval $\Delta\tau$ the proliferation of cells in population $A$ will result in a reindexing $\Delta i$ of the $i^{th}$ cell by an amount

[TABLE]

such that

[TABLE]

where $m\equiv m(\tau)$.

Mathematical modelling of mechanical interactions between cells.

We make the assumption that mechanical interactions between nearest neighbour cells depend on the distance between their centres. We denote the force exerted on the $i^{th}$ cell of population $l$ by its left and right neighbours by $F_{l}(r_{i}-r_{i-1})$ and $F_{l}(r_{i+1}-r_{i})$, respectively, and introduce the parameter $\eta_{l}>0$ to model the damping coefficient of cells in population $l$, where $l=A,B$. With this notation and neglecting cell-cell friction, the dynamics of the positions of the cell centres are described via the following system of differential equations

[TABLE]

where $m\equiv m(\tau)$ and $n\equiv n(\tau)$. We complete system (2.4) with the following differential equations

[TABLE]

2.2. Derivation of the corresponding free-boundary problem

In order to formally derive a continuum version of our individual-based mechanical model (2.4) and (2.1), considering the scenario where the number of cells in both populations is large, we introduce the continuous variable $y\in\mathbb{R}$ so that, for some $\delta>0$ sufficiently small,

[TABLE]

and

[TABLE]

Moreover, we use the notation

[TABLE]

We assume the function $r(\tau,y)$ to be continuously differentiable with respect to the variable $\tau$ and twice continuously differentiable with respect to the variable $y$. Under these assumptions, letting $\Delta\tau$ and $\delta$ be sufficiently small, and using the Taylor expansions

[TABLE]

along with the approximation

[TABLE]

from (2.3) we obtain

[TABLE]

Moreover, using the Taylor expansions

[TABLE]

and making the additional assumption that the functions $F_{A}$ and $F_{B}$ are twice continuously differentiable, we approximate the force terms in (2.4) for $i=2,\ldots,n-1$ as

[TABLE]

Using the approximations (2.7) and (2.8), and combining proliferation (2.3) and mechanical interaction (2.4) processes, we obtain

[TABLE]

and

[TABLE]

Similarly, mechanical interaction processes (2.1) for $i=1$, $i=m$ and $i=n$ yield, respectively,

[TABLE]

and

[TABLE]

Based on the ideas presented in [36, 37] and according to the considerations given in Remark 2.1, we define the cell number densities of populations $A$ and $B$ as

[TABLE]

Remark 2.1**.**

The definitions of the cell densities given by (2.13) are based on the observation that, at any time $\tau$, the quotient of the number of cells in a generic interval $[r_{i},r_{j}]$, with $j>i$, and the length of the interval is

[TABLE]

From the above relation, choosing $j=i+1$ and using the fact that $\delta$ is small, we obtain the following approximate expression for the cell density

[TABLE]

The change of coordinates $(\tau,y)\mapsto(t,r)$, with $t=\tau$ and $r=r(\tau,y)$ yields [36]

[TABLE]

Substituting this relation along with the expressions

[TABLE]

into equations (2.9) and (2.10) yields

[TABLE]

and

[TABLE]

Multiplying equations (2.15) and (2.16) by $\rho_{A}$ and $\rho_{B}$, respectively, differentiating with respect to $r$, using the fact that

[TABLE]

[TABLE]

with the growth rate $G$ of population $A$ defined as

[TABLE]

and renaming $r$ to $x$, we obtain the following equations for the cell densities $\rho_{A}(t,x)$ and $\rho_{B}(t,x)$

[TABLE]

where

[TABLE]

Similarly, the evolution equations for the positions of the free boundaries $s_{1}(t)$ and $s_{2}(t)$ are obtained from equations (2.11) and (2.12), respectively, yielding

[TABLE]

In order to obtain the boundary conditions (i.e. the conditions at $s_{0}$ and $s_{2}(t)$) and the transmission conditions (i.e. the conditions at $s_{1}(t)$), that are needed to complete the problem, we consider the mass balance equations

[TABLE]

Using the fact that $\displaystyle{\frac{ds_{0}}{dt}=0}$, together with equations (2.18) and (2.19), yields

[TABLE]

and

[TABLE]

Hence

[TABLE]

and

[TABLE]

The above equations along with equations (2.21) give

[TABLE]

We complement the free-boundary problem (2.18), (2.19) and (2.22) with the following initial conditions for the moving boundaries $s_{1}(t)$ and $s_{2}(t)$, and the cell densities $\rho_{A}(t,x)$ and $\rho_{B}(t,x)$:

[TABLE]

where $s_{0}<s_{1}^{\ast}<s_{2}^{\ast}$. Letting $\rho^{{\rm eq}}_{l}>0$ denote the equilibrium cell density (i.e. the density below which intercellular forces are zero) and $\rho^{\rm M}>\rho^{{\rm eq}}_{A}$ a critical cell density above which cells stop dividing due to contact inhibition, throughout the rest of the paper we will make the following assumptions:

[TABLE]

[TABLE]

where $l=A,B$, and

[TABLE]

Assumptions (2.25) and (2.26), together with notations (2.17) and (2.20), imply that the nonlinear diffusion coefficient $D_{l}(\rho_{l})$ and the growth rate $G(\rho_{A})$ are such that

[TABLE]

and

[TABLE]

3. An existence result for the free-boundary problem

Due to the specific structure of our boundary and transmission conditions, the existing well-posedness results for one-dimensional free-boundary problems, such as those presented in [21, 22, 23, 51], are not directly applicable to our problem. Therefore, in this section we prove an existence result for the free-boundary problem (2.18), (2.19), (2.22) and (2.23).

Assumption 3.1**.**

We make the following assumptions on the force terms $F_{A}$ and $F_{B}$, the diffusion coefficients $D_{A}$ and $D_{B}$, the growth rate $G$, and the initial conditions $\rho_{A}^{0}$ and $\rho_{B}^{0}$.

• (i)

The force terms $F_{l}\in H^{1}(0,\infty)\cap C^{3}(\rho_{l}^{\rm eq},\infty)$, with $l=A,B$, and satisfy (2.25).

• (ii)

The diffusion coefficients $D_{l}\in C^{2}(\rho_{l}^{\rm eq},\infty)\cap L^{\infty}(0,R)$, for any $R>0$, with $l=A,B$, are defined by (2.20), satisfy assumptions (2.27), and $D_{l}(\xi)\geq d_{l}>0$ for $\xi>\rho_{l}^{\rm eq}$.

• (iii)

The growth rate $G\in C^{3}(\mathbb{R})$ and satisfies assumptions (2.28).

• (iv)

The initial conditions $\rho_{A}^{0},\rho_{B}^{0}\in C^{3}_{0}(\mathbb{R})$ and satisfy assumptions (2.24).

Throughout this section we use the notation

[TABLE]

for $t\in[0,T]$, with $T>0$, and consider solutions in the sense specified by the following definition.

Definition 3.2**.**

A solution of the free-boundary problem (2.18), (2.19), (2.22), (2.23) is given by functions $s_{1},s_{2}\in W^{1,3}(0,T)$ and $\rho_{l}\in H^{1}(0,T;L^{2}(\Omega_{l}(t)))\cap L^{2}(0,T;H^{2}(\Omega_{l}(t)))$, with $\rho_{l}\geq\rho^{{\rm eq}}_{l}$ and $\rho_{l}\in L^{\infty}(0,T;L^{\infty}(\Omega_{l}(t)))$ for $l=A,B$, that satisfy equations (2.18) and (2.19), the following boundary and transmission conditions

[TABLE]

the equations for the free boundaries

[TABLE]

and the initial conditions (2.23).

Theorem 3.3**.**

Under Assumptions 3.1 there exists a solution of the free-boundary problem (2.18), (2.19), (2.22) and (2.23).

Proof.

In order to prove the existence of a solution of the free-boundary problem (2.18), (2.19), (2.22) we consider iterations over successive time intervals and use a fixed point argument. In particular, we first show the existence of a solution on a time interval $[0,T_{1}]$ such that

[TABLE]

Subsequently, the boundedness of $s_{1}^{\prime}$ and $s_{2}^{\prime}$, shown at the end of the proof, will allow iteration over successive time intervals in order to obtain an existence result for $t\in(0,T]$.

We begin by making the change of variables

[TABLE]

with $\zeta,\xi\in C^{2}_{0}(\mathbb{R})$ such that $\zeta(y)=1$ for $|y-s_{1}^{\ast}|<\alpha$ and $\zeta(y)=0$ for $|y-s_{1}^{\ast}|>2\alpha$, while $\xi(y)=1$ for $|y-s_{2}^{\ast}|<\alpha$ and $\xi(y)=0$ for $|y-s_{2}^{\ast}|>2\alpha$, where $\alpha=\min\{(s_{1}^{\ast}-s_{0})/4,(s_{2}^{\ast}-s_{1}^{\ast})/4\}$. The change of variables (3.4) transforms the time-dependent domains $\Omega_{A}(t)=(s_{0},s_{1}(t))$ and $\Omega_{B}(t)=(s_{1}(t),s_{2}(t))$ into the fixed intervals $\Omega_{A}^{\ast}=(s_{0},s_{1}^{\ast})$ and $\Omega_{B}^{\ast}=(s_{1}^{\ast},s_{2}^{\ast})$, respectively. A similar change of variables was considered in [21]. Notice that, for $s_{1}(t)$ and $s_{2}(t)$ satisfying conditions (3.3), such a change of variables defines a diffeomorphism from $[0,+\infty)$ into $[0,+\infty)$. Hence we obtain

[TABLE]

where $w_{A}$ and $w_{B}$ satisfy the reaction-diffusion-convection equations

[TABLE]

complemented with the nonlinear transmission and boundary conditions

[TABLE]

and the equations for the velocities $s^{\prime}_{1}$ and $s^{\prime}_{2}$

[TABLE]

Notice that for ease of notation we denote

[TABLE]

$G_{A}(w_{A})=G(w_{A})w_{A}$, and the functions $D_{l}(w_{l})$ are defined in terms of $F_{l}(w_{l})$ according to (2.20).

In equations (3.5), since $\xi(y)=0$ for $y<s_{2}^{\ast}-2\alpha$ and $s_{1}(t)<s_{2}^{\ast}-2\alpha$ for $t\in[0,T_{1}]$, we have

[TABLE]

The assumptions on $F_{l}$ and $D_{l}$, for $l=A,B$, ensure that

[TABLE]

Notice that, without loss of generality, we can focus on the case where $\rho_{l}^{0}>\rho_{l}^{\rm eq}$ for $l=A,B$. In fact, if $\rho_{l}^{0}=\rho_{l}^{\rm eq}$ the growth term in the equation for $w_{A}$ would result into $w_{A}(t,y)>\rho_{A}^{\rm eq}$ and $F_{A}(w_{A})>0$, thus ensuring that $w_{B}(t,y)>\rho_{B}^{\rm eq}$ due to the transmission conditions at $s_{1}^{\ast}$ and the convection term in the equation for $w_{B}$. In the case where $\rho_{l}^{0}>\rho_{l}^{\rm eq}$ for $l=A,B$, using the maximum principle and relations (3.8) we obtain that $w_{l}(t,x)>\rho_{l}^{\rm eq}$ for $(t,x)\in(0,T)\times\Omega_{l}^{\ast}$. Therefore, we conclude that system (3.6), (3.7), or equivalently system (2.18), (2.19), is nondegenerate. Notice also that assuming $s_{j}(t)\in H^{2}(0,T)$ with $s_{j}^{\prime}(t)\geq 0$ for $t\in[0,T]$ and $j=1,2$, and considering $F_{l}(w_{l})$, $\partial_{y}^{2}\partial_{t}F_{l}(w_{l})$ and $\partial_{t}^{2}w_{l}$ as test functions in equations (3.6) and (3.7), one can prove that $w_{l}$ is continuous in $\overline{\Omega}_{l,T}^{\ast}$, while $\partial_{t}w_{l}$ and $\partial_{y}^{2}F_{l}(w_{l})$ are continuous in $(0,T)\times\Omega_{l}^{\ast}$ for $l=A,B$, which is the regularity required to apply the maximum principle. A similar approach was previously used in the analysis of the porous medium equation [50].

The assumptions on $F_{B}$ imply that

[TABLE]

and, applying the maximum principle to the equation for $w_{B}$, we find that $w_{B}$ has a minimum at $s_{2}^{\ast}$ and a maximum at $s_{1}^{\ast}$. Hence, $\partial_{y}w_{B}(t,y)\leq 0$ at $y=s_{1}^{\ast}$ for $t>0$ and, therefore, $\partial_{y}F_{A}(w_{A})\leq 0$ at $y=s_{1}^{\ast}$. Applying the comparison principle and using the fact that $F^{\prime}_{A}(w_{A})=0$ for $0<w_{A}(t,y)\leq\rho_{A}^{\rm eq}$, along with the assumptions on $G$ and on the initial conditions, we obtain

[TABLE]

where $\rho^{\rm M}_{A}=\max\{\rho^{\rm M},\max\limits_{x\in[s_{0},s_{1}^{\ast}]}\rho_{A}^{0}(x)\}$. Moreover, applying the maximum principle to the equation for $w_{B}$ and using the assumptions on $F_{B}$ and on the initial data, along with the boundedness of $w_{A}$ and the transmission conditions at $y=s_{1}^{\ast}$, yield

[TABLE]

where $\rho_{B}^{\rm M}=\max\{F_{B}^{-1}(F_{A}(\rho^{\rm M}_{A})),\max\limits_{x\in[s_{1}^{\ast},s_{2}^{\ast}]}\rho_{B}^{0}(x)\}$. Using these results along with the change of variables given by equation (3.4) we conclude that

[TABLE]

If $w_{B}$ is nonconstant in $(s_{1}^{\ast},s_{2}^{\ast})$ and $w_{B}(t,s_{j}^{\ast})\neq\rho_{B}^{\rm eq}$ for $j=1,2$ and $t\geq 0$, the maximum principle yields $\partial_{y}w_{B}(t,s_{2}^{\ast})<0$ and $\partial_{y}w_{B}(t,s_{1}^{\ast})<0$. This along with the assumptions on $F_{B}$ ensures the monotonicity of the free boundaries $\{x=s_{1}(t)\}$ and $\{x=s_{2}(t)\}$, i.e.

[TABLE]

To prove the existence of a solution of problem (3.5)-(3.7) we use a fixed point argument. Let

[TABLE]

which are both well-defined quantities due to the assumptions on $F_{l}$ and $\rho_{l}^{0}$, for $l=A,B$. Moreover, consider

[TABLE]

for $l=A,B$, some constant $\mu>0$, and $T_{1}>0$. Notice that for $(s_{1},s_{2})\in\mathcal{W}_{s}$ we have

[TABLE]

Therefore, choosing

[TABLE]

we find that $s_{j}$ satisfies the conditions (3.3), for $j=1,2$, and the change of coordinates (3.4) is well-defined for all $(s_{1},s_{2})\in\mathcal{W}_{s}$.

For some given $(\tilde{s}_{1},\tilde{s}_{2})\in\mathcal{W}_{s}$ and $\widetilde{w}_{l}\in\mathcal{W}_{l}$, with $l=A,B$, we first consider the problem given by the following equations for $w_{A}$ and $w_{B}$

[TABLE]

For $(\tilde{s}_{1},\tilde{s}_{2})\in\mathcal{W}_{s}$ and $\widetilde{w}_{l}\in L^{2}(0,T_{1};H^{1}(\Omega_{l}^{\ast}))$, with $l=A,B$, applying the Rothe-Galerkin method and using the a priori estimates obtained by considering $F_{l}(w_{l})/{\widetilde{w}_{l}}$ as a test function in the equations for $w_{l}$, we obtain the existence of a weak solution $F_{l}(w_{l})\in L^{2}(0,T_{1};H^{1}(\Omega_{l}^{\ast}))$, with $\partial_{t}w_{l}\in L^{2}(0,T_{1};H^{-1}(\Omega_{l}^{\ast}))$, of problem (3.12). Notice that for $\tilde{s}_{1},\tilde{s}_{2}\in H^{2}(0,T_{1})$, in the same way as below, we can show that the solutions of (3.12) satisfy the regularity properties required by the maximum principle, and obtain that the solutions are bounded and satisfy (3.9) and (3.10), and the equations in (3.12) are nondegenerate.

To derive a priori estimates for $\partial_{t}w_{A}$ and $\partial_{t}w_{B}$, we consider $\phi=\partial_{t}F_{A}(w_{A})/\widetilde{w}_{A}$ and $\psi=\partial_{t}F_{B}(w_{B})/\widetilde{w}_{B}$ as test functions for the equations in problem (3.12). In this way, we obtain

[TABLE]

for $\tau\in(0,T_{1}]$, where $J_{A}(\tilde{s})=J_{A}(\tilde{s}_{1})$, $J_{B}(\tilde{s})=J_{B}(\tilde{s}_{1},\tilde{s}_{2})$, $R_{A}(\tilde{s})=R_{A}(\tilde{s}_{1})$, $R_{B}(\tilde{s})=R_{B}(\tilde{s}_{1},\tilde{s}_{2})$, $Q_{A}(\tilde{s}^{\prime})=Q_{A}(\tilde{s}_{1}^{\prime})$, and $Q_{B}(\tilde{s}^{\prime})=Q_{B}(\tilde{s}_{1}^{\prime},\tilde{s}_{2}^{\prime})$.

The transmission conditions in problem (3.12) ensure that the integral at $y=s_{1}^{\ast}$ is equal to zero, while for the integral at $y=s_{2}^{\ast}$ we have

[TABLE]

From the equation for $w_{A}$ in problem (3.12) we obtain

[TABLE]

Using the definition of $Q_{A}$ and Hölder inequality, the third term on the right-hand side is estimated as

[TABLE]

for any fixed $\delta>0$. The assumptions on $G_{A}$ and $F_{A}$, the boundedness of $J_{A}(\tilde{s}_{1})$, $R_{A}$, and $\partial_{y}R_{A}$, and the fact that $R_{A}^{2}(\tilde{s}_{1})\geq 4/9$ and $F_{A}^{\prime}(w_{A})=D_{A}(w_{A})\geq d_{A}>0$ for $w_{A}>\rho_{A}^{\rm eq}$, imply

[TABLE]

for $\tau\in(0,T_{1}]$. Notice that for $s_{1}^{\prime}\in L^{\infty}(0,T_{1})$ we would have the $L^{2}$-norm of $\partial_{y}F_{A}(w_{A})$ on the right-hand side of the last inequality. A similar inequality for $\|\partial_{y}^{2}F_{B}(w_{B})\|^{2}_{L^{2}(\Omega_{B,\tau}^{\ast})}$ follows from the equation for $w_{B}$ in problem (3.12). The Gagliardo-Nirenberg inequality gives

[TABLE]

and, using the fact that $w_{l}$ is uniformly bounded and choosing $\delta>0$ sufficiently small, we obtain

[TABLE]

for $\tau\in(0,T_{1}]$, where $|\tilde{s}^{\prime}|=|\tilde{s}_{1}^{\prime}|$ if $l=A$ and $|\tilde{s}^{\prime}|=|\tilde{s}_{1}^{\prime}|+|\tilde{s}_{2}^{\prime}|$ if $l=B$. In a similar way, we also obtain the following pointwise in the time variable estimate

[TABLE]

for a.e. $t\in[0,T_{1}]$. Additionally, using the Gagliardo-Nirenberg inequality we have

[TABLE]

We shall estimate each term in equation (3.13) separately. Using estimates (3.14) and (3.15) yields

[TABLE]

Notice that assuming the boundedness of $\tilde{s}^{\prime}_{j}$, with $j=1,2$, we would have the $L^{2}(0,T;L^{2}(\Omega_{l}^{\ast}))$-norm instead of the $L^{6}(0,T;L^{2}(\Omega_{l}^{\ast}))$-norm of $\partial_{y}F_{l}(w_{l})$ in the last inequality. Using the results in (3.15) and (3.17), along with the boundedness of $w_{l}$ and $\widetilde{w}_{l}$, we estimate the next term in (3.13) as

[TABLE]

For the fourth integral in (3.13) we have

[TABLE]

for $\tau\in(0,T_{1}]$ and any fixed $\delta>0$. Here we used the following estimate

[TABLE]

along with estimate (3.15). Using (3.17) and the boundedness of $w_{l}$ we also obtain

[TABLE]

The boundedness of $\widetilde{w}_{A}$, along with the assumptions on $G_{A}$, implies

[TABLE]

for $\tau\in(0,T_{1}]$ and any fixed $\delta>0$.

Thus for $\partial_{t}\widetilde{w}_{l}\in L^{2}((0,T_{1})\times\Omega_{l}^{\ast})$ and $\partial_{y}\widetilde{w}_{l}\in L^{6}(0,T_{1};L^{4}(\Omega_{l}^{\ast}))$, combining the estimates from above, choosing $\delta>0$ sufficiently small, and applying the Gronwall inequality yields

[TABLE]

In a similar way we also obtain

[TABLE]

Thus using (3) and (3), along with (3.15) and (3.17), and considering $T_{1}$ sufficiently small, we find that the map $K:\mathcal{W}_{A}\times\mathcal{W}_{B}\to\mathcal{W}_{A}\times\mathcal{W}_{B}$, where $(w_{A},w_{B})=K(\widetilde{w}_{A},\widetilde{w}_{B})$ is defined as a solution of problem (3.12) for a given $(\tilde{s}_{1},\tilde{s}_{2})\in\mathcal{W}_{s}$, is continuous.

Considering $w_{l}$ in equation (3.13) instead of $\widetilde{w}_{l}$ and using the boundedness of $w_{l}$ yield

[TABLE]

and

[TABLE]

for $\tau\in(0,T_{1}]$. Choosing $\delta>0$ sufficiently small, applying the Gronwall inequality, and considering $T_{1}$ such that

[TABLE]

we obtain the following estimates for $w_{A}$ and $w_{B}$

[TABLE]

The estimate for $\|\partial_{t}\partial_{y}w_{l}\|_{L^{2}(0,T_{1};H^{-1}(\Omega_{l}^{\ast}))}$ in terms of $\|\partial_{y}F_{l}(\widetilde{w}_{l})\|_{L^{4}(\Omega_{l,T_{1}}^{\ast})}$ and $\|\partial_{t}\widetilde{w}_{l}\|_{L^{2}(\Omega_{l,T_{1}}^{\ast})}$ follows directly from differentiating the equation for $w_{l}$ with respect to $y$ and using the boundedness of $\|\partial_{y}^{2}w_{l}\|_{L^{2}(\Omega_{l,T_{1}}^{\ast})}$ and $\|\partial_{y}F_{l}(w_{l})\|_{L^{\infty}(0,T_{1};L^{2}(\Omega_{l}^{\ast}))}$, which is ensured by (3) and (3.15). Using (3.20) and (3.15) and differentiating the equation for $w_{l}$ in (3.12) with respect to $y$, while considering $w_{l}$ instead of $\widetilde{w}_{l}$, gives

[TABLE]

Thus for a sufficiently small $T_{1}$, or small initial data, $(w_{A},w_{B})=K(w_{A},w_{B})$ is uniformly bounded in $\mathcal{W}_{A}\times\mathcal{W}_{B}$ and $\partial_{y}w_{l}$ in $\mathcal{V}_{l}$, for $l=A,B$, where

[TABLE]

The Aubin-Lions lemma, along with the Sobolev embedding theorem, ensures that $\mathcal{V}_{l}$ is a compact subset of $L^{2}(\Omega_{l,T_{1}}^{\ast})$ and of $L^{2}(0,T_{1};C(\overline{\Omega}_{l}^{\ast}))$, for $l=A,B$. Using inequality (3.17) we also obtain that the embedding $\mathcal{V}_{l}\subset L^{6}(0,T_{1};L^{4}(\Omega_{l}^{\ast}))$ is compact. Thus applying the Schauder fixed point theorem, see e.g. [45], gives that for a given pair $(s_{1},s_{2})\in\mathcal{W}_{s}$ there exists a solution of problem (3.5) and (3.6) for $t\in(0,T_{1}]$, with an appropriate choice of $T_{1}>0$.

To complete the proof we shall show that $\mathcal{M}:L^{3}(0,T_{1})^{2}\to L^{3}(0,T_{1})^{2}$ given by

[TABLE]

where $s_{j}(t)=s_{j}^{\ast}+\int_{0}^{t}r_{j}(\tau)d\tau$ for $j=1,2$, maps $\mathcal{W}^{\prime}_{s}=\{(r_{1},r_{2})\in L^{3}(0,T_{1})^{2}:\|r_{j}-s_{j}^{\ast,1}\|_{L^{3}(0,T_{1})}\leq 1\}$ into itself and is precompact. Considering $(r_{1},r_{2})\in\mathcal{W}^{\prime}_{s}$ we have

[TABLE]

for an appropriate choice of $\delta>0$ and $T_{1}>0$. To show that $\mathcal{M}:L^{3}(0,T_{1})^{2}\to L^{3}(0,T_{1})^{2}$ is precompact we consider two sequences $\{r_{1}^{n}\}$ and $\{r_{2}^{n}\}$ bounded in $L^{3}(0,T_{1})$ and obtain

[TABLE]

for $l=A,B$, with a constant $C$ independent of $n$. Using the fact that the embedding $H^{1}(\Omega_{l}^{\ast})\subset C(\overline{\Omega}_{l}^{\ast})$ is compact and applying the Aubin-Lions lemma we obtain the strong convergence $w_{l}^{n}\to w_{l}$ in $L^{p}(0,T_{1};C(\overline{\Omega}_{l}^{\ast}))$, for any $1<p<\infty$, and $\partial_{y}F_{l}(w_{l}^{n})\to\partial_{y}F_{l}(w_{l})$ in $L^{2}(0,T_{1};C(\overline{\Omega}_{l}^{\ast}))$ as $n\to\infty$. This combined with the estimate (3.22) ensures that

[TABLE]

as $n\to\infty$, where we used the fact that

[TABLE]

The convergence in (3) implies the strong convergence $\mathcal{M}(r^{n}_{1},r^{n}_{2})\to\mathcal{M}(r_{1},r_{2})$ in $L^{3}(0,T_{1})^{2}$ as $n\to\infty$. Hence, we have proved the existence of a solution of problem (3.5)-(3.7) in $(0,\hat{T})\times\Omega_{l}^{\ast}$ with $\hat{T}=\min\{T_{1},T_{2}\},$ where

[TABLE]

Now we show that $s^{\prime}_{1}(t)$ and $s^{\prime}_{2}(t)$ are uniformly bounded, which will allow us to iterate over successive time intervals and obtain that $\hat{T}\leq T_{2}$. The uniform boundedness of $\rho_{B}$, the assumptions on $F_{B}$, and equations (3.2) ensure that $s^{\prime}_{2}(t)$ is uniformly bounded. To show the boundedness of $s_{1}^{\prime}(t)$, we consider the original problem (2.18), (2.19), (2.22), and (2.23) and apply the comparison principle to the following problem for $v=\eta_{A}^{-1}F_{A}(1/\rho_{A})$:

[TABLE]

where $\tilde{G}_{A}(v)=G_{A}([F_{A}^{-1}(\eta_{A}v)]^{-1})$ and $\tilde{D}_{A}(v)=D_{A}([F_{A}^{-1}(\eta_{A}v)]^{-1})$.

Consider the interval $(s_{1}(t)-\delta,s_{1}(t))$, with $t>0,$ and the function

[TABLE]

for some $\delta>0$. A similar idea was used in [21]. Since $F_{A}$ and $F_{B}$ are monotonically decreasing functions, we obtain

[TABLE]

For the derivatives of $\omega$, since $s^{\prime}_{1}(t)\geq 0$, we have

[TABLE]

Using the assumptions on $G_{A}$, for $\delta>0$ sufficiently small we obtain

[TABLE]

Since $F_{A}$ is continuous and $\eta^{-1}_{A}F_{A}(1/\rho_{A}^{0})=\eta^{-1}_{B}F_{B}(1/\rho_{B}^{0})$ at $x=s_{1}^{\ast}$, there exists a sufficiently small $\delta>0$ such that

[TABLE]

Then applying the comparison principle for parabolic equations gives

[TABLE]

Hence we have

[TABLE]

and for some sufficiently small fixed $\delta>0$

[TABLE]

Therefore, provided that $\partial_{x}F_{l}(1/\rho_{l})$ is uniformly bounded in $L^{\infty}(0,T;L^{2}(\Omega_{l}(t)))$, for $l=A,B$, the uniform boundedness of $s_{1}^{\prime}$ and $s_{2}^{\prime}$ allows us to iterate over successive time intervals and prove the existence of a global solution of the free-boundary problem (2.18), (2.19), (2.22) and (2.23).

Thus, as next we prove the uniform boundedness of $\|\partial_{x}F_{l}(1/\rho_{l})\|_{L^{\infty}(0,T;L^{2}(\Omega_{l}(t)))}$, or equivalently of $\|\partial_{y}F_{l}(w_{l})\|_{L^{\infty}(0,T;L^{2}(\Omega_{l}^{\ast}))}$, for $l=A,B$. First we show higher regularity of the solutions of problem (3.5)-(3.7) by differentiating the equations in (3.5) with respect to the time variable and considering $\phi=\partial_{t}^{2}F_{A}(w_{A})/w_{A}$ and $\psi=\partial_{t}^{2}F_{B}(w_{B})/w_{B}$ as test functions, respectively,

[TABLE]

The second term in the equation above can be estimated as

[TABLE]

for $\tau\in(0,T]$. For the third and fourth terms we have

[TABLE]

Moreover,

[TABLE]

We estimate the next terms as

[TABLE]

The reaction term is estimated by

[TABLE]

For the non-zero contributions from the boundary terms we have

[TABLE]

for $l=A,B$, where

[TABLE]

and

[TABLE]

Hence, applying the Gronwall inequality and using the estimates (3.20) and (3.21) yields

[TABLE]

with constants $C_{j}$, for $j=1,2,3$, depending on $\|s_{1}^{\prime}\|_{L^{\infty}(0,\hat{T})}$ and $\|s_{2}^{\prime}\|_{L^{\infty}(0,\hat{T})}$.

Using the Gagliardo-Nirenberg inequality and boundedness of $\partial_{y}F_{A}(w_{A}(t,s_{1}^{\ast}))$ we obtain

[TABLE]

Differentiating (3.5) with respect to the time variable and using definition of $J_{A}$, $R_{A}$, and $Q_{A}$ yields

[TABLE]

The last term is estimated by

[TABLE]

Similar estimates, using the boundedness of $\partial_{y}F_{B}(w_{B}(t,s_{2}^{\ast}))$, hold for $\partial_{t}\partial_{y}^{2}F_{B}(w_{B})$ and $s^{\prime\prime}_{2}$. Thus choosing $\delta>0$ sufficiently small and using the boundedness of $s^{\prime}_{1}$ and $s^{\prime}_{2}$ we obtain

[TABLE]

and, considering the same calculations as in the derivation of (3.25), conclude

[TABLE]

where the constant $C$ depends on $\|\partial_{y}w_{l}\|_{L^{\infty}(0,\hat{T};L^{2}(\Omega^{\ast}_{l}))}$, $\|\partial_{t}w_{l}\|_{L^{2}(\Omega^{\ast}_{l,\hat{T}})}$ for $l=A,B$, and $\|s_{j}^{\prime}\|_{L^{\infty}(0,\hat{T})}$ for $j=1,2$.

Hence estimates (3.21) and (3.26) imply that $\partial_{y}F_{l}(w_{l})\in C(\overline{\Omega}_{l,\hat{T}}^{\ast})$. Since $\partial_{y}F_{A}(w_{A}(t,s_{1}^{\ast}))$ and $\partial_{y}F_{B}(w_{B}(t,s_{j}^{\ast}))$, for $j=1,2$, are bounded, then $\partial_{y}F_{A}(w_{A})$ is bounded in $[s_{1}^{\ast}-\delta,s_{1}^{\ast}]$ and $\partial_{y}F_{B}(w_{B})$ is bounded in $[s_{1}^{\ast},s_{1}^{\ast}+\delta]$ and $[s_{2}^{\ast}-\delta,s_{2}^{\ast}]$, for $t\in[0,\hat{T}]$ and for a sufficiently small $\delta>0$. Taking $F_{l}(w_{l})$ as test function in (3.5) and using boundedness of $\partial_{y}F_{A}(w_{A}(t,s_{1}^{\ast}))$ and $\partial_{y}F_{B}(w_{B}(t,s_{j}^{\ast}))$, for $j=1,2$, yield

[TABLE]

Considering a cut-off function $\zeta_{l}\in C^{2}(\overline{\Omega}_{l}^{\ast})$, such that $\zeta_{A}(y)=1$ for $y\in[s_{0},s_{1}^{\ast}-\delta/2]$ and $\zeta_{A}(y)=0$ for $y\in[s_{1}^{\ast}-\delta/4,s_{1}^{\ast}]$ and $\zeta_{B}(y)=1$ for $y\in[s_{1}^{\ast}+\delta/2,s_{2}^{\ast}-\delta/2]$ and $\zeta_{B}(y)=0$ for $y\in[s_{1}^{\ast},s_{1}^{\ast}+\delta/4]$ and $y\in[s_{2}^{\ast}-\delta/4,s_{2}^{\ast}]$ and taking $\partial_{t}F_{l}(w_{l})\zeta_{l}^{2}$ as test function in (3.5) gives

[TABLE]

Therefore, $\|\partial_{y}F_{A}(w_{A})\|_{L^{\infty}(0,\hat{T};L^{2}(\Omega^{\ast}_{A}))}$ and $\|\partial_{y}F_{B}(w_{B})\|_{L^{\infty}(0,\hat{T};L^{2}(\Omega^{\ast}_{B}))}$ are uniformly bounded and $\hat{T}$ can be chosen independently of the initial data. Then iterating over $\hat{T}$, we can conclude that there exists a global solution of problem (2.18), (2.19), (2.22), (2.23).

∎

4. Travelling-wave solutions of the free-boundary problem

In this section, we carry out a travelling-wave analysis for the free-boundary problem (2.18), (2.19) and (2.22).

We begin by noting that under assumptions (2.25), (2.27) and (2.28), the free-boundary problem (2.18), (2.19) and (2.22) can be written in terms of the cell pressures $P_{A}$ and $P_{B}$ defined according to the following barotropic relation

[TABLE]

Under assumptions (2.27), the barotropic relation (4.1) is such that

[TABLE]

The monotonicity conditions (4.2) allow one to write both the force terms $F_{l}(1/\rho_{l})$ and the growth rate $G(\rho_{A})$ as functions of the cell pressure $\tilde{F}_{A}(P_{A})$, $\tilde{F}_{B}(P_{B})$ and $\tilde{G}(P_{A})$. Moreover, with the notation

[TABLE]

the monotonicity conditions (4.2) make it possible to reformulate assumptions (2.25) and (2.28) on the functions $F_{l}$ and $G$ as

[TABLE]

and

[TABLE]

Hence, under assumptions (2.25), (2.27) and (2.28), using the barotropic relation (4.1), we can rewrite the free-boundary problem (2.18), (2.19) and (2.22) in the following alternative form:

[TABLE]

Having rewritten the problem in this form allows us to construct travelling-wave solutions using an approach that builds on the method of proof recently presented in [12, 31]. For the sake of brevity, in this section we drop the tildes from all the quantities in problem (4.6).

We construct travelling-wave solutions of the free-boundary problem (4.6) such that both the position of the inner free boundary, $s_{1}(t)$, and the position of the outer free boundary, $s_{2}(t)$, move at a given constant speed $c>0$. Without loss of generality, we let $s_{0}$ go to $-\infty$ and consider the case where

[TABLE]

for some $\ell>0$, so that

[TABLE]

We make the following travelling-wave ansatz for the cell densities $\rho_{A}$ and $\rho_{B}$

[TABLE]

which are related to the cell pressures $P_{A}(z)$ and $P_{B}(z)$ through the barotropic relation (4.1). In this framework, substituting the travelling-wave ansatz (4.8) into problem (4.6) we find

[TABLE]

where $\rho^{\prime}_{l}$ and $P^{\prime}_{l}$ denote the derivatives of $\rho_{l}$ and $P_{l}$ with respect to the variable $z$, with $l=A,B$. We consider the case where the following condition holds

[TABLE]

which implies $P_{A}(z)\xrightarrow[z\rightarrow-\infty]{}P^{\rm M}$. Moreover, we note that the principle of mass conservation ensures that

[TABLE]

for some $M>0$. The results of our travelling-wave analysis are summarised in the following theorem.

Theorem 4.1**.**

Under Assumptions 3.1(i)-(iii), for any $M>0$ given there exist $c>0$ and $\ell>0$ such that the travelling-wave problem defined by system (4.9), complemented with the asymptotic condition (4.10), admits a solution whereby:

• (i)

$\rho_{A}(z)$* is decreasing in $(-\infty,0)$ and satisfies the condition*

[TABLE]

• (ii)

$\rho_{B}(z)$* is decreasing in $(0,\ell)$ and satisfies the condition*

[TABLE]

along with the condition (4.11).

Moreover, in the case where $F_{A}(\cdot)=F_{B}(\cdot)$, the following jump condition holds:

[TABLE]

Proof.

We prove Theorem 4.1 in five steps.

Step 1: existence of a solution of problem (4.9). For $c>0$ given, we have the following problem for $P_{B}(z)$

[TABLE]

Integrating the equation for $P_{B}$ over $(0,z)$, with $z<\ell$, and using the condition at $z=0$ we obtain an ordinary differential equation with final condition at $z=\ell$, that is,

[TABLE]

which can be solved explicitly giving P_{B}(z)=c\,(\ell-z)+F_{B}^{-1}\Big{(}\dfrac{\eta_{B}}{2}c\Big{)}. Notice that since $c>0$ we have $P_{B}(\ell)>P_{B}^{\rm eq}$ and $F_{B}$ is invertible. Knowing $P_{B}(z)$, we have the following problem for $P_{A}(z)$

[TABLE]

Integrating the equation for $P_{A}$ over $(z,0)$, with $z<0$, and using the second condition at $z=0$ we obtain an ordinary differential equation with final condition at $z=0$, that is,

[TABLE]

*Under Assumptions 3.1(i), (iii) on the functions $F_{A}$, $F_{B}$ and $G$, the above problem admits a solution. Hence, for a given $c>0$ there exists a solution of problem (4.9). *

Step 2: monotonicity of $\rho_{B}$ in $(0,\ell)$ and proof of the condition (4.13). Integrating  (4.9)2 between a generic point $z\in[0,\ell)$ and $\ell$, and using the boundary condition (4.9)7, we find

[TABLE]

Moreover, integrating (4.15) between a generic point $z\in[0,\ell)$ and $\ell$ gives

[TABLE]

Therefore,

[TABLE]

Furthermore, note that (4.15) allows us to rewrite the boundary condition (4.9)6 as

[TABLE]

Since under assumptions (4.4) we have that $F_{B}\left(P_{B}(\ell)\right)>0$ if and only if $P_{B}(\ell)>P^{\rm eq}_{B}$, we conclude that

[TABLE]

Finally, since the function $F_{B}\left(P_{B}\right)$ is monotone for $P_{B}>P^{\rm eq}_{B}$, cf. (4.4), the value of $P_{B}(\ell)$ is uniquely determined by (4.18).

Using the results (4.15), (4.17) and (4.19) along with the fact that, under assumptions (4.1) and (4.2), $P_{B}>0$ if and only if $\rho_{B}>\rho^{\rm eq}_{B}$ and $P_{B}$ is a monotonically increasing and continuous function of $\rho_{B}$ for $\rho_{B}>\rho^{\rm eq}_{B}$, we conclude that the function $\rho_{B}$ is continuous in $(0,\ell)$ and satisfies the following conditions

[TABLE]

and

[TABLE]

Step 3: identification of $\ell$. For $M>0$ given, since the value of $\rho_{B}(z)$ is uniquely determined for all $z\in[0,\ell]$, the value of $\ell$ is uniquely defined by the integral identity (4.11).

Step 4: monotonicity of $\rho_{A}$ in $(-\infty,0]$ and proof of the condition (4.12). Recalling that $\dfrac{dP_{A}(\rho_{A})}{d\rho_{A}}>0$ for $\rho_{A}>\rho^{\rm eq}_{A}$, we multiply both sides of (4.9)1 by $\dfrac{dP_{A}}{d\rho_{A}}$ and use assumption (4.10) to obtain the following boundary-value problem for $P_{A}$

[TABLE]

[TABLE]

Let $z^{*}\in(-\infty,0)$ be a critical point of $P_{A}$. Using (4.22) we conclude that

[TABLE]

Therefore, under the conditions (4.5) and (4.23), the strong maximum principle ensures that $P_{A}<P^{\rm M}$ in $(-\infty,0]$ and that $P_{A}$ cannot have a local minimum in $(-\infty,0)$, i.e.

[TABLE]

Hence the solution $P_{A}$ of (4.22), (4.23) is a continuous and nonincreasing function that satisfies

[TABLE]

Using the results (4.24) and (4.25) along with the fact that, under assumptions (4.1) and (4.2), $P_{A}>0$ if and only if $\rho_{A}>\rho^{\rm eq}$ and $P_{A}$ is a monotonically increasing and continuous function of $\rho_{A}$ for $\rho_{A}>\rho^{\rm eq}$, we conclude that the function $\rho_{A}$ is continuous in $(-\infty,0)$ and satisfies the following conditions

[TABLE]

and

[TABLE]

with $\rho^{\rm M}$ being related to $P^{\rm M}$ by (4.3).

Step 5: proof of the jump condition (4.14). The transmission conditions (4.9)4 and (4.9)5 give

[TABLE]

Furthermore, due to the uniqueness of the value of $P_{B}(0)>P^{\rm eq}_{B}$, under the monotonicity assumptions (4.4), the transmission condition (4.9)3 allows one to uniquely determine the value of $P_{A}(0)>P^{\rm eq}_{A}$. In particular, in the case where $F_{A}(\cdot)=F_{B}(\cdot)\equiv F(\cdot)$, the transmission condition (4.9)3 gives

[TABLE]

from which, using the monotonicity assumptions (4.1) and (4.4), one finds the jump condition (4.14). ∎

5. Numerical solutions of the free-boundary problem and computational simulations for the individual-based model

In this section, we illustrate the results established in Theorem 4.1 by presenting a sample of numerical solutions of the free-boundary problem (2.18), (2.19) and (2.22). Moreover, we compare numerical solutions of the continuum model with computational simulations for the individual-based model (2.3)-(2.1). We focus on the case where the force terms $F_{A}(d_{ij})$ and $F_{B}(d_{ij})$ are both given by the following cubic approximation of the JKR force law [27]

[TABLE]

where $d_{ij}=|r_{i}-r_{j}|$ and $d^{\rm eq}$ stands for the equilibrium intercellular distance, which is the distance between cell centres above which cells do not exert any force upon one another (i.e. $F(d_{ij})=0$ for all $d_{ij}\geq d^{\rm eq}$). The equilibrium distance $d^{\rm eq}$ and the coefficients $a_{1}$, $a_{2}$ and $a_{3}$ depend on the biophysical characteristics of the cells and are defined as

[TABLE]

In the formulas (5.2), $R$ is the cell radius, the parameter $\gamma$ measures the strength of cell-cell adhesion and $\tilde{E}$ is an effective Young’s modulus defined as

[TABLE]

with $E$ and $\nu$ being, respectively, the Young’s modulus and the Poisson’s ratio of the cells. We refer the interested reader to Appendix A for a detailed derivation of the approximate representation of the JKR force law given by (5.1)-(5.3).

Using the formal relations between the intercellular distance and the cell density (2.13) and (2.14), we compute the cell densities and the equilibrium cell density as

[TABLE]

and $\rho^{\rm eq}_{A}=\rho^{\rm eq}_{B}=\rho^{\rm eq}=1/{d^{\rm eq}}$. The approximation of the JKR force law (5.1) can be rewritten in terms of the cell densities $\rho_{l}$ and $\rho^{\rm eq}$ as follows

[TABLE]

Inserting the latter expression for $F\left(1/\rho_{l}\right)$ into the definition (2.20) of the nonlinear diffusion coefficient $D_{l}(\rho_{l})$ yields

[TABLE]

from which, using (4.1), we obtain the following barotropic relation for the cell pressure $P_{l}$

[TABLE]

In (5.7), the term $P^{l}_{0}$ is an integration constant such that $P_{l}(\rho^{\rm eq})=0$ and

[TABLE]

Numerical simulations were performed using parameter values chosen in agreement with those used in [20], that is,

[TABLE]

where $\zeta=10^{15}\,\text{m}^{2}$ is the density of cell-cell adhesion molecules in the cell membrane, $k_{B}$ is the Boltzmann constant and $T=298\,\text{K}$ is an absolute temperature. Figure 1 displays the plots of the force $F$ between neighbouring cells, the nonlinear diffusion coefficient $D_{l}$, and the cell pressure $P_{l}$, obtained using the parameter values given by (5.9).

We let the cell damping coefficients of population $A$ be $\eta_{A}=0.5\times 10^{-2}\,\text{kg}\,\text{s}^{-1}$, and considered the cases where $\eta_{A}=\eta_{B}$ or $\eta_{B}=2\,\eta_{A}$ or $\eta_{B}=0.5\,\eta_{A}$. Moreover, for the cell proliferation term we assumed

[TABLE]

where $\tilde{H}$ is a smooth approximation to the Heaviside function.

To construct numerical solutions, the free-boundary problem (2.18), (2.19) and (2.22) was transformed to a Lagrangian reference frame and the method of lines was employed to solve the resultant equations. The resulting system of ordinary differential equations, as well as the ordinary differential equations (2.4) and (2.1) of the individual-based model, were numerically solved using the Matlab routine ode15s.

The plots in Figures 2-4 show sample dynamics of the cell density $\rho$ defined as

[TABLE]

with $\rho_{A}$ and $\rho_{B}$ being either numerical solutions of the free-boundary problem (2.18), (2.19), (2.22) (black lines) or approximate cell densities computed from numerical solutions of the individual-based model (2.3)-(2.1) using (5.4) (red markers). In agreement with the results established in Theorem 4.1, we observe the emergence of travelling-wave solutions, whereby the positions of the inner free boundary $s_{1}(t)$ and the outer free boundary $s_{2}(t)$ move at the same constant speed, and the cell densities $\rho_{A}$ and $\rho_{B}$ are monotonically decreasing.

Moreover, $\rho$ is continuous at $s_{1}$ if $\eta_{A}=\eta_{B}$, cf. the plots in Figure 2, whereas it has a jump discontinuity at $s_{1}$ both for $\eta_{A}<\eta_{B}$ and for $\eta_{A}>\eta_{B}$. The sign of the jump $\rho(s_{1}^{+})-\rho(s_{1}^{-})$ satisfies condition (4.14), cf. the plots in Figures 3 and  4, and, once that the travelling-wave is formed, the size of the jump is constant and such that the transmission condition (4.9)3 is met, see Supplementary Figure B.1(a). As shown by these plots, there is an excellent match between the numerical solutions of the free-boundary problem and the computational simulation results for the individual-based model.

6. Discussion

We presented an off-lattice individual-based model that describes the dynamics of two contiguous cell populations with different proliferative and mechanical characteristics.

We formally showed that this discrete model can be represented in the continuum limit as a free-boundary problem for the cell densities. We proved an existence result for the free-boundary problem and constructed travelling-wave solutions. We performed numerical simulations in the case where the cellular interaction forces are described by the celebrated JKR model of elastic contact, and we found excellent agreement between the computational simulation results for the individual-based model, the numerical solutions of the corresponding free-boundary problem and the travelling-wave analysis. Taken together, the results of numerical simulations demonstrate that the solutions of the free-boundary problem faithfully capture the qualitative and quantitative properties of the outcomes of the off-lattice individual-based model.

In this paper, we focussed on a one-dimensional scenario where the two cell populations do not mix. It would be interesting to extend the individual-based mechanical model presented here, and the related formal method of derivation of the corresponding continuum model as well, to more realistic two-dimensional cases whereby spatial mixing between the two cell populations can occur. In this regard, an additional development of our study would be to formulate probabilistic discrete mechanical models of interacting cell populations and, using asymptotic methods analogous to those employed, for instance, in [11, 24, 38, 39], to perform a rigorous derivation of their continuum counterparts. These are all lines of research that we will be pursuing in the near future.

Appendix A Approximate representation of the JKR force law

The nonlinear function $F^{JKR}(d_{ij})$ that gives the JKR force law between the $i^{th}$ cell and the $j^{th}$ cell, with centres at distance $d_{ij}$, is implicitly defined by the following formulas [20]

[TABLE]

where

[TABLE]

In the formulas (A.1), the parameter $\gamma$ models the strength of cell-cell adhesion, $R_{i}$ stands for the radius of the $i^{th}$ cell, $E_{i}$ is the Young’s modulus of the $i^{th}$ cell, and $\nu_{i}$ denotes the Poisson’s ratio of the $i^{th}$ cell. Analogous considerations hold for the parameters of the $j^{th}$ cell. Moreover, $\delta_{ij}$ is the sum of the deformations undergone by the $i^{th}$ cell and the $j^{th}$ cell.

As the computational cost of numerical simulations carried out by solving implicitly for $F^{JKR}(d_{ij})$ is prohibitive, we derive an approximate representation of this function based on a third degree polynomial of the form

[TABLE]

with

[TABLE]

In the above equations, $d^{\rm eq}_{ij}$ denotes the equilibrium distance between the centres of cell $i$ and cell $j$ (i.e. the distance $d_{ij}$ such that $F(d_{ij})=0$ for all $d_{ij}\geq d^{\rm eq}_{ij}$).

With this goal in mind, we look for explicit expressions of $d^{\rm eq}_{ij}$, $F^{\prime}(d^{\rm eq}_{ij})$, $F^{\prime\prime}(d^{\rm eq}_{ij})$ and $F^{\prime\prime\prime}(d^{\rm eq}_{ij})$ in terms of the cell radii and the mechanical parameters of the cells. In the rest of this appendix we will use the abridged notation $F_{ij}$ for $F(d_{ij})$.

Expression for $d^{\rm eq}_{ij}$.

The equilibrium distance $d^{\rm eq}_{ij}$ can be directly computed from the formulas (A.1). In fact, choosing $d_{ij}=d^{\rm eq}_{ij}$ and using the fact that $F(d^{\rm eq}_{ij})=0$, from the second formula in (A.1) we obtain

[TABLE]

Substituting this expression into the first formula in (A.1) yields

[TABLE]

and noting that $d^{\rm eq}_{ij}=R_{i}+R_{j}-\delta^{\rm eq}_{ij}$ we find the equilibrium distance $d^{\rm eq}_{ij}$.

Expression for $F^{\prime}(d^{\rm eq}_{ij})$.

We substitute the second formula in (A.1) into the first formula to obtain

[TABLE]

where

[TABLE]

Differentiating $f$ with respect to $F_{ij}$ yields

[TABLE]

and

[TABLE]

Hence, for $d_{ij}=d^{\rm eq}_{ij}$ we have

[TABLE]

Differentiating both sides of (A.4) with respect to $d_{ij}$ we find

[TABLE]

where

[TABLE]

Rearranging terms in the latter equation yields

[TABLE]

and, therefore, for $d_{ij}=d^{\rm eq}_{ij}$ we have

[TABLE]

Finally, noting that

[TABLE]

we find

[TABLE]

Expression for $F^{\prime\prime}(d^{\rm eq}_{ij})$.

Differentiating twice both sides of (A.4) with respect to $d_{ij}$ yields

[TABLE]

Rearranging terms in the latter equation gives

[TABLE]

Hence, choosing $d_{ij}=d^{\rm eq}_{ij}$ and using the fact that $F(d^{\rm eq}_{ij})=0$ along with the expressions (A.5) for $f^{\prime}(0)$ and $f^{\prime\prime}(0)$, we find

[TABLE]

Noting that [cf. the expression of $f(0)$ in (A.5)]

[TABLE]

and inserting the expression (A.6) of $F^{\prime}(d^{\rm eq}_{ij})$ into (A.7) yields

[TABLE]

Expression for $F^{\prime\prime\prime}(d^{\rm eq}_{ij})$.

Differentiating thrice both sides of (A.4) with respect to $d_{ij}$ yields

[TABLE]

Rearranging terms in the latter equation we obtain

[TABLE]

Hence, choosing $d_{ij}=d^{\rm eq}_{ij}$ and using the fact that $F(d^{\rm eq}_{ij})=0$, along with the expression (A.7) for $F^{\prime\prime}(d^{\rm eq}_{ij})$ and the expressions (A.5) for $f^{\prime}(0)$, $f^{\prime\prime}(0)$ and $f^{\prime\prime\prime}(0)$, we find

[TABLE]

where

[TABLE]

Finally, inserting the expression (A.6) for $F^{\prime}(d^{\rm eq}_{ij})$ into (A.8) we obtain

[TABLE]

Expressions for $d^{\rm eq}_{ij}$, $a^{ij}_{1}$, $a^{ij}_{2}$ and $a^{ij}_{3}$.

Taken together the results from above give

[TABLE]

Approximate representation of $F^{JKR}(d_{ij})$ used to perform numerical simulations.

To perform numerical simulations, we assumed that all cells have the same radius $R$, Young’s modulus $E$ and Poisson’s ratio $\nu$, i.e.

[TABLE]

Under these assumptions, we have that

[TABLE]

and the approximate representation of the JKR force law given by (A.2) and (A.9) reads as

[TABLE]

with

[TABLE]

Appendix B Supplementary figure