Open Access

Research articles

Discrete and continuous mathematical models of sharp-fronted collective cell migration and invasion

School of Mathematical Sciences, Queensland University of Technology, Brisbane, Queensland, Australia

Contribution: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Visualization, Writing – original draft

Google Scholar

Find this author on PubMed

Search for more papers by this author

Keeley M. Murphy

School of Mathematical Sciences, Queensland University of Technology, Brisbane, Queensland, Australia

Contribution: Conceptualization, Formal analysis, Investigation, Methodology, Writing – review and editing

Google Scholar

Find this author on PubMed

Search for more papers by this author

Scott W. McCue

School of Mathematical Sciences, Queensland University of Technology, Brisbane, Queensland, Australia

Contribution: Conceptualization, Supervision, Writing – review and editing

Google Scholar

Find this author on PubMed

Search for more papers by this author

and

Pascal R. Buenzli

https://orcid.org/0000-0003-3962-5393

School of Mathematical Sciences, Queensland University of Technology, Brisbane, Queensland, Australia

Contribution: Conceptualization, Supervision, Writing – review and editing

Google Scholar

Find this author on PubMed

Search for more papers by this author

Matthew J. Simpson

https://orcid.org/0000-0001-6254-313X

School of Mathematical Sciences, Queensland University of Technology, Brisbane, Queensland, Australia

[email protected]

Google Scholar

Find this author on PubMed

Search for more papers by this author

Keeley M. Murphy

School of Mathematical Sciences, Queensland University of Technology, Brisbane, Queensland, Australia

Contribution: Conceptualization, Formal analysis, Investigation, Methodology, Writing – review and editing

Google Scholar

Find this author on PubMed

Search for more papers by this author

Scott W. McCue

School of Mathematical Sciences, Queensland University of Technology, Brisbane, Queensland, Australia

Contribution: Conceptualization, Supervision, Writing – review and editing

Google Scholar

Find this author on PubMed

Search for more papers by this author

and

Pascal R. Buenzli

https://orcid.org/0000-0003-3962-5393

School of Mathematical Sciences, Queensland University of Technology, Brisbane, Queensland, Australia

Contribution: Conceptualization, Supervision, Writing – review and editing

Google Scholar

Find this author on PubMed

Search for more papers by this author

Published:15 May 2024https://doi.org/10.1098/rsos.240126

Review history
- Review history is available via Web of Science at https://www.webofscience.com/api/gateway/wos/peer-review/10.1098/rsos.240126

Abstract

Mathematical models describing the spatial spreading and invasion of populations of biological cells are often developed in a continuum modelling framework using reaction–diffusion equations. While continuum models based on linear diffusion are routinely employed and known to capture key experimental observations, linear diffusion fails to predict well-defined sharp fronts that are often observed experimentally. This observation has motivated the use of nonlinear degenerate diffusion; however, these nonlinear models and the associated parameters lack a clear biological motivation and interpretation. Here, we take a different approach by developing a stochastic discrete lattice-based model incorporating biologically inspired mechanisms and then deriving the reaction–diffusion continuum limit. Inspired by experimental observations, agents in the simulation deposit extracellular material, which we call a substrate, locally onto the lattice, and the motility of agents is taken to be proportional to the substrate density. Discrete simulations that mimic a two-dimensional circular barrier assay illustrate how the discrete model supports both smooth and sharp-fronted density profiles depending on the rate of substrate deposition. Coarse-graining the discrete model leads to a novel partial differential equation (PDE) model whose solution accurately approximates averaged data from the discrete model. The new discrete model and PDE approximation provide a simple, biologically motivated framework for modelling the spreading, growth and invasion of cell populations with well-defined sharp fronts. Open-source Julia code to replicate all results in this work is available on GitHub.

1. Introduction

Continuum partial differential equation (PDE) models have been used for over 40 years to model and interpret the spatial spreading, growth and invasion of populations of cells [1–3]. PDE models have been used to improve our understanding of various biological processes including wound healing [4–9], embryonic development [10–13], tissue growth [14–16] as well as disease progression, such as cancer [17–24]. For a homogeneous population of cells with density $u \geq 0$ , a typical PDE model can be written as

\frac{\partial u}{\partial t} = - \nabla \cdot J + S,

(1.1)

where $J$ is the flux of cells and $S$ is a source term that can be used to model proliferation and/or cell death. Different PDE models are specified by choosing different forms of $J$ and $S$ . Within the context of modelling homogenous cell populations, the most common choice for the flux term is based on the assumption that cells move randomly [25], giving rise to linear diffusion with a flux term given by Fick’s law, $J = - D \nabla u$ , where $D > 0$ is the cell diffusivity [3,7,8]. A standard choice for the source term is to specify a logistic term to represent carrying capacity-limited proliferation, $S = λ u (1 - u / K)$ , where $λ > 0$ is the proliferation rate and $K > 0$ is the carrying capacity density [3,7,8]. These choices of $J$ and $S$ mean that equation (1.1) is a multi-dimensional generalization of the well-known Fisher–Kolmogorov model [26–29], which has been successfully used to interpret a number of applications including in vivo tumour progression [18], in vivo embryonic development [10], in vitro wound healing [7,8] and tissue growth [15,16].

Figure 1a shows experimental images of a simple two-dimensional in vitro cell migration experiment, called a barrier assay [30,32]. These experiments are initiated by uniformly placing approximately 30 000 fibroblast cells as a monolayer inside a circular barrier of radius 3 mm. In these experiments, cells are pre-treated with an anti-mitotic drug that prevents proliferation [33], and there is no observed cell death [30]. Accordingly, we model this experiment by setting $S = 0$ in equation (1.1). The experiment proceeds by lifting the barrier at $t = 0$ and observing how the population of cells spreads over time, with the rightmost image in figure 1a showing the extent to which the population has spread after $t = 3$ days. Two key features of this experiment are immediately clear from these images: (i) the population of cells spreads symmetrically with time and (ii) the experimental image at $t = 3$ days shows a clear well-defined sharp front at the leading edge of the population as it spreads. Images in figure 1b show a numerical solution of equation (1.1) with $S = 0$ and the standard choice of linear diffusion, $J = - D \nabla u$ , for a typical choice of $D$ [30]. Consistent with the experiments in figure 1a , we see that the simulated population spreads symmetrically, but plotting the density along the line $y = 0$ in the rightmost panel of figure 1b shows that we have $u > 0$ for all $x$ which is inconsistent with the well-defined sharp fronts at the leading edge in the experimental images. This property of having $u > 0$ for all $x$ persists for all $t > 0$ , which is a well-known deficiency of linear diffusion [34]. Figure 1c shows a numerical solution of equation (1.1) with $S = 0$ and a nonlinear degenerate diffusive flux, $J = - D u \nabla u$ , for a typical choice of $D$ in this model [7,8,16,35–37]. Consistent with the experiments we see that the simulated population spreads symmetrically, and plotting the solution along the line $y = 0$ in the rightmost panel of figure 1c shows that we have a well-defined sharp front: $u > 0$ for $| x | < X (t)$ and $u = 0$ for $| x | \geq X (t)$ , where $X (t)$ is the front location at time $t$ . Full details of our numerical method for solving equation (1.1) are given in appendix A.

Figure 1. ( a ) Experimental images showing a population of non-proliferative fibroblast cells spreading in a two-dimensional barrier assay. The image at $t = 0$ shows the population just as the barrier is lifted, and the image at $t = 3$ days shows the population of migrating cells spreading symmetrically with a sharp front. Images reproduced from Simpson *et al.* [30] with permission. Numerical solutions of equation (1.1) with ( b ) $J = - D \nabla u$ and ( c ) $J = - D u \nabla u$ . Both numerical solutions have $S = 0$ and $u (x, y, 0) = 1$ inside a disc of radius 3 mm, and $u (x, y, 0) = 0$ elsewhere to match the initial distribution of cells in the experiments shown in ( a ). The numerical domain is a square of side length 10, and equation (1.1) is discretized on a $201 \times 201$ uniform mesh. The numerical solution of equation (1.1) at $t = 3$ days is given in the middle panel of ( b , c ), and the details of the density profile are shown in the rightmost panels where $u (x, 0, t)$ is plotted at $t = 0$ (red) and $t = 3$ days (blue). Details of the leading edge of the profiles are highlighted in the green rectangle near $x = 4$ illustrating that the density profile in ( b ) has $u > 0$ at all locations, whereas the density profile in ( c ) has compact support, which is consistent with the experimental images in ( a ). All values of $x$ and $y$ in ( b , c ) measure location in terms of mm to be consistent with the experimental images in ( a ). The numerical solution in ( b ) corresponds to a typical value of $D = D_{1} = 2100 μ m^{2} h^{- 1}$ for linear diffusion [30], and in ( c ) we set $D = D_{2} = 4200 μ m^{2} h^{- 1}$ to satisfy $\int_{0}^{1} D_{1} d u = \int_{0}^{1} D_{2} u d u$ , which ensures that both the linear and nonlinear diffusion models lead to a similar amount of spreading over the experimental time scale [31].

