Bridging the gap between individual-based and continuum models of growing cell populations

TL;DR

This paper establishes a formal link between individual-based stochastic models and continuum PDE models for growing cell populations, showing that the latter can be derived from the former and that simulations align well.

Contribution

It introduces a simple stochastic individual-based model and demonstrates how nonlinear PDEs for cell growth can be derived from it, bridging a gap in modeling approaches.

Findings

Derivation of PDEs from the individual-based model

Qualitative and quantitative agreement between models

Emergence of complex patterns from simple rules

Abstract

Continuum models for the spatial dynamics of growing cell populations have been widely used to investigate the mechanisms underpinning tissue development and tumour invasion. These models consist of nonlinear partial differential equations that describe the evolution of cellular densities in response to pressure gradients generated by population growth. Little prior work has explored the relation between such continuum models and related single-cell-based models. We present here a simple stochastic individual-based model for the spatial dynamics of multicellular systems whereby cells undergo pressure-driven movement and pressure-dependent proliferation.We show that nonlinear partial differential equations commonly used to model the spatial dynamics of growing cell populations can be formally derived from the branching random walk that underlies our discrete model. Moreover, we carry out a systematic comparison between the individual-based model and its continuum counterparts, both in the case of one single cell population and in the case of multiple cell populations with different biophysical properties. The outcomes of our comparative study demonstrate that the results of computational simulations of the individual-based model faithfully mirror the qualitative and quantitative properties of the solutions to the corresponding nonlinear partial differential equations. Ultimately, these results illustrate how the simple rules governing the dynamics of single cells in our individual-based model can lead to the emergence of complex spatial patterns of population growth observed in continuum models.

Figures

Equations

$${\partial }_t\rho -\mu \mathrm{div}\left(\rho \nabla p\right)=G(p)\rho ,$$

$$\frac{\mathrm{d}G}{\mathrm{d}p}<0,G(P)=0.$$

$$\Pi (0)=0,\frac{\mathrm{d}\Pi }{\mathrm{d}\rho }>0\text{ for }\rho >0.$$

$$\Pi (\rho )=K_{\gamma }{\rho }^{\gamma }\text{with}\gamma >1\text{and}{K}_{\gamma }>0.$$

$${\partial }_t{\rho }_h-{\mu }_h\mathrm{div}\left({\rho }_h\nabla p\right)=G_h(p){\rho }_h,h=1,\dots ,M.$$

$$\rho (t,x)=\sum _{h=1}^M{\rho }_h(t,x),$$