The qualitative comparison between the solution of the linear diffusion equation, the nonlinear degenerate diffusion equation and the experimental images in figure 1 has been made with $S = 0$ so that the continuum PDE model is consistent with the experiments where proliferation is suppressed. However, the difference between spreading cell fronts having sharp or smooth fronts is also relevant for models with $S \neq 0$ [3,7,8]. Throughout the first part of this work, we set $S = 0$ , noting that the difference between smooth and sharp-fronted solutions of equation (1.1) is, in general, determined by the choice of $J$ rather than $S$ . We will come back to this point in §2.5 and provide evidence to support this claim.

As mentioned above, many continuum models of homogeneous cell populations adopt a simple linear diffusive flux, $J = - D \nabla u$ , and this approximation is often made with the implicit or explicit acknowledgement that solutions of this PDE model fail to predict a well-defined front as observed experimentally. In contrast, working with the degenerate nonlinear diffusion model by setting $J = - D u \nabla u$ can lead to a better match with experimental data with well-defined sharp fronts [7,8,16,36,38,39]. With this choice of flux and $S = 0$ , equation (1.1) is also known as the porous medium equation [40–44]. Working with the degenerate nonlinear diffusion model is complicated by the fact that this model is one member of a family of models obtained by setting $J = - D u^{n} \nabla u$ , where $n > 0$ is some constant. Solutions of equation (1.1) with this more general choice of nonlinear flux also lead to symmetric spreading with a well-defined sharp front like we saw in figure 1c for all values of $n > 0$ . These sharp-fronted solutions with compact support are similar to moving boundary problems in the sense that there is a well-defined front location with zero density, and the position of this front evolves with time which we can interpret as a model of the position of the cell front in an experiment [40–44]. The question of how to choose the value of the exponent $n$ remains unclear. For example, Sherratt & Murray [2] studied an in vivo wound healing experiment with $n = 0, 1$ and 4, and showed that all three choices of exponent could be used to make their reaction–diffusion PDE model match their experimental data. Later, Jin et al. [36] studied a series of in vitro scratch assays by setting $n = 0, 0.5, 1, 2, 3$ and $4$ and concluded that $n = 1$ led to the best match to their experimental data without attempting to provide a biological motivation or interpretation of this choice of $n$ . Similarly, McCue et al. [45] studied a series of two-dimensional in vitro wound closure experiments and also found that $n = 1$ provided the best match to their experimental data. Other continuum modelling studies have simply worked with $n = 1$ without explicitly considering other choices of the exponent [14,16,46]. In summary, a key challenge in using continuum PDE models with this generalized nonlinear degenerate diffusivity is that the exponent $n$ often acts as a fitting parameter [38], and lacks a clear biological interpretation. In addition to using these kinds of degenerate diffusion models to interpret biological observations, there is also a great deal of inherent mathematical interest in these models and their solutions [35,47].

An alternative to working with a continuum model to understand the collective spatial spreading, growth and invasion of cell populations is to work with a discrete modelling framework that considers the stochastic motion of individual cells [19,25]. Many kinds of discrete models of cell populations have been implemented to interpret experimental observations ranging from simple lattice-based models [48,49] to more complicated lattice-free [50,51] and vertex-based models [52–55]. An attractive feature of working with discrete models is that experimental images and time-lapse movies showing individual cellular-level behaviours can be translated into a set of individual rules that can be implemented with a stochastic framework to provide a high-fidelity simulation-based model capturing the key biological processes of interest [56,57]. Discrete models can be implemented to visualize snapshots of the spreading population in a way that is directly analogous to performing and imaging an experiment to reveal the positions of individual cells within the population. Another advantage of working with discrete stochastic models is that the discrete mechanism can be coarse-grained into an approximate continuum model, which means that we can encode different individual-level rules into a simulation-based model, and then convert these rules into approximate continuum PDE models, and the solution of these coarse-grained models can be compared with averaged discrete data obtained by repeated simulation [48,58–66]. As described previously, there has been a great deal of effort devoted to understanding how different forms of continuum PDE models predict smooth or sharp-fronted solution profiles; however, far less attention has been devoted to understanding what individual-level mechanisms lead to smooth or sharp fronts in discrete models of cell migration.

All mathematical models discussed so far are simple in the sense that they involve a single PDE or a single population of agents in a discrete framework that can be used to describe the spreading of a homogeneous population of cells. Of course, there are many other more complicated models of collective cell spreading that can lead to sharp-fronted solution profiles. These models include coupled reaction–diffusion models of multiple interacting cell populations [67] as well as discrete models describing multiple populations [68]. Other families of mathematical models include models that describe cell migration that involves biased movement along chemical gradients, such as chemotaxis or haptotaxis [69,70]. Here, we will focus on more fundamental mathematical models of simple homogeneous populations composed of one cell type only, and we do not explicitly consider any biased migration mechanism, such as chemotaxis or haptotaxis.

In this work, we propose a simple, biologically motivated, lattice-based discrete model of collective cell migration and proliferation. The discrete model explicitly models how individual cells in a two-dimensional in vitro experiment produce a biological substrate (e.g. biological macromolecules and extracellular material) that is deposited onto the surface of the tissue culture plate [14,71,72]. Substrate is produced at a particular rate, and deposited locally by individuals within the simulated population. Individual agents within the stochastic model undergo an unbiased random walk at a rate that is proportional to local substrate concentration, and crowding effects are incorporated by ensuring that each lattice site can be occupied by no more than a single agent. As we will demonstrate, this simple biologically inspired mechanism allows us to simulate cell spreading experiments similar to those in figure 1a . Through simulation, we first show that altering the rate of substrate deposition visually impacts the sharpness of the agent density front. A deeper mathematical understanding of these observations is obtained by coarse-graining the discrete mechanism to give a novel PDE model whose solution describes the average behaviour of the stochastic model. One way to interpret this new PDE model is that it naturally describes a linear diffusion mechanism at spatial locations well behind the leading edge of the population, as well as more complicated transport mechanisms at the leading edge of the spreading population that give rise to sharp-fronted solution profiles consistent with experimental observations. We show that averaged data from the discrete model can be very well approximated by numerical solutions of the new continuum limit PDE. In particular, both the continuum and discrete models predict the formation of sharp-fronted density profiles. A careful examination of the new continuum limit PDE model allows us to interpret how the different terms in the model lead to the formation of sharp, sometimes non-monotonic fronts. We conclude this study by incorporating a minimal model of cell proliferation into the discrete model, coarse-graining the proliferative discrete mechanism, and comparing averaged data from the discrete model with proliferation to numerical solutions of the new PDE model.

2. Results and discussion

2.1. Stochastic model and simulations

To account for crowding effects, we implement a lattice-based exclusion process where each lattice site can be either vacant or occupied by, at most, a single agent [48,58]. From this point forward, we will use the word agent to refer to individuals within the simulated population and the word cell to refer to individuals within an experimental population of biological cells. For simplicity, we implement the model on a two-dimensional square lattice with lattice spacing $Δ$ . Each site is indexed by $(i, j)$ , where $i, j \in Z_{+}$ , and each site has position $(x, y) = (i Δ, j Δ)$ . The lattice spacing is taken to be the size of a typical cell diameter [30,48]. In any single realization of the stochastic model, the occupancy of each site $(i, j)$ is a binary variable $U_{i, j}$ , with $U_{i, j} = 1$ if the site is occupied, and $U_{i, j} = 0$ if the site is vacant. Each site is also associated with a substrate concentration, which is a continuous function of time, ${\bar{S}}_{i, j} (t) \in [0, {\bar{S}}_{max}]$ , where ${\bar{S}}_{max}$ is the maximum amount of substrate that can be accommodated at each lattice site. For simplicity, we write $S_{i, j} (t) \in [0, 1]$ , where $S_{i, j} (t) = {\bar{S}}_{i, j} (t) / {\bar{S}}_{max}$ is the non-dimensional substrate density.

A random sequential update method is used to advance the stochastic simulations through time. If there are $N$ agents on the lattice, during the next time step of duration $τ$ , $N$ agents are selected independently, at random, one at a time with replacement, and given the opportunity to move. If the chosen agent is at site $(i, j)$ , the agent will attempt to move with probability $P S_{i, j}$ , where $P \in [0, 1]$ is the probability that an isolated agent will attempt to move during a time interval of duration $τ$ . The target site for all potential motility events is selected at random from one of the four nearest neighbour lattice sites, and the potential motility event will be successful if the target site is vacant [30,48]. Once $N$ potential movement events have been attempted, the density of substrate is updated by assuming that agents deposit substrate at a rate of $Γ$ per time step, so that the amount of substrate at each occupied lattice site increased by an amount $Γ τ$ , taking care to ensure that the maximum non-dimensional substrate density at each site is 1. In addition to specifying initial conditions for the distribution of agents and the initial density of substrate, we must specify values of two parameters to implement the stochastic simulation algorithm: $P \in [0, 1]$ determines the motility of agents, and $Γ > 0$ determines the rate of substrate deposition. With this framework, a typical cell diffusivity is given by $D = P Δ^{2} / (4 τ)$ [48].

To illustrate how the discrete model can be used to model the barrier assay in figure 1a , we perform a suite of simulations summarized in figure 2. The radius of the barrier assay is 3 mm, and a typical cell diameter is approximately 20 µm [30,48]. This means we can simulate the initial placement of cells within the barrier by taking a circular region of radius $3000 / 20 = 150$ lattice sites to represent the disc enclosed by the barrier. The experiments in figure 1a are initiated by placing approximately 30 000 cells uniformly, as a monolayer, within the circular barrier. In the simulations, we have $⌊ π 150^{2} ⌉ = 70 686$ lattice sites within the simulated barrier, and we initialize the simulations by randomly populating each lattice site within the barrier with probability $30 000 / 70 686 \approx 0.42$ . With the discrete model, we can simulate a population of cells with a typical cell diffusivity of $D = 2100 μ$ m² h⁻¹ [30] by choosing $P = 1$ and $τ = 24 / 500$ h. This means that simulating 500 time steps of duration $τ = 24 / 500$ h is equivalent to 1 day in the experiment.

Figure 2. Discrete simulations illustrating the role of the substrate deposition rate $Γ$ . All simulations are performed on a $500 \times 500$ square lattice where the lattice spacing corresponds to 20 μm making the diameter of the simulated population distribution in the left column equal to the diameter of the barrier assay at $t = 0$ in figure 1a . Simulations are initiated by randomly occupying sites within a circular region of radius 150 lattice sites so that the expected number of agents at the beginning of the simulation is 30 000. Simulations are performed by setting $P = 1$ and $τ = 24 / 500$ h, with values of $Γ$ as indicated. Results in ( a ) correspond to initializing $S_{i, j} (0) = 1$ at all lattice sites, whereas results in ( b – d ) correspond to initializing $S_{i, j} (0) = 0$ at all lattice sites. Each day of simulation corresponds to 500 time steps in the discrete model, and snapshots are reported in terms of the $(i, j)$ index of the lattice, which can be re-scaled to give the dimensional coordinates noting that $(x, y) = (i Δ, j Δ)$ . Each panel shows a dashed line indicating the initial placement of the circular barrier. Comparing the extent of the spreading populations with this dashed line gives a visual indication of the extent of spreading.

Results in figure 2a show a preliminary simulation with this initial condition where we set $S_{i, j} (0) = 1$ at all lattice sites at the beginning of the simulation. This first simulation corresponds to the simplest possible case where all lattice sites have the maximum amount of substrate present at the beginning of the experiment which means that the simulation does not depend upon the rate of deposition, $Γ$ . In figure 2a , we see that the population of agents spreads symmetrically, and after 3 days, we have a symmetric distribution of individuals without any clear front at the leading edge of the population. In fact, by $t = 3$ days, we see that some agents within the simulated population become completely isolated, having spread far away from the bulk of the population as a result of chance alone. This situation is inconsistent with the experimental images in figure 1a , where we see a clear front at the leading edge of the population and a complete absence of individuals that become separated from the bulk population. Open-source Julia code to replicate these stochastic simulations is available on GitHub.

Additional simulation results in figure 2b–d involve setting up the same initial distribution of agents as in figure 2a except that we set $S_{i, j} (0) = 0$ at all lattice sites at the beginning of the simulation. These simulations in figure 2b–d are more biologically realistic than the simulations in figure 2a , because in the real experiment, cells are placed into the barriers at the beginning of the experiment without having any chance to deposit significant amounts of substrate onto the surface of the tissue culture plate before the barrier is lifted. Simulations in figure 2b–d are shown for different substrate deposition rates, $Γ$ . If the substrate is deposited sufficiently fast, as in figure 2b , the distribution of individual agents at $t = 3$ days is visually indistinguishable from the case in figure 2a , where we do not observe a clear front in the spreading population. As $Γ$ is reduced, results in figure 2c,d show that the populations spread symmetrically with time, and now we see an increasingly well-defined sharp front as the population of agents spreads. The snapshot of individuals in figure 2d shows that after $t = 3$ days, there are very few individual agents that are isolated away from the bulk of the population, and this distribution is consistent with the experimental observations in figure 1a .

2.2. Continuum limit partial differential equation model

We now provide a greater mathematical understanding and interpretation of the discrete simulation results in figure 2 by coarse-graining the discrete mechanism to give an approximate continuum limit description in terms of a PDE model in the form of equation (1.1). We begin by considering the average occupancy of site $(i, j)$ , where the average is constructed by considering a suite of $M$ identically prepared simulations to give

⟨ U_{i, j} ⟩ (t) = \frac{1}{M} \sum_{m = 1}^{M} U_{i, j}^{m} (t),

(2.1)

where $U_{i, j}^{m} (t)$ is the binary occupancy of lattice site $(i, j)$ at time $t$ in the $m$ th identically prepared realization. With this definition, we treat $U_{i, j} (t) \in [0, 1]$ as a smooth function of time, and for notational convenience, we will simply refer to this quantity as $U_{i, j}$ .

Under these conditions, we can write down an approximate conservation statement describing the change in average occupancy of site $(i, j)$ during the time interval from time $t$ to time $t + τ$ [48,58]:

\begin{aligned} δ ⟨ U_{i, j} ⟩ = & \frac{P}{4} [\underset{migration onto site (i, j)}{\underset{⏟}{(1 - ⟨ U_{i, j} ⟩) \sum [S_{i, j} ⟨ U_{i, j} ⟩]}} - \underset{migration out of site (i, j)}{\underset{⏟}{S_{i, j} ⟨ U_{i, j} ⟩ (4 - \sum ⟨ U_{i, j} ⟩)}}], \end{aligned}

(2.2)