$$\left\{\begin{array}[]{ll}\partial_{t}\rho_{1}-\mu_{1}\,{\rm div}\left(\rho_{1}\nabla p\right)=G(p)\rho_{1},\\ \\ \partial_{t}\rho_{2}-\mu_{2}\,{\rm div}\left(\rho_{2}\nabla p\right)=0,\end{array}\right.$$

$${\rho }_h(t_k,x_i)={\rho }_{hi}^k=n_{hi}^k{\chi }^{-1}\text{ and }\rho (t_k,x_i)={\rho }_i^k=\sum _{h=1}^M{\rho }_h(t_k,x_i).$$

$$\tau G_h(p_i^k{)}_{+}\text{mbox}where{G}_h(p_i^k{)}_{+}=\max \left(0,G_h(p_i^k)\right)$$

$$\tau G_h(p_i^k{)}_{-}\text{mbox}where{G}_h(p_i^k{)}_{-}=-\min \left(0,G_h(p_i^k)\right)$$

$$1-\left(\tau G_h(p_i^k{)}_{+}+\tau G_h(p_i^k{)}_{-}\right)=1-\tau |G_h(p_i^k)|.$$

$$p_i^k\le \overline{p}\text{mbox}forall(k,i)\in {\mathbb{N}}_0\times \mathbb{Z},\text{ with }\overline{p}=\max \left(\underset{i\in \mathbb{Z}}{\max }{p}_i^0,P\right).$$

$$J_h^L(p_i^k-p_{i-1}^k)={\nu }_h\frac{(p_i^k-p_{i-1}^k{)}_{+}}{2\overline{p}},$$

$$J_h^R(p_i^k-p_{i+1}^k)={\nu }_h\frac{(p_i^k-p_{i+1}^k{)}_{+}}{2\overline{p}},$$

$$1-J_h^L(p_i^k-p_{i-1}^k)-J_h^R(p_i^k-p_{i+1}^k).$$

$${\rho }_{hi}^{k+1}$$

$${\rho }_{hi}^{k+1}-{\rho }_{hi}^k$$

$$t_k\approx t,t_{k+1}\approx t+\tau ,x_i\approx x,x_{i\pm 1}\approx x\pm \chi ,$$

$${\rho }_{hi}^k\approx {\rho }_h(t,x),{\rho }_{hi}^{k+1}\approx {\rho }_h(t+\tau ,x),{\rho }_{hi\pm 1}^k\approx {\rho }_h(t,x\pm \chi ),$$

$$p_i^k\approx p(t,x),p_{i\pm 1}^k\approx p(t,x\pm \chi ),$$

$${\rho }_h(t+\tau ,x)-{\rho }_h(t,x)$$

$${\rho }_h\in C^2({\mathbb{R}}^{+}\times \mathbb{R}),h=1,\dots ,M,$$

$$\tau {\partial }_t{\rho }_h(t,x)+\frac{{\tau }^2}{2}{\partial }_{tt}^2{\rho }_h(t,x)$$

$$F(t,x)$$

$${\partial }_t{\rho }_h(t,x)+\frac{\tau }{2}{\partial }_{tt}^2{\rho }_h(t,x)$$

$$\frac{{\nu }_h{\chi }^2}{2\overline{p}\tau }\to {\mu }_h\text{as }\tau \to 0\text{ and }\chi \to 0,\text{for}h=1,\dots ,M,$$

$${\partial }_t{\rho }_h=G_h(p){\rho }_h+\mu \left({\rho }_h{\partial }_{xx}^2{p}_h+{\partial }_x{\rho }_h{\partial }_x{p}_h\right),h=1,\dots ,M,$$

$${\partial }_t{\rho }_h-{\mu }_h{\partial }_x\left({\rho }_h{\partial }_{x}p\right)=G_h(p){\rho }_h,h=1,\dots ,M.$$

Full text

1 Introduction

Brief overview of continuum models of growing cell populations.

Continuum models for the spatial dynamics of growing cell populations have been widely used to complement empirical research in developmental biology and cancer research. These models consist of nonlinear partial differential equations that describe the evolution of cellular densities in response to pressure gradients generated by population growth, which can be mechanically regulated, nutrient-limited or pressure-dependent ambrosi2002mechanics ; ambrosi2002closure ; araujo2004history ; bresch2010computational ; byrne2010dissecting ; byrne1995growth ; byrne1996growth ; byrne1997free ; byrne2009individual ; byrne2003two ; byrne2003modelling ; chaplain2006mathematical ; chen2001influence ; ciarletta2011radial ; greenspan1976growth ; lowengrub2009nonlinear ; perthame2014some ; preziosi2003cancer ; ranft2010fluidization ; roose2007mathematical ; sherratt2001new ; ward1999mathematical ; ward1997mathematical .

For a population of cells with growth rate that depends on the local pressure bru2003universal ; byrne2003modelling ; drasdo2012modeling ; ranft2010fluidization , a prototypical example of such models is given by the following equation

[TABLE]

which was proposed by Byrne & Drasdo byrne2009individual . Equation (1) is a conservation equation for the function $\rho(t,x)\geq 0$ that represents the density of cells at position $x\in\mathbb{R}^{d}$, with $d=1,2,3$ in the biologically relevant cases, and time $t\in\mathbb{R}^{+}$. The function $p$ stands for the cell pressure and the term $G$ is the net growth rate of the cell density.

The second term on the left-hand side of equation (1) represents the rate of change of the cellular density due to pressure-driven cell movement (i.e. the tendency of cells to move towards regions where they feel less compressed). The definition of this term builds on the seminal paper of Greenspan greenspan1976growth and subsequent work of Byrne & Chaplain byrne1997free . In analogy with the classical Darcy’s law for fluids, the parameter $\mu>0$ can be seen as the cell mobility coefficient, which is defined as the quotient between the permeability of the porous medium in which the cells are embedded (e.g. the extracellular matrix) and the cellular viscosity.

The term on the right-hand side of equation (1) represents the rate of change of the cellular density due to cell proliferation (i.e. cell division and cell death). In the mathematical framework of equation (1), the effect of pressure-limited growth can be captured in a first approximation by assuming the net growth rate $G$ to be a smooth function of $p$ that satisfies the following assumptions

[TABLE]

In assumptions (2), the parameter $P>0$ stands for the pressure at which cell death exactly compensates cell division. The term homeostatic pressure has been coined to indicate such a critical pressure basan2009homeostatic .

In order to close equation (1) the pressure $p$ can be defined according to a barotropic relation $p\equiv\Pi(\rho)$. Typically, the function $\Pi(\rho)$ is identically zero for $\rho~{}\leq~{}\rho^{*}$ and is monotonically increasing for $\rho>\rho^{*}$, with $0<\rho^{*}<\Pi^{-1}(P)$ being the density below which cells do not exert any force upon one another tang2014composite . For simplicity, in this paper we let $\Pi$ be a smooth function of $\rho$ that satisfies the following assumptions

[TABLE]

For instance, in the attempt to reduce the biological problem to its essentials while ensuring analytical tractability of the mathematical model, Perthame et al. perthame2014hele have proposed the following definition of $\Pi(\rho)$, which satisfies assumptions (3):

[TABLE]

In definition (4), the parameter $\gamma$ provides a measure of the stiffness of the barotropic relation and $K_{\gamma}$ is a scale factor. The limit $\gamma\to\infty$ corresponds to the case where cells behave like an incompressible fluid. In this asymptotic regime, it has been proven that models of the form of equation (1) converge to free-boundary problems of Hele-Shaw type kim2018porous ; kim2016free ; mellet2017hele ; perthame2014derivation .

The model (1) can be generalised to the case of multiple cell populations with different biophysical properties (i.e. different mobilities and growth rates) as follows

[TABLE]

The system of coupled equations (5) relies on notation and assumptions analogous to those underlying equation (1). In particular, the coefficient $\mu_{h}>0$ represents the mobility of cells in the $h^{th}$ population and the pressure $p$ is defined by a barotropic relation $p\equiv\Pi(\rho)$. Here $\rho$ stands for the total cell density, i.e.

[TABLE]

and the function $\Pi$ satisfies conditions (3). Moreover, the net growth rate of the cell density $G_{h}(p)$ can be assumed to be a smooth function of the cell pressure that satisfies assumptions (2) for all $h~{}=~{}1,\dots,M$.

For example, building on the computational results presented by Drasdo & Hoehme drasdo2012modeling , Lorenzi et al. lorenzi2017interfaces considered the following variant of the system of equations (5)

[TABLE]

complemented with the barotropic relation (4). The system of equations (7) describes the interaction dynamics between a population of proliferating cells (i.e. population $1$) and a population of nonproliferating cells (i.e. population $2$) with different mobilities.

Further to the biological and clinical insights into the underpinnings of tissue development and tumour growth they can provide, these continuum models exhibit a range of interesting qualitative behaviours. For instance, it was shown that models in the form of equation (1) admit travelling-wave solutions with composite shapes and discontinuities tang2014composite . Moreover, in analogy with reaction-diffusion systems arising in the mathematical modelling of other biological and ecological problems dancer1999spatial ; mimura2000reaction , models like the system of equations (7) can give rise to sharp interfaces, which bring about spatial segregation between cell populations with different biophysical properties lorenzi2017interfaces .

Derivation of continuum models of growing cell populations from individual-based models.

A key advantage of continuum models for the spatial dynamics of growing cell populations over their individual-based counterparts (i.e. discrete models that track the dynamics of individual cells) drasdo2005coarse ; van2015simulating is that they are amenable to mathematical analysis. This enables a complete exploration of the model parameter space, which ultimately allows more robust conclusions to be drawn. Moreover, compared to individual-based models, continuum models offer the possibility to carry out numerical simulations at the level of larger portions of tissues or even of whole organs, while keeping computational costs within acceptable bounds.

Since continuum models are defined at the scale of whole cell populations, they are usually formulated on the basis of phenomenological considerations, which can hinder a precise mathematical description of crucial biological and physical aspects. On the contrary, stochastic individual-based models that describe the dynamics of single cells in terms of algorithmic rules can be more easily tailored to capture fine details of cellular dynamics, thus making it possible to achieve a more accurate mathematical representation of multicellular systems. Furthermore, individual-based models are able to reproduce the emergence of population-level phenomena that are induced by stochastic fluctuations in single-cell biophysical properties, which are relevant in the regime of low cellular densities and cannot easily be captured by continuum models. Therefore, it is desirable to derive continuum models for the spatial dynamics of cell populations as the appropriate limit of individual-based models for spatial cell movement and proliferation, in order to have a clearer picture of the modelling assumptions that are made and ensure that they correctly reflect the salient features of the underlying application problem.

For this reason, the derivation of continuum models formulated in terms of partial differential equations or partial integrodifferential equations from underlying individual-based models has attracted the attention of a considerable number of mathematicians and physicists. Examples in this active field of research include the derivation of continuum models of chemotaxis from velocity-jump process hillen2009user ; othmer1988models ; othmer2000diffusion ; painter2003modelling or from self-attracting reinforced random walks stevens2000derivation ; stevens1997aggregation ; the derivation of diffusion and nonlinear diffusion equations from underlying random walks champagnat2007invasion ; inoue1991derivation ; oelschlager1989derivation ; othmer2002diffusion ; penington2011building ; penington2014interacting , from systems of discrete equations of motion oelschlager1990large ; murray2009discrete ; murray2012classifying , from discrete lattice-based exclusion processes binder2009exclusion ; dyson2012macroscopic ; fernando2010nonlinear ; johnston2017co ; johnston2012mean ; landman2011myopic ; lushnikov2008macroscopic ; simpson2010cell or from cellular automata deroulers2009modeling ; drasdo2005coarse ; simpson2007simulating ; and, most recently, the derivation of nonlocal models of cell-cell adhesion from position-jump processes buttenschoen2018space . However, with the exception of the results presented in byrne2009individual ; engblom2018scalable ; motsch2018short , little prior work has investigated the relation between single-cell-based models and continuum models in the form of equation (1) and the system of equations (5). In particular, a derivation of these continuum models of growing cell populations from underlying individual-based models of spatial cell movement and proliferation has remained elusive.

Contents of the paper.

In the present paper we aim to bridge such a gap in the existing literature. In Section 2, we develop a simple stochastic individual-based model for the dynamics of growing cell populations. Our model is based on a branching random walk that takes into account the effects of pressure-driven cell movement and pressure-dependent cell proliferation. In Section 3, we show that equation (1) and the system of equations (5) can be formally derived from the branching random walk that underlies our discrete model. In Section 4, we carry out a systematic quantitative comparison between the individual-based model and its continuum counterparts, both in the case of one single cell population and in the case of multiple cell populations with different biophysical properties. In summary, we construct travelling-wave solutions both for equation (1) and for the system of equations (7) (Section 4.1), we present numerical solutions that illustrate the results of the travelling-wave analysis, and we compare such numerical solutions with the results of computational simulations of the individual-based model (Section 4.2). Section 5 concludes the paper and provides a brief overview of possible research perspectives.

2 An individual-based model for growing cell populations

We consider a multicellular system composed of $M$ cell populations. We represent each cell within the system as an agent that occupies a position on a lattice. Cells can move and proliferate according to a set of simple rules that result in a discrete-time branching random walk. For ease of presentation, we let cells be arranged along the real line $\mathbb{R}$, but there would be no additional difficulty in considering the case of branching random walks in higher spatial dimensions.

We discretise the time variable $t\in\mathbb{R}^{+}$ and the space variable $x\in\mathbb{R}$ as $t_{k}=k\tau$ with $k\in\mathbb{N}_{0}$ and $x_{i}=i\chi$ with $i\in\mathbb{Z}$, respectively, where $0<\tau,\chi\ll 1$. We introduce the dependent variable $n^{k}_{hi}\in\mathbb{N}_{0}$ to model the number of cells of population $h=1,\ldots,M$ on the lattice site $i$ and at the time-step $k$, and we compute the density of cells of population $h$ and the total cell density, respectively, as

[TABLE]

Moreover, for each lattice site $i$ and time-step $k$ we assume the cell pressure $p(t_{k},x_{i})=p^{k}_{i}$ to be given by a barotropic relation $p^{k}_{i}\equiv\Pi(\rho^{k}_{i})$ with $\Pi$ being a function of the total cell density that satisfies conditions (3). At each time-step, we allow every cell to undergo pressure-dependent proliferation and pressure-driven movement according to the following algorithmic rules, which are schematised in Figure 1.

Pressure-dependent cell proliferation.

We allow every cell to divide, die or remain quiescent with probabilities that depend on the local pressure, and we assume that a dividing cell is replaced by two identical daughter cells that are placed on the original lattice site of the parent cell. In order to model pressure-limited growth, given the net growth rate $G_{h}$ that satisfies assumptions (2) for all values of $h$, we assume that at the $k^{th}$ time-step a focal cell of population $h$ on the lattice site $i$ can divide with probability

[TABLE]

or die with probability

[TABLE]

or remain quiescent with probability

[TABLE]

We assume the time-step $\tau$ to be sufficiently small so that the quantities (9)-(11) are all between 0 and 1. Under assumptions (2), the definitions (9)-(11) are such that if $p^{k}_{i}>P$ then every cell on the $i^{th}$ lattice site can only die or remain quiescent at the $k^{th}$ time-step. Therefore, we have that

[TABLE]

Pressure-driven cell movement.

We model pressure-driven cell movement (i.e. the tendency of cells to move down pressure gradients) as a biased random walk whereby the movement probabilities depend on the difference between the pressure at the site occupied by a cell and the pressure at the neighbouring sites. In particular, for a focal cell of population $h$ on the lattice site $i$ at the time-step $k$, we define the probability of moving to the lattice site $i-1$ (i.e. the probability of moving left) as

[TABLE]

where $(p^{k}_{i}-p^{k}_{i-1})_{+}=\max\left(0,p^{k}_{i}-p^{k}_{i-1}\right)$, the probability of moving to the lattice site $i+1$ (i.e. the probability of moving right) as

[TABLE]

where $(p^{k}_{i}-p^{k}_{i+1})_{+}=\max\left(0,p^{k}_{i}-p^{k}_{i+1}\right)$, and the probability of remaining stationary on the lattice site $i$ as

[TABLE]

In the above equations, the coefficient $0<\nu_{h}\leq 1$ is directly proportional to the mobility of cells in population $h$ and the parameter $\overline{p}$ is defined in (12). Notice that the definitions (13)-(15) are such that the cells will move down pressure gradients. Moreover, the a priori estimate (12) ensures that the quantities defined according to (13)-(15) are all between 0 and 1.

3 Formal derivation of continuum models

In this section, we show how continuum models of growing cell populations of the form of equation (1) and of the system of equations (5) and (7) can be derived as formal limits of the branching random walk that underlies our stochastic individual-based model.

For a multicellular system the dynamic of which is governed by the algorithmic rules presented in Section 2, the principle of mass balance gives

[TABLE]

and after a little algebra we find

[TABLE]

The above equation simplifies to

[TABLE]

Using the fact that the following relations hold for $\tau$ and $\chi$ sufficiently small

[TABLE]

[TABLE]

[TABLE]

we rewrite equation (16) in the approximate form

[TABLE]

Assuming

[TABLE]

we approximate the terms $\rho_{h}(t+\tau,x)$, $\rho_{h}(t,x-\chi)$ and $\rho_{h}(t,x+\chi)$ in equation (17) by their second order Taylor expansions about the point $(t,x)$. Moreover, since $p\equiv\Pi(\rho)$ and $\Pi$ is a smooth function of $\rho$, under assumption (18) the pressure $p(t,x)$ is twice continuously differentiable with respect to the variable $x$. Hence we also approximate the terms $p(t,x-\chi)$ and $p(t,x+\chi)$ in equation (17) by their second order Taylor expansions about the point $(t,x)$. In so doing, after a little algebra we find

[TABLE]

which implies

[TABLE]

with

[TABLE]

Dividing both sides of the resulting equation by $\tau$ we obtain

[TABLE]

Letting both $\tau\to 0$ and $\chi\to 0$ in such a way that

[TABLE]

from (20) we formally obtain the following system of coupled conservation equations

[TABLE]

which can be rewritten as the system of equations (5), that is,

[TABLE]

In a similar way, in the case of one single cell population, letting $M=1$ and dropping the index $h$ we formally obtain equation (1). Moreover, we can formally obtain the system of equations (7) by choosing $M=2$, labelling the two populations by $h=1$ and $h=2$, and setting $G_{1}\equiv G$ and $G_{2}\equiv 0$.

Remark 1

We note that condition (21) is a natural counterpart of the usual parabolic scaling of Brownian motion. Hence, our formal derivation does not impose any additional assumptions than those commonly employed for the asymptotic investigation of random walks.

4 Comparison between individual-based and continuum models

In this section, we carry out a systematic quantitative comparison between the outcomes of our individual-based model, both in the case of one cell population and in the case of two cell populations, and the solutions of the corresponding continuum models. We first establish the existence of travelling-wave solutions for the continuum models (1) and (7) (Section 4.1). We then construct numerical solutions of the model equations which illustrate the results of the travelling-wave analysis and we compare such numerical solutions with the results of computational simulations of our individual-based model (Section 4.2).

4.1 Travelling-wave analysis of the continuum models

We first consider the continuum model (1) and we look for one-dimensional travelling-wave solutions of the form

[TABLE]

that satisfy the following asymptotic conditions

[TABLE]

Therefore, we study the existence of pairs $(\rho,c)$ that satisfy the problem defined by the differential equation

[TABLE]

subject to conditions (23). Our main results are summarised by the following theorem.

Theorem 4.1

Under assumptions (2) and (3), there exists $c>0$ such that the travelling-wave problem defined by the differential equation (24) subject to conditions (23) admits a nonnegative and nonincreasing solution $\rho$.

Proof

We divide the proof of Theorem 4.1 into two steps. First we prove that for $c>0$ fixed the differential equation (24) subject to the asymptotic conditions (23) admits a nonnegative and nonincreasing solution $\rho$ (Step 1). Then we show that there exists a unique value of the wave speed $c$ that satisfies the travelling-wave problem (Step 2).

Step 1. Multiplying both sides of the differential equation (24) by $\displaystyle{\frac{{\mathrm{d}}p}{{\mathrm{d}}\rho}}$ we obtain the following boundary-value problem for $p$

[TABLE]

[TABLE]

Let $z^{*}$ be a critical point of $p$ in $\mathbb{R}$. Using the differential equation (25) we see that

[TABLE]

and, under assumptions (2) and conditions (26), using the strong maximum principle we conclude that $p<P$ in $\mathbb{R}$ and that $p$ cannot have a local minimum in $\mathbb{R}$, i.e.

[TABLE]

Hence the solution $p$ of the differential equation (25) subject to conditions (26) is a nonnegative and nonincreasing function that satisfies

[TABLE]

Since $p\equiv\Pi(\rho)$ and $\Pi$ is a smooth monotonically increasing function of $\rho$, we can then conclude that the cell density $\rho$ is a nonnegative and nonincreasing function as well.

Step 2. Using a method similar to that used in Step 5 of the proof of Theorem 4.2 one can prove that $p$ is a monotonically decreasing function of the parameter $c$. This ensures that given the pressure $p$ or, equivalently, the cell density $\rho$ the wave speed $c$ can be uniquely identified through a monotonicity argument. ∎

We then turn to the travelling-wave analysis of the system of equations (7). On the basis of the results presented in lorenzi2017interfaces , we consider one-dimensional travelling-wave solutions of the form

[TABLE]

that satisfy the following conditions

[TABLE]

for some $\ell>0$, along with the asymptotic condition

[TABLE]

Hence, we study the existence of triples $(\rho_{1},\rho_{2},c)$ along with $\ell>0$ that satisfy the problem defined by the following system of differential equations

[TABLE]

subject to conditions (29) and (30). Notice that the principle of mass conservation gives

[TABLE]

for some $N_{2}>0$ that represents the number of cells in population $2$. Our main results are summarised by the following theorem.

Theorem 4.2

Under assumptions (2) and (3), for any $N_{2}>0$ given, there exists $c>0$ and $\ell>0$ such that the system of differential equations (31) subject to conditions (29) and (30) admits a component-wise nonnegative solution $(\rho_{1},\rho_{2})$ with $\rho_{1}$ nonincreasing, and $\rho_{2}$ nonincreasing and satisfying condition (32). Moreover, the pressure $p$ has a kink in $z=0$ with

[TABLE]

Proof

Building upon the method of proof presented by Lorenzi et al lorenzi2017interfaces for the case of the barotropic relation (4), we prove Theorem 4.2 in five steps. We fix the parameter $c>0$ and first prove that, for $N_{2}>0$ given, the problem under study admits a component-wise nonnegative solution $(\rho_{1},\rho_{2})$ with $\rho_{2}$ nonincreasing and with the value of $\ell$ being determined by condition (32) (Step 1), and with $\rho_{1}$ nonincreasing (Step 2). Then we prove that the total cell density $\rho$ is continuous on $(-\infty,\ell)$ (Step 3) and the jump condition (33) holds (Step 4). Finally, we show that there exists a unique value of the wave speed $c$ that satisfies the travelling-wave problem (Step 5).

Step 1. Integrating the differential equation (31)2 between a generic point $z\in[0,\ell]$ and $\infty$, and using the fact that both $p^{\prime}(z)\to 0$ and $\rho_{2}(z)\to 0$ as $z\to\infty$, we find

[TABLE]

Integrating equation (34) between a generic point $z\in[0,\ell)$ and $\ell$, and using the fact that $p(\ell)=0$, gives

[TABLE]

which implies that

[TABLE]

Since $\rho_{1}\equiv 0$ on $[0,\ell]$, under assumptions (3), we have that $p$ is a monotonically decreasing function of $\rho_{1}$ in $[0,\ell]$. Hence, the results (34) and (35) allow us to conclude that $\rho_{2}$ is decreasing in $(0,\ell)$. Moreover, for $N_{2}>0$ given, since the value of $\rho_{2}(z)$ is uniquely determined for all $z\in[0,\ell]$, the value of $\ell$ is uniquely fixed by the mass conservation condition (32).

Step 2. Since $\rho_{2}\equiv 0$ on $(-\infty,0)$ and, therefore, $\rho_{1}\equiv\rho$ on $(-\infty,0)$, multiplying both sides of the differential equation (31)1 by $\displaystyle{\frac{{\mathrm{d}}p}{{\mathrm{d}}\rho}}$ we obtain the following boundary-value problem for $p$

[TABLE]

[TABLE]

Hence, using a method similar to that used in Step 1 of the proof of Theorem 4.1 one can prove that $\rho_{1}$ is decreasing in $(-\infty,0)$.

Step 3. The results proved in Step 1 and Step 2 ensure that the total cell density $\rho$ is nonincreasing and continuous in $(-\infty,0)$ and in $(0,\ell)$. We now prove that $\rho$ is continuous in $z=0$. Adding together the differential equations (31)1 and (31)2 gives

[TABLE]

Multiplying both sides of the above differential equation by $p$ and using the fact that

[TABLE]

we find that

[TABLE]

Integrating both sides of the latter differential equation between a generic point $z^{*}<0$ and $\ell$, and estimating the right-hand side from above by using the fact that $-\infty<\left(\rho_{1}+\rho_{2}\right)^{\prime}(z)\leq 0$ for all $z\in[z^{*},\ell)$, $p(\ell)=0$ and $\rho_{2}(z^{*})=0$, yields

[TABLE]

The above integral inequality ensures that $p^{\prime}\in L^{2}_{loc}(\mathbb{R})$. This result along with the fact that $p\in L^{\infty}(\mathbb{R})$ allow us to conclude that $p$ is continuous in $z=0$. Since $p\equiv\Pi(\rho)$ and $\Pi$ is a smooth monotonically increasing function of $\rho$, we have that the total cell density $\rho$ is continuous in $z=0$ as well.

Step 4. Integrating the differential equation (39) between a generic point $z<\ell$ and $\ell$ and using the fact that $\rho_{1}(\ell)=\rho_{2}(\ell)=0$ yields

[TABLE]

Letting $z\to 0^{-}$ and using the fact that $\rho_{1}(0^{-})=\rho(0^{-})$ and $\rho_{1}\equiv 0$ on $[0,\ell]$ we find that

[TABLE]

Similarly, letting $z\to 0^{+}$ and using the fact that $\rho_{2}(0^{+})=\rho(0^{+})$ and $\rho_{1}\equiv 0$ on $[0,\ell]$ gives

[TABLE]

Since $\rho(z)$ is continuous in $z=0$, combining equations (40) and (41) we obtain

[TABLE]

This result along with the expression (34) for $p^{\prime}(0^{+})$ gives

[TABLE]

From the latter equation we deduce the jump condition (33). Moreover, substituting the expression of $p^{\prime}(0^{-})$ into the differential equation (37) gives $p^{\prime\prime}(0^{-})<0$.

Step 5. We prove that $p$ is a monotonically decreasing function of the parameter $c$ on $(-\infty,0)$. This ensures that given the cell density $\rho_{1}$ the wave speed $c$ can be uniquely identified through a monotonicity argument. Recalling that $\rho_{2}\equiv 0$ on $(-\infty,0)$ and, therefore, $\rho_{1}\equiv\rho$ on $(-\infty,0)$, differentiating equation (31)1 with respect to $z$ we find

[TABLE]

with

[TABLE]

On the other hand, differentiating equation (31)1 with respect to $c$ gives

[TABLE]

with

[TABLE]

Using the fact that $\displaystyle{p^{\prime}=\frac{{\mathrm{d}}p}{{\mathrm{d}}\rho}\,\rho^{\prime}}$, we rewrite the differential equations (42) and (43), respectively, as

[TABLE]

and

[TABLE]

Since $\displaystyle{\frac{{\mathrm{d}}p}{{\mathrm{d}}\rho}>0}$, we have

[TABLE]

Hence, introducing the notation $f=p^{\prime}$, $\displaystyle{g=\frac{\partial p}{\partial c}}$ and

[TABLE]

we rewrite the differential equations (44) and (45), respectively, as

[TABLE]

with

[TABLE]

and

[TABLE]

with

[TABLE]

Since $f(z)<0$ for all $z\in(-\infty,0)$ and $f^{\prime}(0)<0$, noting that both $f(z)\to 0$ and $g(z)\to 0$ as $z\to-\infty$ and the right-hand side of the differential equation (48) contains the additional term $\rho^{\prime}<0$ compared to the right-hand side of the differential equation (46), we deduce that $\displaystyle{g=\frac{\partial p}{\partial c}<0}$, which concludes the proof of Theorem 4.2. ∎

Remark 2

Based on the jump condition (33), we expect the travelling-wave solution of Theorem 4.2 to be unstable if $\mu_{1}>\mu_{2}$. In fact, a small perturbation of $\rho_{1}(z)$ that is greater than zero on $[0,\ell)$ will propagate with approximate speed $-\mu_{1}\,p^{\prime}(0^{+})$. Noting that when $\mu_{1}>\mu_{2}$ the jump condition (33) gives $-\mu_{1}\,p^{\prime}(0^{+})>-\mu_{1}\,p^{\prime}(0^{-})$, we deduce that such a perturbation will separate from the rest of the travelling wave $\rho_{1}(z)$.

4.2 Quantitative comparison between individual-based and continuum models

4.2.1 One cell population

We begin by comparing the computational simulation results of our individual-based model in the case of one cell population with the numerical solutions of equation (1). For consistency with equation (1), we set $M=1$ and we drop the index $h=1$ both from the functions and from the parameters of the individual-based model. A complete description of the setup of numerical simulations is given in Appendix A.1 and Appendix B.1. In particular, we use the following definition of the net growth rate $G$

[TABLE]

so that assumptions (2) are satisfied.

Travelling fronts.

Figure 2, along with the video accompanying it (vid. Online Resource 1), demonstrates that there is an excellent quantitative match between the numerical solutions of equation (1) and the computational simulation results of our individual-based model. In agreement with the results established by Theorem 4.1, the cell density is nonincreasing and connects the homogeneous steady state $\rho\equiv\Pi^{-1}(P)$ to the homogeneous steady state $\rho\equiv 0$. Accordingly, the cell pressure is nonincreasing and connects the homogeneous steady state $p\equiv P$ to the homogeneous steady state $p\equiv 0$.

Higher values of $\beta$ lead to higher speed of invasion.

Figure 3, along with the video accompanying it (vid. Online Resource 2), indicates that, as one would expect, increasing the value of the parameter $\beta$ in the definition (50) of $G(p)$ accelerates the growth of the cell population, thus leading to a higher speed of invasion.

Differences between the outcomes of individual-based and continuum models in the presence of sharp transitions from high to low cell densities.

The results presented so far indicate that there is an excellent agreement between the computational simulation results of our individual-based model and the solutions of the corresponding continuum models. However, due to extinction phenomena related to stochasticity effect that occur in the individual-based model for low cell densities, we expect differences between the outcomes of the two modelling approaches to emerge in the presence of sharp transitions from high to low total cell densities. In order to verify this hypothesis, exploiting the asymptotic results of Perthame et al.  perthame2014hele – who have shown that under the barotropic relation (4) sufficiently high values of the parameter $\gamma$ lead equation (1) to develop sharper invasion fronts – we compare the computational simulation results of our individual-based model in the case of one cell population with the numerical solutions of equation (1) under the barotropic relation (4) for increasing values of $\gamma$. A complete description of the setup of numerical simulations is given in Appendix A.1 and Appendix B.1. The results obtained are summarised by Figure 4 which shows that larger values of the parameter $\gamma$ can bring about sharper invasion fronts, which come along with more abrupt variations in the cell density, thus leading to more evident differences between the computational simulation results of the stochastic individual-based model and the numerical solutions of equation (1) at the front of invasion. Ultimately, this causes the invasion front of the individual-based model to travel at the same speed but behind the front of the corresponding continuum model.

4.2.2 Two cell populations

We now turn to the case of two cell populations and we compare computational simulation results of our individual-based model with numerical solutions of the system of equations (7). For consistency with the system of equations (7), we choose $M=2$, and we set $G_{1}\equiv G$ and $G_{2}\equiv 0$ in the individual-based model. Full details of the setup of numerical simulations can be found in Appendix A.2 and Appendix B.2 for the results reported in Figures 5 and 6, and Appendix A.3 and Appendix B.3 for the results reported in Figures 7 and 8. In particular, we use the definition (50) of the net growth rate $G(p)$ with the homeostatic pressure $P=10\times 10^{4}$ and we let the population of nonproliferating cells (i.e. population $2$) be ahead of the population of proliferating cells (i.e. population $1$) at the initial time $t=0$.

Figures 5 and 6, along with the videos accompanying them (vid. Online Resource 3 and Online Resource 4), demonstrate that there is an excellent quantitative match between the numerical solutions of the system of equations (7) and the computational simulation results of our individual-based model, both in the case where $\mu_{1}\leq\mu_{2}$ and when $\mu_{1}>\mu_{2}$. Over time, the pressure $p$ converges to the homeostatic pressure $P$ while the cell density $\rho_{1}$ converges to the corresponding value $\Pi^{-1}(P)$.

Travelling fronts and spatial segregation between the two cell populations.

If $\mu_{1}\leq\mu_{2}$ (i.e. if $\nu_{1}\leq\nu_{2}$), in agreement with the results established by Theorem 4.2, spatial segregation occurs and the two cell populations remain separated by a sharp interface (vid. Figure 5 and Online Resource 3). The population of nonproliferating cells stays ahead of the population of proliferating cells and, over the regions where they are greater than zero, the cell densities are nonincreasing. The pressure itself is continuous across the interface between the two cell populations, whereas its first derivative jumps from a lower negative value to a larger negative value, i.e. the sign of the jump coincides with ${\rm sgn}(\mu_{2}-\mu_{1})$ (cf. the jump condition (33)).

Mixing between the two cell populations.

If $\mu_{1}>\mu_{2}$ (i.e. if $\nu_{1}>\nu_{2}$) the cell population 2 is left behind by the cell population 1, which ultimately propagates alone (vid. Figure 6 and Online Resource 4). This is consistent with the heuristic argument provided in Remark 2, which suggests that the travelling-wave solutions of Theorem 4.2 are unstable in the case where $\mu_{1}>\mu_{2}$.

Numerical simulations for more realistic barotropic relations.

In this paper, we have focused on barotropic relations that satisfy assumptions (3). However, as mentioned earlier in Section 1, more realistic barotropic relations satisfy the following conditions

[TABLE]

with $0<\rho^{*}<\Pi^{-1}(P)$. Figures 7 and 8, along with the videos accompanying them (vid. Online Resource 5 and Online Resource 6), display the computational simulation results of our stochastic individual-based model and the numerical solutions of the system of equations (7) obtained under a barotropic relation that satisfies the more general assumptions (51) (vid. Appendix A.3 and Appendix B.3 for a complete description of the numerical simulation setup). These computational results and numerical solutions clearly share the same properties as those of Figures 5 and 6. This supports the conclusion that the essentials of the results obtained using barotropic relations of the form (3) remain intact even under the more realistic assumptions (51).

5 Conclusions and research perspectives

In summary, we have developed a simple, yet effective, stochastic individual-based model for the spatial dynamics of multicellular systems whereby cells undergo pressure-driven movement and pressure-dependent proliferation. We have shown that nonlinear partial differential equations commonly used to model the spatial dynamics of growing cell populations can be formally derived from the branching random walk that underlies our discrete model. Moreover, we have carried out a systematic comparison between the individual-based model and its continuum counterparts, both in the case of one single cell population and in the case of multiple cell populations with different biophysical properties. The outcomes of our comparative study demonstrate that the results of computational simulations of the individual-based model faithfully mirror the qualitative and quantitative properties of the solutions to the corresponding nonlinear partial differential equations. Ultimately, these results illustrate how the simple rules governing the dynamics of single cells in our individual-based model can lead to the emergence of complex spatial patterns of population growth observed in continuum models – e.g. travelling waves with composite shapes and sharp interfaces corresponding to spatial segregation between cell populations with different biophysical properties.

We have focussed on the case of one spatial dimension and considered barotropic relations that satisfy assumptions (3). However, the model presented here, as well as the related formal method to derive corresponding continuum models, can be adapted to higher spatial dimensions and more realistic barotropic relations of the form (51). In this regard, it would be interesting to use our stochastic individual-based model to further investigate the formation of finger-like patterns of invasion observed for the system of equations (7) posed on a two dimensional spatial domain lorenzi2017interfaces . Such spatial patterns resemble infiltrating patterns of cancer-cell invasion commonly observed in breast tumours wang2012adipose . An additional development of our study would be to compare the results presented here with those obtained from equivalent models defined on irregular lattices, as well as to investigate how our modelling approach could be related to off-lattice models of growing cell populations drasdo2005coarse ; motsch2018short ; van2015simulating .

Our modelling framework can be easily extended to incorporate additional layers of biological complexity, such as nutrient-limited proliferation, undirected random cell movement, chemotaxis and haptotaxis. These are all lines of research that we will be pursuing in the near future.

Appendix A Details of numerical simulations of the individual-based model

We use a uniform discretisation of the interval $[0,100]$ that consists of 1001 points as the spatial domain (i.e. the grid-step is $\chi=0.1$) and we choose the time-step $\tau=2\times 10^{-3}$. We implement zero-flux boundary conditions by letting the attempted move of a cell be aborted if it requires moving out of the spatial domain. For all simulations, we use the definition (50) of $G(p)$ and we perform numerical computations in Matlab. Further details specific either to the case of one single cell population or to the case of two cell populations are provided in the next subsections.

A.1 Setup of numerical simulations for the case of one single cell population

For consistency with equation (1), we set $M=1$ and we drop the index $h=1$ both from the functions and from the parameters of the individual-based model. We define the rate $G$ in equations (9)-(11) according to equation (50). We set the homeostatic pressure $P=120\times 10^{4}$ for the simulation results reported in Figure 2 and Figure 3, while we choose $P=120\times 10^{5}$ for the simulation results of Figure 4. Moreover, we choose $\beta=4\times 10^{-6}$ for the simulation results reported in Figure 2, $\beta=4\times 10^{-5}$ for the simulation results of Figure 4, and $\beta\in\left\{1.5\times 10^{-6},4\times 10^{-6},4\times 10^{-5}\right\}$ for the simulation results of Figure 3. We set $\nu=0.02$ in equations (13)-(15) and we define the pressure $p^{k}_{i}$ according to the following barotropic relation

[TABLE]

which satisfies conditions (3). We let $\gamma\in\left\{1.2,1.5,2\right\}$ for the simulation results of Figure 4, while we choose $\gamma=1.2$ for the simulation results reported in Figure 2 and Figure 3. We use the initial cell density

[TABLE]

The results presented in Figure 2 and Figure 3 correspond to the average over three simulations of our individual-based model, while the results in Figure 4 correspond to one single simulation when $\gamma=1.5$ or $\gamma=2$ and the average over two simulations when $\gamma=1.2$.

A.2 Setup of numerical simulations for the case of two cell populations – Figure 5 and Figure 6

For consistency with the system of equations (7), we choose $M=2$, and we set $G_{1}\equiv G$ and $G_{2}\equiv 0$ in equations (9)-(11), where $G$ is defined according to equation (50) with the homeostatic pressure $P=10\times 10^{4}$ and the factor $\beta=4\times 10^{-5}$. We set $\nu_{1}=0.01$ and $\nu_{2}=0.5$ in equations (13)-(15) for the simulation results reported in Figure 5, while we consider $\nu_{1}=0.5$ and $\nu_{2}=0.01$ for the simulation results of Figure 6. We define the pressure $p^{k}_{i}$ according to the following simplified barotropic relation

[TABLE]

which satisfies conditions (3). We make use of the initial cell densities

[TABLE]

where

[TABLE]

The results presented in Figure 5 and Figure 6 correspond to one single simulation of our individual-based model.

A.3 Setup of numerical simulations for the case of two populations – Figure 7 and Figure 8

For consistency with the system of equations (7), we choose $M=2$, and we set $G_{1}\equiv G$ and $G_{2}\equiv 0$ in equations (9)-(11), where $G$ is defined according to equation (50) with the homeostatic pressure $P=10\times 10^{4}$ and the factor $\beta=4\times 10^{-5}$. We set $\nu_{1}=0.01$ and $\nu_{2}=0.5$ in equations (13)-(15) for the simulation results reported in Figure 7, while we consider $\nu_{1}=0.5$ and $\nu_{2}=0.01$ for the simulation results of Figure 8. We define the pressure $p^{k}_{i}$ according to the following barotropic relation

[TABLE]

which satisfies conditions (51). We make use of the initial cell densities (52) with

[TABLE]

The results presented in Figure 7 and Figure 8 correspond to one single simulation of our individual-based model.

Appendix B Details of numerical simulations of the continuum models

We let $x\in[0,100]$ and we construct numerical solutions for equation (1) and for the system of equations (7) complemented with zero Neumann boundary conditions. We use a finite volume method based on a time-splitting between the conservative and nonconservative parts. For the conservative parts, transport terms are approximated through an upwind scheme whereby the cell edge states are calculated by means of a high-order extrapolation procedure leveque2002finite , while the forward Euler method is used to approximate the nonconservative parts. We consider a uniform discretisation of the interval $[0,100]$ that consists of 1001 points and we perform numerical computations in Matlab. For all simulations, we use the definition (50) of $G(p)$. Further details specific either to the case of one cell population or to the case of two cell populations are provided in the next subsections.

B.1 Setup of numerical simulations for equation (1)

The rate $G$ is defined according to equation (50) with the homeostatic pressure $P=120\times 10^{4}$ for the numerical solutions reported in Figure 2 and Figure 3, while we choose $P=120\times 10^{5}$ for the numerical solutions of Figure 4. Moreover, we choose $\beta=4\times 10^{-6}$ for the numerical solutions reported in Figure 2, $\beta=4\times 10^{-5}$ for the numerical solutions of Figure 4, and $\beta\in\left\{1.5\times 10^{-6},4\times 10^{-6},4\times 10^{-5}\right\}$ for the numerical solutions reported in Figure 3. We define the pressure $p$ according to the following barotropic relation

[TABLE]

which satisfies conditions (3). We let $\gamma\in\left\{1.2,1.5,2\right\}$ for the numerical solutions of Figure 4, while we choose $\gamma=1.2$ for the numerical solutions reported in Figure 2 and Figure 3. Given the parameter values used for the individual-based model in the case of one single cell population, we choose the mobility $\mu=4.166\times 10^{-7}$ for the numerical solutions reported in Figure 2 and Figure 3, while we set $\mu=4.166\times 10^{-8}$ for the numerical solutions of Figure 4. In this way, both values of $\mu$ satisfy condition (21) for $h=1$. We impose the initial condition

[TABLE]

B.2 Setup of numerical simulations for the system of equations (7) – Figure 5 and Figure 6

The rate $G$ is defined according to equation (50) with the homeostatic pressure $P=10\times 10^{4}$ and $\beta=4\times 10^{-5}$. Given the parameter values used for the individual-based model in the case of two cell populations, we choose the mobilities $\mu_{1}=2.5\times 10^{-7}$ and $\mu_{2}=1.25\times 10^{-5}$ for the numerical solutions reported in Figure 5, and the mobilities $\mu_{1}=1.25\times 10^{-5}$ and $\mu_{2}=2.5\times 10^{-7}$ for the numerical solutions of Figure 6. This ensures that conditions (21) for $h=1,2$ are satisfied. We define the pressure $p$ according to the following simplified barotropic relation

[TABLE]

which satisfies conditions (3). We impose the initial conditions

[TABLE]

where

[TABLE]

B.3 Setup of numerical simulations for the system of equations (7) – Figure 7 and Figure 8

The rate $G$ is defined according to equation (50) with the homeostatic pressure $P=10\times 10^{4}$ and $\beta=4\times 10^{-5}$. Given the parameter values used for the individual-based model in the case of two cell populations, we choose the mobilities $\mu_{1}=2.5\times 10^{-7}$ and $\mu_{2}=1.25\times 10^{-5}$ for the numerical solutions reported in Figure 7, and the mobilities $\mu_{1}=1.25\times 10^{-5}$ and $\mu_{2}=2.5\times 10^{-7}$ for the numerical solutions of Figure 8. This ensures that conditions (21) for $h=1,2$ are satisfied. We define the pressure $p$ according to the following barotropic relation

[TABLE]

which satisfies conditions (51). We impose the initial conditions (53) with

[TABLE]

Acknowledgements.

FRM is funded by the Engineering and Physical Sciences Research Council (EPSRC). TL and FRM gratefully acknowledge Dirk Drasdo and Luís Neves de Almeida for insightful discussions.