\begin{aligned} δ S_{i, j} = & {\begin{cases} Γ ⟨ U_{i, j} ⟩ for S_{i, j} < 1, \\ 0 for S_{i, j} = 1, \end{cases} \end{aligned}

(2.3)

where, for notational convenience, we write

\begin{aligned} \sum [S_{i, j} ⟨ U_{i, j} ⟩] & = S_{i + 1, j} ⟨ U_{i + 1, j} ⟩ + S_{i - 1, j} ⟨ U_{i - 1, j} ⟩ + S_{i, j + 1} ⟨ U_{i, j + 1} ⟩ + S_{i, j - 1} ⟨ U_{i, j - 1} ⟩, \end{aligned}

(2.4)

\begin{aligned} \sum ⟨ U_{i, j} ⟩ & = ⟨ U_{i + 1, j} ⟩ + ⟨ U_{i - 1, j} ⟩ + ⟨ U_{i, j + 1} ⟩ + ⟨ U_{i, j - 1} ⟩ . \end{aligned}

(2.5)

The first term on the right of equation (2.2) approximately describes the increase in expected occupancy of site $(i, j)$ owing to motility events that would place agents on that site. Similarly, the second term on the right of equation (2.2) approximately describes the decrease in expected occupancy of site $(i, j)$ owing to motility events associated with agents leaving site $(i, j)$ . We describe these terms as approximate as we have invoked the mean-field assumption that the average occupancies of lattice sites are independent [58]. While this assumption is clearly questionable for any particular realization of a discrete model, when we consider the expected behaviour of an ensemble of identically prepared simulations this approximation turns out to be quite accurate [48,58]. Note that setting $S_{i, j} = 1$ at all lattice sites means that this conservation statement simplifies to previous discrete conservation statements that have neglected the role of the substrate [48].

To proceed to the continuum limit, we identify $U_{i, j}$ and $S_{i, j}$ with smooth functions $u (x, y, t)$ and $s (x, y, t)$ , respectively. Throughout this work, we associate uppercase variables with the stochastic model and lowercase variables with the continuum limit model. We expand all terms in equation (2.2) in a Taylor series about $(x, y) = (i Δ, j Δ)$ , neglecting terms of $O (Δ^{3})$ and smaller. Dividing the resulting expressions by $τ$ , we take limits as $Δ \to 0$ and $τ \to 0$ , with the ratio $Δ^{2} / τ$ held constant [25] to give

\begin{aligned} \frac{\partial u}{\partial t} = & D \nabla \cdot [s \nabla u + u (1 - u) \nabla s], \end{aligned}

(2.6)

\begin{aligned} \frac{\partial s}{\partial t} = & {\begin{cases} γ u for s < 1, \\ 0 for s = 1, \end{cases} \end{aligned}

(2.7)

where

D = lim_{\begin{matrix} Δ \to 0 \\ τ \to 0 \end{matrix}} (\frac{P Δ^{2}}{4 τ}), γ = lim_{\begin{matrix} Δ \to 0 \\ τ \to 0 \end{matrix}} (\frac{Γ}{τ}),

(2.8)

which relates parameters in the discrete model: $Δ, τ, P$ and $Γ$ , to parameters in the continuum model: $D$ and $γ$ . In practice, for given values of ∆, τ, P and Γ, we will use the approximation $D = P Δ^{2} / (4 τ)$ and $γ = Γ / τ$

The evolution equation for $s$ , equation (2.7), arises directly from our discrete model where we assume that each lattice site can occupy a maximum amount of substrate. This leads to a mechanism that is very similar to an approach that has been recently adopted to study a generalization of the well-known Fisher–KPP model where the nonlinear logistic source term is replaced with a linear saturation mechanism [73–75]. Solutions of these saturation-type models of invasion involve moving boundaries that form as a result of the saturation mechanism since this provides a natural moving boundary between regions where $s = 1$ and $s < 1$ . Later in §§2.3 and 2.4, we will show that equations (2.6) and (2.7) can also be interpreted as moving boundary problems in exactly the same way as [73–75].

The form of equations (2.6) and (2.7) provides insight into the population-level mechanisms encoded with the discrete model. To see this, we write the flux encoded within equation (2.6) as

J = \underset{diffusive flux}{\underset{⏟}{- D s \nabla u}} - \underset{advective flux}{\underset{⏟}{D u (1 - u) \nabla s}} .

(2.9)

Written in this way, we can now interpret how these two components of the cell flux impact the population-level outcomes. One way to interpret these terms is to note that the first term on the right of equation (2.9) is proportional to $\nabla u$ which is similar to a diffusive flux, and the second term on the right of equation (2.9) is proportional to $u (1 - u)$ which acts like a nonlinear advective flux. In particular, this non-linear advective flux is similar to fluxes often encountered in mathematical models of traffic flow [76].

We can also interpret how the two components of $J$ in equation (2.9) give rise to different features in the solution of the model depending on the location within a spreading population of individuals, such as the discrete populations shown in figure 2. For example, in regions that have been occupied by agents for a sufficiently long period of time, such as regions near the centre of the spreading populations in figure 2 where $u > 0$ , locally we will eventually have $s = 1$ and $\nabla s = 0$ . This means that equations (2.6) and (2.7) simplify the linear diffusion equation since the nonlinear advective flux vanishes and the diffusive-like flux simplifies to Fick’s law of diffusion. In contrast, for regions that have been recently occupied by agents, such as near the leading edge of a population, we have $s < 1$ and $\nabla s \neq 0$ . Under these conditions, the diffusive flux is similar to a nonlinear diffusion term where the diffusive flux of $u$ is proportional to $s$ , which reflects the fact that agent motility in the discrete model is directly proportional to the local density of substrate. The advective-like component of the flux acts like a nonlinear advection term since the flux is proportional to $u (1 - u)$ [76], meaning that the advective flux vanishes when $u = 0$ and $u = 1$ , and is a maximum when $u = 1 / 2$ . The direction of the nonlinear advective flux is opposite to $\nabla s$ . The nonlinear advective flux explicitly includes crowding effects encoded into the discrete model by enforcing that each lattice site can be occupied by, at most, a single agent.

2.3. Continuum–discrete comparison

We now examine how well the numerical solution of equations (2.6) and (2.7) matches averaged data from the discrete model. The experimental images and stochastic simulations in figures 1 and 2 correspond to a radially symmetric polar coordinate system, which can be described by writing equations (2.6) and (2.7) in terms of a radial coordinate system. Instead, we consider a second set of discrete simulations, shown in figure 3, where we consider a rectangular domain with a width of 300 lattice sites, and height of 20 lattice sites. Simulations are initialized by setting $S_{i, j} (0) = 0$ at all lattice sites, and uniformly occupying all sites within $i \leq 150$ with agents. Reflecting boundary conditions are imposed along all boundaries to ensure that the distribution of agents remains, on average, independent of vertical position [48,59], and simulations are performed for $Γ = 10^{2}, 10^{1}, 10^{0}, 10^{- 1}$ and $10^{- 2}$ per time step, as shown in figure 3. Simulation results are consistent with previous results in figure 2, where we see that simulations performed with sufficiently large substrate deposition rates lead population spreading with a smooth front, without any obvious well-defined front position. Simulations with larger $Γ$ lead to population spreading with a visually noticeable defined front. Our main motivation for performing simulations on a rectangular-shaped lattice is that we can work with equations (2.6) and (2.7) in a one-dimensional Cartesian coordinate system [48].

Discrete simulations with, h, and (a–f) various values of, as indicated. All simulations are performed on a rectangular lattice of width and height — Figure 3. Discrete simulations with $P = 1$ , $τ = 24 / 500$ h, and (a–f) various values of $Γ$ , as indicated. All simulations are performed on a rectangular lattice of width $W = 300$ and height $H = 20$ . Simulations are initialized by setting $S_{i, j} = 0$ at all lattice sites, and all sites with $i \leq 150$ are occupied by agents. Snapshots are shown at $t = 1, 2, 3$ and $4$ days. Each day of simulation corresponds to 500 time steps in the discrete model, and snapshots are reported in terms of the $(i, j)$ index of the lattice, which can be re-scaled to give the dimensional coordinates noting that $(x, y) = (i Δ, j Δ)$ .

Averaged agent density data are extracted from the simulations illustrated in figure 3 by considering $M$ identically prepared realizations of the discrete model, averaging the occupancy of each lattice site across these realizations and then further averaging the occupancy along each column of the lattice to give [48]

⟨ U_{i} ⟩ = \frac{1}{H M} \sum_{m = 1}^{M} \sum_{j = 1}^{H} U_{i, j}^{m},

(2.10)

where $H$ is the height of the lattice. Numerical solutions of equations (2.6) and (2.7) in a one-dimensional Cartesian coordinate system are obtained for parameter values and initial data consistent with the discrete simulations. Details of the numerical method used to solve the continuum PDE model are given in appendix A. Results in figure 4 compare numerical solutions of equations (2.6) and (2.7) with averaged data from the discrete simulations, given by equation (2.10) for various values of $Γ$ , as indicated.

Figure 4. Averaged discrete data (dots) superimposed on numerical solutions of equations (2.6) and (2.7) (solid). Each panel compares averaged discrete data, constructed using equation (2.10) with $H = 20$ and $M = 100$ , with a numerical solution of equations (2.6) and (2.7). Four sets of solutions are shown for (a–d) $Γ = 10^{1}, 10^{0}, 10^{- 1}$ and $10^{- 2}$ per time step, as indicated. Discrete simulations are initialized by occupying all lattice sites with $i \leq 150$ , and with $P = 1$ and $τ = 24 / 500$ h. Within each panel, a comparison is made at $t = 0, 1, 2, 3$ and $4$ days shown in blue, green orange and yellow, respectively, as indicated. Each day of simulation corresponds to 500 time steps in the discrete model, and snapshots are reported in terms of the $(i, j)$ index of the lattice, which can be re-scaled to give the dimensional coordinates noting that $(x, y) = (i Δ, j Δ)$ . The arrows within each panel show the direction of increasing time.

Results in figure 4 indicate that the quality of the continuum–discrete match is very good for all values of $Γ$ considered. For sufficiently large values of the substrate deposition rate in figure 4a,b we see that the density profiles are smooth, with no clear well-defined front location at the low-density leading edge. These results are consistent with the preliminary numerical results in figure 1a,b for the barrier assay geometry. In contrast, for sufficiently small values of the substrate deposition rate, density profiles in figure 4c,d show that we have a well-defined sharp front at the leading edge of the spreading populations. The density profiles in figure 4d indicate that the solution of equations (2.6) and (2.7) for $u (x, t)$ has compact support, and the density profiles are non-monotonic with a small dip in density just behind the leading edge. Interestingly, we see the small dip in density behind the leading edge in the averaged discrete data. This indicates that the continuum limit PDE model provides an accurate approximation of the average densities from the discrete simulations. As far as we are aware this dip in density just behind the leading edge has not been measured experimentally.

2.4. Front structure

Now that we have confirmed that averaged data from the discrete model can be approximated by numerical solutions of equations (2.6) and (2.7), we will briefly describe and summarize the general features of the front structure in a simple one-dimensional Cartesian geometry analogous to the results in figure 4. This discussion of the front structure is relevant for initial conditions of the form $s (x, 0) = 0$ for all $x$ , and $u (x, 0) = 1 - H (X)$ , where H is the usual Heaviside step function so that initially we have $u = 1$ for $x < X$ and $u = 0$ for $x > X$ . In all cases considered we impose zero flux boundaries on $u (x, t)$ at both boundaries of the one-dimensional domain. Figure 5 shows a typical solution of equations (2.6) and (2.7). This schematic solution corresponds to the most interesting case with sufficiently small $γ$ that we see a clear sharp-fronted solution, and both $\partial u / \partial x$ and $\partial s / \partial x$ are discontinuous at some moving location $x = η (t)$ . As discussed in §2.2, the moving boundary at $x = η (t)$ arises because of the saturation mechanism governing the dynamics of $s$ in equation (2.7). This kind of moving boundary problem has been previously studied in the case of a generalized Fisher–KPP model [73–75], except that these previous investigations have not involved any discrete stochastic models, or any kind of coarse-graining to arrive at an approximate PDE model.

Schematic front structure. (a) Profiles for (green) and (yellow) within Region 1 for and Region 2 for — Figure 5. Schematic front structure. ( a ) Profiles for $u (x, t)$ (green) and $s (x, t)$ (yellow) within Region 1 for $x < η (t)$ and Region 2 for $η (t) < ξ (t)$ . ( b ) Using the profiles for $u (x, t)$ and $s (x, t)$ in ( a ), we show the corresponding spatial distribution of $J_{d} (x, t) = - D s \partial u / \partial x$ (red) and $J_{a} (x, t) = - D u (1 - u) \partial s / \partial x$ (blue). ( c ) Using the profiles for $u (x, t)$ and $s (x, t)$ in ( a ), we show the spatial distribution of the total flux $J (x, t) = J_{d} (x, t) + J_{a} (x, t)$ (black). The locations of $x = η (t)$ and $x = ξ (t)$ are shown with vertical dashed lines with $ξ (t) > η (t)$ . In ( b ), we include a horizontal dashed line at $J = 0$ to emphasize the point that the diffusive flux changes sign at $x = η (t)$ .

The schematic showing $u (x, t)$ and $s (x, t)$ in figure 5a motivates us to consider two regions within the solution:

—	Region 1: $x < η (t)$ , where $s (x, t) = 1$ .
—	Region 2: $η (t) < x < ξ (t)$ , where $0 < s (x, t) < 1$ .

Ahead of Region 2 where $x > ξ (t)$ we have $u = s = 0$ , and so we consider $x = ξ (t)$ to be the front of the solution. In Region 1, we have $s (x, t) = 1$ and $\partial s / \partial x = 0$ , which means that the evolution of $u (x, t)$ in Region 1 is governed by the linear diffusion equation and the flux of $u$ simplifies to $J = - D \partial u / \partial x$ . This simplification explains why, for this initial condition, $u (x, t)$ is a monotonically decreasing function of $x$ within Region 1 because solutions of the linear diffusion equation obey a maximum principle [77].

Region 2 is characterized by having $s (x, t) < 1$ with $\partial s / \partial x < 0$ . The interface between Region 1 and Region 2 has $s (η (t), t) = 1$ and $u (η (t), t) = u^{*}$ , for some value $0 < u^{*} < 1$ . Within Region 2, the flux of $u$ is given by $J = - D s \partial u / \partial x - D u (1 - u) \partial s / \partial x$ . The advective component of the flux, $- D u (1 - u) \partial s / \partial x$ , is directed in the positive $x$ direction, which means that the flux of $u$ entering Region 2 across the interface at $x = η (t)$ is partly advected in the positive $x$ direction due to the advective flux term that acts within Region 2 only. This additional advective flux in the positive $x$ direction within Region 2 explains why there can be a local minima in $u$ at $x = η (t)$ . The diffusive component of the flux in Region 2, $- D s \partial u / \partial x$ , can act in either the positive $x$ direction when $\partial u / \partial x < 0$ or in the negative $x$ direction when $\partial u / \partial x > 0$ . The schematic in figure 5a shows $u (x, t)$ and $s (x, t)$ across Regions 1 and 2. Associated schematic in figure 5b shows $J_{d} = - D s \partial u / \partial x$ and $J_{a} = - D u (1 - u) \partial s / \partial x$ , and the schematic in figure 5c shows $J = J_{d} + J_{a}$ for the $u$ and $s$ profiles in figure 5a . These plots of the fluxes show that while the total flux $J > 0$ across both Regions 1 and 2, we see that $J_{a}$ vanishes everywhere except within Region 2, and $J_{d} > 0$ within Region 1, but $J_{d}$ changes sign within Region 2 in this case.

Exploring numerical solutions of equations (2.6) and (2.7) indicates that the width of Region 2, $w (t) = ξ (t) - η (t)$ , decreases with $γ$ . This is both intuitively reasonable and consistent with the observations in figure 2 regarding how the structure of the front appeared to vary with the deposition rate in the discrete model. Numerical solutions of equations (2.6) and (2.7) indicate that as $γ \to \infty$ , we have $s (x, t) \to 1 - H (η (t))$ and $w (t) \to 0^{+}$ . Since the width of Region 2 vanishes for sufficiently large $γ$ , the solution of equations (2.6) and (2.7) can be accurately approximated by the solution of the linear diffusion equation, which is independent of $γ$ . Again, this outcome is consistent with the discrete simulations in figure 2, where we observed that simulations with large deposition rates were visually indistinguishable from simulations where all lattice sites were initialized with the maximum substrate concentration where the continuum limit of the discrete model is the linear diffusion equation [48].

The schematic profiles of $u (x, t)$ and $s (x, t)$ in figure 5 can also be interpreted in terms of the mechanisms acting in the discrete model. When agents within Region 2, close to the front, move in the positive $x$ direction to a lattice site that has never been previously occupied, that agent will experience $S_{i, j} = 0$ at the new site. This means that the agent will be stationary for a period of time until that agent deposits substrate, which means that $S_{i, j}$ increases. While there is empty space behind that agent, for example, at site $(i - 1, j)$ where it was previously located, the agent cannot easily move back until a sufficient amount of time has passed to build up the amount of substrate. Therefore, Region 2 within the discrete model acts as a low-motility zone where agents become momentarily stationary until sufficient substrate is produced to enable the agents to continue to move.

In addition to presenting a physical interpretation of the front structure, from both the continuum and discrete points of view summarized here, we also attempted to examine the structure of the front more formally using an interior layer analysis by identifying $γ^{- 1}$ as a small parameter in the system of governing equations. Unfortunately, this leads to a nonlinear PDE as the $O (1)$ problem that determines the shape of the interior layer. Since we are unable to solve the $O (1)$ problem we did not proceed any further with this approach.

2.5. Proliferation

The experimental image in figure 1a shows a barrier assay describing the spatial spreading of a population of fibroblast cells that are pre-treated to prevent proliferation [30,33]. All subsequent discrete and continuum modelling in this work so far has focused on conservative populations without any death or proliferation mechanisms so that these simulations are consistent with the preliminary experimental observations in figure 1. In the discrete model, this is achieved by simulating a population of $N$ agents, where $N$ is a constant. In the continuum model, this is achieved by working with PDE models like equation (1.1) with $S = 0$ . To conclude this study, we now re-examine all discrete and continuum models by incorporating a minimal proliferation mechanism motivated by the additional experimental results summarized in figure 6. The proliferation mechanism involves agent division only without any death process because there is no evidence of cell death in the experiments that motivate these simulations [30]. The leftmost image in figure 6 shows a barrier assay initialized with approximately 30 000 fibroblast cells just after the barrier is lifted at $t = 0$ . The central image in figure 6 shows the outcome of a barrier assay where the fibroblast cells are pre-treated to suppress proliferation [30,33], and the rightmost image shows the outcome of a barrier assay that is initialized in the same way except that the fibroblast cells are not pre-treated to suppress proliferation. This means that the rightmost image in figure 6 shows the outcome of a barrier assay in which fibroblast cells are free to move and proliferate [30]. The motile and proliferative population expands symmetrically, and the leading edge of the population remains sharp. The main difference between the outcome of the barrier assays for the motile and proliferative population compared with the population where proliferation is suppressed is that cell proliferation leads to more rapid spatial expansion of the population.

Figure 6. Circular barrier assay images comparing the spatial spreading of motile and non-proliferative populations with the spatial spreading of a motile and proliferative population of fibroblast cells. The leftmost image shows a barrier assay at $t = 0$ days just after the circular barrier is lifted. This experiment is initiated by placing approximately 30 000 fibroblast cells uniformly inside a barrier of radius 3 mm. The central image shows the spatial extent of the population at $t = 3$ days where the cells are pre-treated to suppress proliferation. The rightmost image shows the spatial extent of the population at $t = 3$ days where the cells are motile and proliferative. All images reproduced from Simpson *et al.* [30] with permission.

A minimal model of proliferation is now incorporated into the discrete model described previously in §2.1. The key difference is that previously the number of agents $N$ remained fixed during the stochastic simulations, whereas now $N (t)$ is a non-decreasing function of time. Within each time step of the discrete model, after giving $N (t)$ randomly selected agents an opportunity to move, we then select another $N (t)$ agents at random, one at a time with replacement, and given the selected agents an opportunity to proliferate with probability $Q \in [0, 1]$ . We take a simple approach and assume that the proliferation is independent of the local substrate density. If a selected agent is going to attempt to proliferate, the target site for the placement of the daughter agent is randomly selected from one of the four nearest neighbour lattice sites [78]. If the target site is occupied then the proliferation event is aborted owing to crowding effects, whereas if the target site is vacant a new daughter agent is placed on the target site. At the end of every time step, we update $N (t)$ to reflect the change in total population owing to proliferation during that time step [48,78]. A set of preliminary simulations comparing results with $Q = 0$ and $Q > 0$ are given in figure 7. In these simulations, we compare the spatial spreading of 30 000 agents uniformly distributed within a circular region of diameter 3 mm. Results in figure 7a are made in the case where we set $S_{i, j} (0) = 1$ at all lattice sites at the beginning of the experiment, and we repeat the comparison for simulations with $S_{i, j} (0) = 0$ with $Γ = 10^{0}, 10^{- 1}$ and $10^{- 2}$ in figure 7b–d , respectively.

Discrete simulations illustrating the role of the substrate deposition rate in the spatial spreading of a population of motile agents without proliferation with the spatial spreading of a population of motile and proliferative agents — Figure 7. Discrete simulations illustrating the role of the substrate deposition rate $Γ$ in the spatial spreading of a population of motile agents without proliferation with the spatial spreading of a population of motile and proliferative agents. All simulations are performed on a $500 \times 500$ square lattice where the lattice spacing corresponds to 20 μm making the diameter of the simulated populations in the left column equal to the diameter of the populations at $t = 0$ in figures 1 and 6. Simulations are initiated by randomly occupying sites within a circular region of radius 150 lattice sites so that the expected number of agents at the beginning of the simulation is 30 000. Simulations of motile and non-proliferative populations correspond to $P = 1$ , $Q = 0$ and $τ = 24 / 500$ days, with values of $Γ$ as indicated. Simulations of motile and proliferative populations correspond to $P = 1$ , $Q = 1 / 500$ and $τ = 24 / 500$ days, with values of $Γ$ as indicated. Results in ( a ) correspond to initializing $S_{i, j} (0) = 1$ at all lattice sites, whereas results in ( b – d ) correspond to initializing $S_{i, j} (0) = 0$ at all lattice sites. Each day of simulation corresponds to 500 time steps in the discrete model, and snapshots are reported in terms of the $(i, j)$ index of the lattice, which can be re-scaled to give the dimensional coordinates noting that $(x, y) = (i Δ, j Δ)$ . Each panel shows a dashed line indicating the initial placement of the circular barrier. Comparing the extent of the spreading populations with this dashed line gives a visual indication of the extent of spreading.

Similar to the experiments in figure 6, our simulations in figure 7 show that incorporating proliferation increases the rate at which the growing populations spread and invade the surrounding area. As for the non-proliferative simulations in figure 2, we see that the front of the spreading populations is poorly defined when $Γ$ is sufficiently large, with a relatively diffuse distribution of agents that includes many isolated individuals that have migrated well ahead of the bulk population. In contrast, reducing $Γ$ leads to visually well-defined sharp fronts with a clearer boundary at the leading edge of the proliferative population. These sharper fronts contain very few isolated individual agents. Visual comparison of the proliferative and non-proliferative snapshots in figure 7c,d indicates that incorporating proliferation leads to an increasingly sharp and well-defined sharp front. These simulations indicate that having substrate-dependent motility and substrate-independent proliferation are sufficient to produce sharp and well-defined fronts in the discrete simulations.

To interpret the differences between the motile populations and the motile and proliferative populations in figure 7, we coarse-grain the discrete model by following a similar approach taken in §2.2. To proceed we write down an approximate conservation statement describing the change in average occupancy of site $(i, j)$ during the time interval from time $t$ to time $t + τ$ :

δ ⟨ U_{i, j} ⟩ = \frac{P}{4} [\begin{aligned} \underset{migration onto site (i, j)}{\underset{⏟}{(1 - ⟨ U_{i, j} ⟩) \sum [S_{i, j} ⟨ U_{i, j} ⟩]}} - \underset{migration out of site (i, j)}{\underset{⏟}{S_{i, j} ⟨ U_{i, j} ⟩ (4 - \sum ⟨ U_{i, j} ⟩)}} \end{aligned}] \begin{aligned} + & \frac{Q}{4} \underset{proliferation onto site (i, j)}{\underset{⏟}{(1 - ⟨ U_{i, j} ⟩) \sum ⟨ U_{i, j} ⟩,}} \end{aligned}

(2.11)

\begin{aligned} δ S_{i, j} = & {\begin{cases} Γ ⟨ U_{i, j} ⟩ for S_{i, j} < 1, \\ 0 for S_{i, j} = 1. \end{cases} \end{aligned}

(2.12)

The new term on the right of equation (2.11) approximately describes the increase in expected density of site $(i, j)$ owing to proliferation events that would place an agent on that site provided that the target site is vacant [48]. To proceed to the continuum limit, we again identify $U_{i, j}$ and $S_{i, j}$ with smooth functions $u (x, y, t)$ and $s (x, y, t)$ , respectively, and expand all terms in (2.2) in a Taylor series about $(x, y) = (i Δ, j Δ)$ , neglecting terms of $O (Δ^{3})$ and smaller. Dividing the resulting expression by $τ$ , we take limits as $Δ \to 0$ and $τ \to 0$ , with the ratio $Δ^{2} / τ$ held constant [25] to give

\begin{aligned} \frac{\partial u}{\partial t} = & D \nabla \cdot [s \nabla u + u (1 - u) \nabla s] + λ u (1 - u), \end{aligned}

(2.13)

\begin{aligned} \frac{\partial s}{\partial t} = & {\begin{cases} γ u for s < 1 \\ 0 for s = 1, \end{cases} \end{aligned}

(2.14)

where

D = lim_{\begin{matrix} Δ \to 0 \\ τ \to 0 \end{matrix}} (\frac{P Δ^{2}}{4 τ}), λ = lim_{\begin{matrix} Δ \to 0 \\ τ \to 0 \end{matrix}} (\frac{Q}{τ}), γ = lim_{\begin{matrix} Δ \to 0 \\ τ \to 0 \end{matrix}} (\frac{Γ}{τ}),

(2.15)

which provides relationships between parameters in the discrete model: $Δ, τ, P, Q$ and $Γ$ , to parameters in the continuum model: $D$ , $λ$ and $γ$ . The additional term in equation (2.13) is simply a logistic source term with carrying capacity of unity, which reflects the fact that the occupancy of any lattice site is limited to a single agent. The numerical method we use to solve equations (2.13) and (2.14) is given in appendix A.

It is straightforward to choose parameters to mimic known biological observations. An important parameter for applying these models to biological experiments is the ratio $P / Q$ , which compares the relative frequency of motility to proliferation events for isolated agents in regions where $S_{i, j} = 1$ . Key parameters in an experiment are the cell diameter $Δ$ , the cell diffusivity $D$ and the proliferation rate $λ$ , which is related to the cell doubling time, $t_{d}$ , by $λ = \log_{e} 2 / t_{d}$ . Using equation (2.15), we have $(Q / P) = Δ^{2} \log_{e} 2 / (4 D t_{d})$ , noting that this ratio is independent of $τ$ . With typical values of $Δ = 20 μ m$ , $D = 2100 μ m^{2} h^{- 1}$ and $t_{d} = 16$ h, we have $Q / P \approx 1 / 500$ , which means setting $P = 1$ and $Q = 1 / 500$ correspond to biologically relevant parameter values of the discrete model. One way of interpreting this choice of parameters is that the average time between proliferation events for an isolated agent is 500 times longer than the average time between motility events in regions where $S_{i, j} = 1$ .

We now repeat the comparison of averaged discrete data with the solution of equations (2.13) and (2.14) in a one-dimensional Cartesian coordinate system for the same domain, initial conditions and parameter values considered previously in figure 4 except now we consider simulations that include proliferation with $Q = 1 / 500$ . The quality of the continuum–discrete match in figure 8 is very good for all values of $γ$ considered. Comparing the solution profiles in figures 4 and 8 shows that the presence of proliferation over a period of 4 days increases the distance that the population front moves in the positive $x$ direction, just as we demonstrated using the stochastic model in figure 7. In addition to noting that the numerical solution of equations (2.13) and (2.14) provides a reasonable match to averaged discrete data, it is important to note that the presence of proliferation in figure 8 does not alter the trends established previously in figure 4 regarding how $Γ$ affects the sharpness of the front, namely that sufficiently large substrate deposition rates lead to smooth-fronted profiles whereas reduced substrate deposition rates lead to sharp-fronted profiles, with the possibility of having a non-monotonic shape.

Figure 8. Averaged discrete data (dots) superimposed on numerical solutions of equations (2.13) and (2.14) (solid). Each panel compares averaged discrete data, constructed using equation (2.10) with $H = 20$ , $M = 100$ , $P = 1$ and $Q = 1 / 500$ with a numerical solution of equations (2.13) and (2.14). Four sets of solutions are shown for (a–d) $Γ = 10^{1}, 10^{0}, 10^{- 1}$ and $10^{- 2}$ , as indicated. Within each panel, a comparison is made at $t = 0, 1, 2, 3$ and $4$ days shown in blue, green, orange and yellow, respectively, as indicated. Each day of simulation corresponds to 500 time steps in the discrete model, and snapshots are reported in terms of the $(i, j)$ index of the lattice, which can be re-scaled to give the dimensional coordinates noting that $(x, y) = (i Δ, j Δ)$ .

Visually comparing results in figures 2 and 7 provides a simple and easy-to-interpret qualitative comparison of the impact of proliferation. In contrast, comparing results in figures 4 and 8 provides a quantitative comparison of the impact of proliferation since these plots involve identical initial conditions, geometry and visualization of agent densities over the same time scales. The only difference is that figure 4 involves motility only whereas figure 8 involves combined motility and proliferation. Results in figure 8 correspond to biologically relevant parameter estimates where $Q / P = 1 / 500 ≪ 1$ where the mean-field approximation is known to be accurate [48,79,80]. Other parameter choices $Q / P$ are not sufficiently small lead to situations where the solution of the mean-field model does not match averaged data from the discrete model as explored in [48,79,80]. Should the reader wish to explore the implications of parameter choices where $Q / P$ is not sufficiently small they may use the open-source code on GitHub to make additional continuum–discrete comparisons. In this work, we restrict our attention to biologically relevant parameter estimates where the mean-field approximation is reasonably accurate.

3. Conclusion and outlook

In this work, we have revisited the question of using continuum PDE models to study spatial spreading and invasion of populations of cells. While many continuum PDE models involve linear diffusion, solutions of these models do not have compact support, and do not replicate clearly defined fronts that are often observed experimentally. Previously, this issue has been addressed by generalizing the linear diffusion flux, $J = - D \nabla u$ , to a degenerate nonlinear diffusion flux, $J = - D u^{n} \nabla u$ , where $n > 0$ . The motivation for working with degenerate nonlinear diffusion is that the flux vanishes when $u = 0$ and the solution of the PDE model has a well-defined sharp front that can match experimental observations [40–43]. While PDE models with this kind of degenerate nonlinear diffusion flux lead to solutions with well-defined sharp fronts, the biological motivation for these models and a biological interpretation of the exponent $n$ remain unclear. In this work, we have revisited the question of modelling spatial spreading and cellular invasion from the point of view of developing a simple lattice-based discrete model. In the discrete model, we assume that agents produce an external substrate (e.g. biomacromolecules and extracellular matrix) that is deposited locally on the lattice, and the rate of randomly directed agent migration is taken to be proportional to the density of substrate at each lattice site. We explicitly incorporate crowding effects in the discrete model by allowing each lattice site to be occupied by, at most, one single agent. This simple, biologically motivated mechanism allows us to model collective spreading and invasion with well-defined sharp fronts provided that the rate of substrate deposition is sufficiently small. Stochastic simulations that mimic the spatial spreading of cells in a two-dimensional circular barrier assay illustrate that our discrete model is capable of replicating key features of the experiment, namely symmetric spreading of the population with a well-defined sharp front at the leading edge of the population.

Coarse-graining the discrete mechanisms leads to a PDE model with a novel flux term that simplifies to linear diffusion in the bulk of the population, and has features similar to a degenerate nonlinear diffusion flux at the leading edge of the population. Importantly, these features arise within the context of a simple, biologically motivated discrete mechanism that is capable of replicating sharp-fronted density profiles and our approach does not involve specifying a degenerate nonlinear diffusivity function that is difficult to relate to biological mechanisms.

Numerical solutions of the new PDE model provide us with a computationally efficient, accurate approximation of averaged data from the stochastic model. Careful examination of the solutions of the PDE indicates that the structure of the leading edge depends upon the rate of substrate deposition. For sufficiently fast substrate deposition the substrate profile approaches a step function at the leading edge of the spreading population, and the nonlinear PDE model simplifies the linear diffusion equation. In contrast, for sufficiently slow substrate deposition the leading edge of the population behaves like a moving boundary problems where the density profile has compact support, and the shape of the density profile at the leading edge can be non-monotonic. The first set of stochastic simulations and coarse-grained PDE models presented in this work focus on conservative populations where cell proliferation and cell death are absent. To understand how the shape of the front could change when considering a proliferative population we present a second set of simulations and coarse-grained PDE models that incorporate a minimal proliferation mechanism. In this case, the coarse-grained PDE model takes the form of a reaction–diffusion model. We solve the new PDE numerically, using biologically motivated parameter values which show that numerical solutions of the PDE model match averaged data from the discrete simulations very well, and confirm that sharp-fronted density profiles occur in the presence of proliferation. In fact, for the biologically motivated parameter values considered in this work, we find that incorporating proliferation tends to sharpen the density fronts at the leading edge relative to non-proliferative stochastic simulations.

There are many options for extending the work presented in this study. One obvious avenue for exploration is to introduce additional details into the discrete model since our approach in this work is to introduce very simple mechanisms only. An interesting option for further examination would be to generalize the transition probabilities in the following way. In the current model, the transition probability for an agent undergoing a motility event from site $(i, j)$ to site $(i + 1, j)$ is proportional to $S_{i, j} U_{i, j} (1 - U_{i + 1, j})$ , which indicates that the transition probability is a linearly increasing function of local substrate density $S_{i, j}$ . An interesting extension would be to generalize the transition probability to be proportional to $g (S_{i, j}) U_{i, j} (1 - U_{i + 1, j})$ , where $0 \leq g (S) \leq 1$ is a smooth function describing how the motility probability for individual agents depends upon the substrate density. In the context of modelling cell migration, it is natural to assume that $g (S)$ is an increasing function. Taking the continuum limit of the discrete mechanism under these circumstances leads to

\begin{aligned} \frac{\partial u}{\partial t} = & D \nabla \cdot [g (s) \nabla u + \frac{d g (s)}{d s} u (1 - u) \nabla s] + λ u (1 - u), \end{aligned}

(3.1)

\begin{aligned} \frac{\partial s}{\partial t} = & {\begin{cases} γ u for s < 1 \\ 0 for s = 1, \end{cases} \end{aligned}

(3.2)

which is a generalization of setting $g (s) = s$ . Returning to our initial discussions in §1, choosing $g (s) = s^{n}$ , for $n > 0$ means that the diffusive flux term in equations (3.1) and (3.2) is analogous to the flux term in the generalized porous medium equation [40,41,44]. All results presented in this study involve working with the simple choice of $g (s) = s$ ; however, generating and comparing averaged discrete data with numerical solutions of equations (3.1) and (3.2) would be very interesting to explore how different choices of $g (s)$ might impact the quality of the discrete–continuum match and the shape of the front. Another extension would be to couple the probability of proliferation in the discrete model to the substrate density. This would, in effect, introduce a substrate-dependent proliferation rate $λ (s)$ into equation (3.1). Again, the question of generating and comparing averaged discrete density data for this generalization would be interesting and a relatively straightforward extension of the current discrete and continuum modelling frameworks established in this work.

Another extension would be to examine long-time travelling wave solutions of equations (2.13) and (2.14) [3]. In the current work, we have limited our examination of this model to relatively short-time simulations of the discrete model and relatively short-time numerical solutions of the continuum limit PDE which is relevant when using these models to mimic standard experimental protocols. Standard experimental protocols examining collective cell migration and proliferation are typically limited to durations of 24 or 48 h [32,36]. This means that for a typical cell line with a doubling time of 12–24 h, these standard experimental protocols last for approximately 1–4 times the cell doubling time. This means that standard experimental protocols are perfectly suited to examine the effects of proliferation that will be evident over these typical time scales. Our theoretical comparison of averaged discrete data and the solution of the continuum limit PDE in figure 8 is relevant for such typical experimental durations since we compare the evolution of the front position over 4 days for a population with a doubling time of 18 h, which is just over five times the doubling time. Despite the fact that we have considered numerical solutions of equations (2.13) and (2.14) over time scales that are five times the doubling time. It is clear that the numerical solutions in figure 8 have not had sufficient time to approach a constant speed, constant shape travelling wave solution [3]. Therefore, taking a more theoretical point of view, it would be mathematically interesting to examine time-dependent numerical solutions of equations (2.13) and (2.14) over much longer time scales and study the resulting travelling wave behaviour as $t \to \infty$ . This could be achieved by transforming the time-dependent PDE model into the travelling wave coordinate, $z = x - c t$ , where $c$ is the long-time asymptotic speed of the travelling wave solutions. Properties of the solution of the resulting dynamical system could then be studied in phase space to provide information about the relationship between parameters in the continuum PDE model and the travelling wave speed $c$ and the shape of the travelling wave profile [3,72]. We leave both of these potential extensions for future consideration.

Ethics

This work did not require ethical approval from a human subject or animal welfare committee.

Data accessibility

All code and data are freely available on GitHub [81].

Declaration of AI use

We have not used AI-assisted technologies in creating this article.

Authors’ contributions

M.J.S.: conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, project administration, resources, software, supervision, visualization, writing—original draft; K.M.M.: conceptualization, formal analysis, investigation, methodology, writing—review and editing; S.W.M.: conceptualization, supervision, writing—review and editing; P.R.B.: conceptualization, supervision, writing—review and editing.

All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

Conflict of interest declaration

We declare we have no competing interests.

Funding

M.J.S. and P.R.B. are supported by the Australian Research Council (DP230100025).

Acknowledgements

We thank the Faculty of Science at QUT for providing K.M.M. with a mid-year research fellowship to support this project.

Appendix A. Numerical methods

Results in figure 1 involve generating numerical solutions of

\frac{\partial u}{\partial t} = \frac{\partial}{\partial x} [D (u) \frac{\partial u}{\partial x}] + \frac{\partial}{\partial y} [D (u) \frac{\partial u}{\partial y}]

(A 1)

on a square domain centred at the origin with side length $L$ . To solve equation (A 1) we discretize all spatial derivative terms on a uniform square mesh with mesh spacing $h$ so that the mesh point with index $(i, j)$ is associated with location $(- L / 2 + (i - 1) h, - L / 2 + (j - 1) h)$ . Applying a standard central difference approximation to the spatial derivative terms in equation (1.1) at the central nodes leads to

\begin{aligned} \frac{d u_{i, j}}{d t} = \frac{1}{2 h^{2}} [(D (u_{i, j}) + D (u_{i + 1, j})) (u_{i + 1, j} - u_{i, j}) - (D (u_{i, j}) + D (u_{i - 1, j})) (u_{i, j} - u_{i - 1, j}) \\ + (D (u_{i, j}) + D (u_{i, j + 1})) (u_{i, j + 1} - u_{i, j}) - (D (u_{i, j}) + D (u_{i, j - 1})) (u_{i, j} - u_{i, j - 1})] . \end{aligned}

(A 2)

This central difference formula is adjusted along the domain boundaries to enforce no-flux boundaries. When we apply this discretization to simulate linear diffusion we set $D (u) = D$ , and when simulating nonlinear degenerate diffusion we set $D (u) = D u^{n}$ . This system of coupled nonlinear ordinary differential equations is solved using the DifferentialEquation.jl package in Julia, which uses automatic time-stepping routines to minimize truncation error. Results in figure 1 are obtained with $h = 0.05$ , which is sufficiently small to ensure that these numerical results are grid independent.

Results in the main text include numerical solutions of equations (2.6) and (2.7) in a one-dimensional Cartesian geometry,

\begin{aligned} \frac{\partial u}{\partial t} = & D \frac{\partial}{\partial x} [s \frac{\partial u}{\partial x} + u (1 - u) \frac{\partial s}{\partial x}] + λ u (1 - u), \end{aligned}

(A 3)

\begin{aligned} \frac{\partial s}{\partial t} = & {\begin{cases} γ u for s < 1 \\ 0 for s \geq 1 \end{cases} \end{aligned}

(A 4)

on $0 < x < L$ . To solve equations (A 3) and (A 4), we discretize all spatial derivative terms on a uniform mesh with mesh spacing $h$ so that the $i$ th mesh point is associated with position $x_{i} = (i - 1) h$ . Applying a standard central difference approximation to the spatial derivative terms in equation (A 3) gives the following system of coupled nonlinear ordinary differential equations at the $i$ th node:

\begin{aligned} \frac{d u_{i}}{d t} = \frac{D}{2 h^{2}} [(s_{i + 1} + s_{i}) (u_{i + 1} - u_{i}) - (s_{i - 1} + s_{i}) (u_{i} - u_{i - 1}) \\ + (u_{i + 1} [1 - u_{i + 1}] + u_{i} [1 - u_{i}]) (s_{i + 1} - s_{i}) - (u_{i - 1} [1 - u_{i - 1}] + u_{i} [1 - u_{i}]) (s_{i} - s_{i - 1})] \\ + λ u_{i} (1 - u_{i}), \end{aligned}

(A 5)

\begin{aligned} \frac{d s_{i}}{d t} = & {\begin{cases} γ u for s < 1 \\ 0 for s \geq 1. \end{cases} \end{aligned}

(A 6)

The discrete equation for $s$ , equation (A 6), holds for all mesh points $i = 1, 2, \dots, I$ because there are no spatial derivative terms in equation (A 4) and no boundary conditions need to be imposed. In contrast, the discrete equation for $u$ , equation (A 5), holds only on the interior mesh points $i = 2, 3, \dots, I - 1$ . Applying no-flux boundary conditions at $i = 1$ and $i = I$ means that we impose the constraints $u_{1} = u_{2}$ and $u_{I - 1} = u_{I}$ , respectively. This system of coupled nonlinear ordinary differential equations is solved using the DifferentialEquation.jl package in Julia, which implements automatic time-stepping routines to control temporal truncation error. All numerical results in this work correspond to $h = 0.1$ , which is sufficiently small to ensure that our numerical results are grid independent for the problems that we consider. Open source Julia code to solve equations (A 5) and (A 6) is available on GitHub. Results in this study are obtained using a Runge–Kutta method within the Tsit5 algorithm within the DifferentialEquations.jl package.

Footnotes

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.

Comments

VIEW FULL TEXT DOWNLOAD PDF

This Issue

May 2024

Volume 11Issue 5

Article Information

Published by:Royal Society
Online ISSN:2054-5703

History:

Manuscript received22/01/2024
Manuscript accepted22/02/2024
Published online15/05/2024

License:

Citations and impact

Keywords

Subjects

Your password must have 8 characters or more and contain 3 of the following:

Discrete and continuous mathematical models of sharp-fronted collective cell migration and invasion

Abstract

1. Introduction

2. Results and discussion

2.1. Stochastic model and simulations

2.2. Continuum limit partial differential equation model

2.3. Continuum–discrete comparison

2.4. Front structure

2.5. Proliferation

3. Conclusion and outlook

Ethics

Data accessibility

Declaration of AI use

Authors’ contributions

Conflict of interest declaration

Funding

Acknowledgements

Appendix A. Numerical methods

Footnotes

Comments

This Issue