1. Introduction
1.1. Biological background
Phenomena of collective cell migration have received a significant amount of interest in recent years, and a particularly large literature has been devoted to their role during cancer invasion processes [Reference Friedl, Locker, Sahai and Segall27]. Histopathological analysis of tissue specimens reveals a plurality of patterns within invading cell fronts, ranging from individually migrating cells to collective strands and clusters that infiltrate the surrounding healthy tissue. Increased invasiveness forms one of the key traits of cancer metastasis, and consequently presents a significant impediment to successful treatment. Understanding the underlying biological processes is therefore of manifest interest.
Invading cell fronts frequently present significant phenotypic heterogeneity, and particular attention has focussed on the extent to which a separation into ‘leader’ and ‘follower’ cells contributes to invasive spread [Reference Vilchez Mercedes, Bocci, Levine, Onuchic, Jolly and Wong68]. Leader cells are those that (seemingly) drive the invasion process: positioned at or near the front, modifying the extracellular matrix (ECM) – i.e. the network of macromolecules in which cells reside, and coordinating with follower cells to facilitate their collective invasion into healthy tissue. These leader cells can be derived from various sources, both from tumour cells that have undergone a phenotypic transition (e.g. following genetic mutations or epimutations) or from surrounding stromal cells, such as fibroblasts, that have been activated and co-opted through factors in the tumour microenvironment [Reference Vilchez Mercedes, Bocci, Levine, Onuchic, Jolly and Wong68].
For carcinomas – cancers of epithelial origin – tumour-derived leader cells will typically have undergone an epithelial to mesenchymal transition (EMT), i.e. undergone a loss of their epithelial nature and acquired mesenchymal characteristics [Reference Brabletz, Kalluri, Nieto and Weinberg11]. Epithelial cells are often tightly bound through strong cell-to-cell adhesion, therefore downregulation of adhesion can release leaders from the bonds that control their normal position in the tissue. Mesenchymal characteristics involve (significantly) enhanced motility and strong interactions with the ECM. These interactions include a capacity to modify the matrix and microenvironment, both mechanically and chemically, in a manner that can facilitate the invasion of other leader and follower cells. First, increased production of fibronectin can increase matrix adhesivity, allowing cells to gain more traction [Reference Konen, Summerbell and Dwivedi38, Reference Schwager, Taufalele and Reinhart-King60, Reference Vilchez Mercedes, Bocci, Levine, Onuchic, Jolly and Wong68]; the directed movement that results from migration up a matrix adhesivity gradient is referred to as haptotaxis [Reference Carter14]. Cells may also realign fibres, leading to an oriented matrix that directs invasion along certain paths [Reference Ray and Provenzano58]. Leader cells may also start to secrete (or increase the secretion of) matrix-degrading enzymes (MDEs) [Reference Chen, Wu, Tang, Tang and Liang17, Reference Kessenbrock, Plaks and Werb36, Reference Schwager, Taufalele and Reinhart-King60, Reference Vilchez Mercedes, Bocci, Levine, Onuchic, Jolly and Wong68, Reference Winkler, Abisoye-Ogunniyan, Metcalf and Werb73]. In turn, this can reduce the volume fraction occupied by the ECM and liberate space into which cells can migrate (or proliferate). Beyond these physical alterations to the microenvironment, leader cells may also drive invasion through chemical means, e.g. through secreted factors leading to chemoattractant gradients that direct follower cell movement [Reference Chen, Wu, Tang, Tang and Liang17, Reference Helvert, Storm and Friedl32, Reference Konen, Summerbell and Dwivedi38, Reference Vilchez Mercedes, Bocci, Levine, Onuchic, Jolly and Wong68].
While leader cells may have enhanced motility and greater capacity to alter the surrounding microenvironment, trade-offs may arise in the form of reduced proliferative potential. Analyses into heterogeneous groups composed from distinct follower and leader subpopulations indicate that the former can be significantly more proliferative [Reference Konen, Summerbell and Dwivedi38]. Furthermore, followers can produce factors that promote a level of growth within the leader cells, hence maintaining the balance between leaders and followers across the overall cell population [Reference Konen, Summerbell and Dwivedi38]. Thus, collective invasion in certain leader-follower cancer cell populations may be a cooperative process, with each subpopulation playing a key role in promoting tumour expansion.
Dichotomising invading cells into leader and follower subtypes can be notionally convenient, but masks the possibility that cells may fall into (significantly) more than two fundamental phenotypic states. Unlike the ‘regulated’ EMTs that occur during embryogenesis or wound healing, EMT within cancer cells can be highly variable, ranging from ‘partial’ to ‘full’ [Reference Brabletz, Kalluri, Nieto and Weinberg11, Reference Chen, Wu, Tang, Tang and Liang17, Reference Kroger39, Reference Nieto, Huang, Jackson and Thiery48, Reference Revenu and Gilmour59, Reference Vilchez Mercedes, Bocci, Levine, Onuchic, Jolly and Wong68]. Partial here indicates a cell has only acquired a subset of mesenchymal characteristics, for example, resulting in only a slight increase in motility, or a part reduction of proliferative potential. Moreover, it has been shown that cells can undergo changes to these phenotypes over time [Reference Pastushenko52, Reference Vilchez Mercedes, Bocci, Levine, Onuchic, Jolly and Wong68, Reference Williams, Gao, Redfern and Thompson72, Reference Yu74]. Consequently, a more accurate picture of leader–follower heterogeneity in invading cancers would be that each cell occupies a fluid position within a quasi-continuous phenotype space. Investigating the role of trade-offs in leader–follower collective migration may form a beneficial treatment target. In fact, it has been shown that leader cells are generally more resistant to treatment [Reference Bocci, Levine, Onuchic and Jolly9, Reference Vilchez Mercedes, Bocci, Levine, Onuchic, Jolly and Wong68]; however, if leader cells are removed then the invasion of the tumour stops [Reference Konen, Summerbell and Dwivedi38]. Furthermore, through inhibiting MDE secretion by leader cells then further invasion can be slowed or prevented [Reference Carey, Starchenko, McGregor and Reinhart-King13].
1.2. Mathematical modelling background
The use of mathematical modelling to investigate the invasion processes of cancer cells forms a well-developed area of research, and for more details we refer to the extensive reviews detailed in [Reference Alfonso, Talkenberger and Seifert2, Reference Franssen, Lorenzi, Burgess and Chaplain26, Reference Sfakianakis and Chaplain61]. A significant part of the literature on this subject has focussed on haptotaxis-fuelled invasion, in which a cancer cell population secretes proteolytic factors that alter the surrounding matrix to generate adhesivity gradients. Early models in this field have been formulated as coupled systems of partial differential equations (PDEs) and infinite-dimensional ordinary differential equations (ODEs), which rely on the assumption that cell migration results from the superposition of random motion, modelled as linear diffusion, and haptotaxis, modelled as advection according to the ECM gradient, e.g. [Reference Anderson, Chaplain, Newman, Steele and Thompson5, Reference Perumpanani and Byrne56]. Focussing on a 1D spatial scenario, as a prototypical example of these haptotaxis models of cancer invasion, we can consider the following system
\begin{equation} \begin{cases} \displaystyle{\partial _t \rho } - \partial _x \Big (D \, \partial _x \rho - \rho \, \mu \, \partial _x E\Big ) = R(\rho ) \, \rho,\\[5pt] \displaystyle{\partial _t M} - D_M \, \partial ^2_{xx} M = \gamma \, \rho - \kappa _M \, M,\\[5pt] \displaystyle{\partial _t E}=-\kappa _E \, E \, M, \end{cases} \quad x \in \mathbb{R}, \end{equation}
where the functions
$\rho \equiv \rho (t,x)$
,
$M \equiv M(t,x)$
and
$E \equiv E(t,x)$
model, respectively, the cancer cell density, the MDE concentration and the ECM density at time
$t \in \mathbb{R}^+$
and position
$x \in \mathbb{R}$
. In the system (1), the parameters
$D$
and
$\mu$
are the random motility coefficient and the coefficient of haptotaxis sensitivity (i.e. sensitivity to matrix adhesivity gradients) of cancer cells, respectively, while the function
$R(\rho )$
is the net growth rate of the cell population due to the proliferation (i.e. division and death) of cancer cells. The dependence of this function on cell density
$\rho$
takes into account density-dependent inhibition of growth (i.e. the fact that cessation of cell division occurs once a critical value of the cell density is reached). Moreover, the parameter
$D_M$
is the diffusion coefficient of MDEs, the parameter
$\gamma$
is the rate of MDE production by cancer cells and the parameter
$\kappa _M$
is the rate of natural decay of MDEs. Finally, the parameter
$\kappa _E$
is the rate at which the ECM is degraded by MDEs upon contact.
Haptotaxis models of cancer invasion of the form of (1) have been studied analytically through proof of local existence, global existence, boundedness and uniqueness of solutions [Reference Tao and Cui66], investigations into blow-up of solutions [Reference Shangerganesh, Nyamoradi, Sathishkumar and Karthikeyan62] and analysis of 2D radially symmetric solutions [Reference Bortuli, Freire and Maidana10]. Numerical simulations have also been used to investigate models of the form of (1) [Reference Ganesan and Lingeshwaran28]. Further models extend systems of this form, for example, to include more detailed mechanisms of ECM remodelling and enzyme activities [Reference Andasari, Gerisch, Lolas, South and Chaplain3, Reference Chaplain and Lolas15, Reference Nguyen Edalgo and Ford Versypt47], or through the inclusion of non-local terms to incorporate the effects of cell–cell adhesion [Reference Gerisch and Chaplain29, Reference Painter, Armstrong and Sherratt51]. Models have also been formulated to include detailed cell mechanics – e.g. stress, strain, elasticity, adhesion, transport by velocity fields and other interaction forces – in the context of invasive melanoma growth through different skin layers [Reference Ciarletta, Foret and Amar18].
While model (1) has been restricted to a homogeneous cancer population, cognate models have also been developed that explore the consequences of phenotypic heterogeneity on invasion. Binary state models consider a division of cancer cells into two phenotypic states. These include models in which cells of the invading population switch between a proliferating state and a migrating state, to investigate ramifications of the “go-or-grow” hypothesis for glioma growth [Reference Kolbe, Sfakianakis, Stinner, Surulescu and Lenz37, Reference Pham, Chauviere and Hatzikirou57, Reference Stepien, Rutter and Kuang63]. Invasion models that feature two competing phenotypes with distinct migratory and proliferation properties have also been formulated in the context of acid-mediated invasion, where the heterogeneity extends to distinct acid-resistance and matrix-altering behaviour [Reference Strobl, Krause, Damaghi, Gillies, Anderson and Maini65].
Greater phenotypic heterogeneity can be accounted for through the inclusion of more discrete states, but this becomes impractical if the phenotypic space becomes almost continuous. As such, an alternative approach is to extend a model like (1) to include a continuous structuring variable representing intercellular variability in certain phenotypic characteristics [Reference Perthame54]. In the context of cell invasion type dynamics, though not specifically in cancer, models of this nature have been developed to describe how a trade-off between chemotactic ability and proliferation may shape the phenotypic structuring of chemotaxis-driven growth processes [Reference Lorenzi and Painter40], and how trade-offs between mobility and proliferation may impact on density- or pressure-driven growth processes [Reference Lorenzi, Perthame and Ruan41, Reference Macfarlane, Ruan and Lorenzi46]; directly relevant to cancer, structured phenotype models of this type have also been developed to explore avascular tumour growth [Reference Fiandaca, Bernardi, Scianna and Delitala23] and the evolutionary dynamics underpinning the emergence of intra-tumour phenotypic heterogeneity [Reference Fiandaca, Delitala and Lorenzi24, Reference Lorenzi, Venkataraman, Lorz and Chaplain42, Reference Villa, Chaplain and Lorenzi69]. In recent work by Guilberteau et al. [Reference Guilberteau, Jain, Jolly, Duteil and Pouchol31], the authors presented a PIDE model that captures transitions between fully-epithelial, hybrid epithelial/mesenchymal and fully-mesenchymal cell states and demonstrated this model to be capable of reproducing experimental observations into the dynamics of EMT. This model, however, does not account for spatial dynamics or invasion processes of cells.
While the above discussions have focussed on continuum models, which provide a population-level description of cell dynamics, it is important to note that a very large number of modelling studies have explored haptotaxis-driven cancer invasion using individual-based models (i.e. models that track the dynamics of single cells) [Reference Wang, Butner, Kerketta, Cristini and Deisboeck70, Reference West, Robertson-Tessi and Anderson71]. Advantages lie in the ability of these models to capture the dynamics and stochasticity of single-cell movement, and a notable early example within the context of cancer invasion was developed in [Reference Anderson, Chaplain, Newman, Steele and Thompson5]. Here, a model of the type of (1) was discretised in space, with the discretised terms subsequently used to specify probabilities of movement in different directions, according to the ECM distribution. Extensions were introduced in the model considered in [Reference Anderson, Weaver, Cummings and Quaranta6], where each cell was allowed to undergo random movement, haptotaxis up the gradient of the ECM, produce MDEs, consume oxygen and undergo cell cycle controlled proliferation depending on the availability of oxygen and free space. Moreover, this model comprised cells of different discrete phenotypic states, controlling aspects such as each cell’s adhesion, oxygen consumption, haptotactic ability, secretion rate of MDEs and proliferative potential. Recently, this underlying modelling framework has been further extended to investigate the role of two specific phenotypes, namely epithelial and mesenchymal phenotypes, in cancer invasion and metastasis [Reference Franssen, Lorenzi, Burgess and Chaplain26].
The original framework in the above method constituted of starting with a continuum model and subsequently discretising it to derive the governing rules for cell movement in a corresponding individual-based model [Reference Anderson, Chaplain, Newman, Steele and Thompson5]. An alternative approach for transitioning between discrete and continuum descriptions of cell motion is to first postulate a model at the single-cell level and then employ coarse-graining procedures to derive a continuous description; these methods have been extensively adopted in recent decades, particularly in the context of motivating PDE models to describe taxis-like behaviours, e.g. [Reference Painter50, Reference Penington, Hughes and Landman53, Reference Stevens and Othmer64].
1.3. Outline
In this paper, we consider the following generalisation of the haptotaxis model of cancer invasion (1), where the continuous structuring variable
$y \in [0,Y] \subset \mathbb{R}^+$
, with
$Y \gt 0$
, represents the cell phenotype (i.e. the position of the cells in phenotypic space) and captures intercellular variability in haptotactic response, proliferative potential and production of MDEs:
\begin{equation} \begin{cases} \displaystyle{\partial _t n - \partial _x\Big (D \, \partial _x n - n \, \chi (y) \, \partial _x E\Big ) = R(y,\rho ) \, n + \lambda \, \partial _{yy} n},\\[5pt] \displaystyle{\rho (t,x) \;:\!=\; \int _0^Y n(t,x,y)\ \mathrm{d}y},\\[5pt] \displaystyle{\partial _t M - D_M\partial _{xx}M = \int _0^Y p(y) \, n(t,x,y) \ \mathrm{d}y- \kappa _M M},\\[5pt] \displaystyle{\partial _t E=-\kappa _E \, E \, M}, \end{cases} \quad (x,y) \in \mathbb{R} \times (0,Y). \end{equation}
Compared to model (1), here the PDE for the cell density,
$\rho (t,x)$
, is replaced by the partial integro-differential equation (PIDE) (2)
$_1$
for the local cell population density function,
$n \equiv n(t,x,y)$
, which is linked to the cell density through the relation (2)
$_2$
. Moreover, the functions
$\chi (y)$
,
$p(y)$
and
$R(y,\rho )$
are, respectively, the haptotaxis sensitivity coefficient, the MDE production rate and the net growth rate of the cell population density due to the proliferation (i.e. division and death) of cancer cells with phenotype
$y$
. Finally, the diffusion term on the right-hand side of the PIDE (2)
$_1$
models the effect of phenotypic changes, which occur at rate
$\lambda$
.
We first formulate a phenotype-structured individual-based haptotaxis model of cancer invasion (cf. Section 2), where the dynamics of individual cells are governed by a set of rules that result in a branching random walk over a lattice [Reference Hughes34], which represents both physical and phenotype spaces. In this model, the rules governing cell dynamics are coupled with a balance equation for the MDE concentration and a balance equation for the density of ECM. Then, using an extension of the limiting procedure that we previously employed in [Reference Bubba, Lorenzi and Macfarlane12, Reference Chaplain, Lorenzi and Macfarlane16, Reference Macfarlane, Chaplain and Lorenzi45, Reference Macfarlane, Ruan and Lorenzi46], we formally derive the model (2) as the continuum limit of this individual-based model (cf. Section 3 and Appendix A.1). After that, building upon the formal asymptotic method that we developed in [Reference Lorenzi and Painter40, Reference Lorenzi, Perthame and Ruan41], we carry out travelling wave analysis of an appropriately rescaled version of the model (2) (cf. Section 4). The results obtained provide a mathematical formalisation for the idea that trade-offs between proliferative potential and the ability to sense spatial gradients of ECM and produce MDEs may promote the emergence of phenotypically structured invading cell fronts. Specifically, wherein leader cells (i.e. cells with a higher haptotactic and MDE production ability but a lower proliferative potential) are localised at the leading edge of the front, while follower cells (i.e. cells with a higher proliferative potential but a lower haptotactic and MDE production ability) occupy the region behind the leading edge. Finally, we report on numerical solutions of such a rescaled continuum model and numerical simulations of the corresponding rescaled version of the individual-based model, and we compare them with the results of travelling wave analysis (cf. Section 5). We conclude with a discussion of our findings and propose some future research directions (cf. Section 6).
2. The individual-based model
In this section, integrating the modelling approaches that we developed in [Reference Bubba, Lorenzi and Macfarlane12, Reference Macfarlane, Ruan and Lorenzi46], we formulate a phenotype-structured individual-based haptotaxis model of cancer invasion. In this model, individual cells are represented as agents, while the density of ECM and the concentration of MDEs are described by non-negative functions. We allow cells to undergo undirected random movement, phenotype-dependent haptotactic movement in response to the ECM, heritable spontaneous phenotypic changes (i.e. heritable phenotypic changes that occur randomly and are not biased by the cell microenvironment) and phenotype-density-dependent proliferation (i.e. division and death). We consider the scenario where cells also perform phenotype-dependent secretion of MDEs, which then diffuse throughout the spatial domain according to Fick’s first Law of diffusion and undergo natural decay. Furthermore, the MDEs break down the ECM to create a gradient that affects the haptotaxis of cancer cells.
Focussing on a 1D spatial scenario, we let the cells, the density of ECM and the concentration of MDEs be distributed along the real line
$\mathbb{R}$
. Furthermore, we describe the phenotypic state of each cell by means of a structuring variable
$y\in [0,Y]\subset \mathbb{R}^+$
, which takes into account the intercellular variability in haptotactic sensitivity, MDE secretion rate and proliferation rate.
In particular, we consider the case where larger values of the structuring variable
$y$
correspond to a higher ability to sense spatial gradients of ECM and produce MDEs but a lower proliferative potential (cf. Figure 1). This choice is motivated by the energetic costs associated with enhanced motility and greater capacity to alter the surrounding microenvironment, which lead to trade-offs in the form of reduced proliferative potential [Reference Konen, Summerbell and Dwivedi38].
Figure 1. Schematic overview of the trade-offs between haptotactic and MDE production ability and proliferative potential incorporated into the individual-based model.
Hence, cells in phenotypic states characterised by values of
$y$
closer to
$0$
display a more epithelial-like phenotype (i.e. they behave more like follower cells), whereas cells in phenotypic states characterised by values of
$y$
closer to
$Y$
display a more mesenchymal-like phenotype (i.e. they behave more like leader cells).
We discretise the time variable
$t\in \mathbb{R}^+$
and the space variable
$x\in \mathbb{R}$
, respectively, as
$t_k=k\tau$
and
$x_i=i\Delta _x$
with
$k\in \mathbb{N}_0$
,
$\tau \in \mathbb{R}^+_*$
,
$i\in \mathbb{Z}$
and
$\Delta _x \in \mathbb{R}^+_*$
, where
$\mathbb{R}^+_*$
denotes the set of positive real numbers. Moreover, we discretise the phenotype variable via
$y_j=j\Delta _y\in [0,Y]$
with
$j\in \mathbb{N}_0$
and
$\Delta _y \in \mathbb{R}^+_*$
. Here,
$\tau$
,
$\Delta _x$
and
$\Delta _y$
are the time-step, space-step and phenotype-step, respectively.
Each individual cell is represented as an agent that occupies a position on the lattice
$\{x_i\}_{i\in \mathbb{Z}}\times \{y_j\}_{j\in \mathbb{N}_0}$
and we introduce the dependent variable
$N_{i,j}^k\in \mathbb{N}_0$
to model the number of cells in the phenotypic state
$y_j$
at position
$x_i$
at time
$t_k$
. The cell population density and the corresponding cell density are then defined, respectively, as
\begin{equation} n_{i,j}^k\equiv n(t_k,x_i,y_j)\;:\!=\;\frac{N_{i,j}^k}{\Delta _x \Delta _y} \quad \text{and} \quad \rho _i^k\equiv \rho (t_k,x_i)\;:\!=\;\Delta _y \sum _j n_{i,j}^k. \end{equation}
The concentration of MDEs and the density of ECM at position
$x_i$
at time
$t_k$
are denoted by
$M_i^k \equiv M(t_k,x_i)$
and
$E_i^k \equiv E(t_k,x_i)$
, respectively.
The biological mechanisms incorporated into the model and the corresponding modelling strategies are summarised by the schematics in Figure 2 and are described in the remainder of this section.
Figure 2. Schematic overview of the mechanisms incorporated into the individual-based model along with the corresponding modelling strategies. Between time-steps
$k$
and
$k+1$
, each cell in phenotypic state
$y_j\in (0,Y)$
at spatial position
$x_i\in \mathbb{R}$
may: a. move via random motion to either of the positions
$x_{i-1}$
and
$x_{i+1}$
with probabilities
${P}_{L_{i,j}}^k$
and
${P}_{R_{i,j}}^k$
defined via (4) or do not undergo random movement with probability
${P}_{S_{i,j}}^k=1-({P}_{L_{i,j}}^k+{P}_{R_{i,j}}^k)$
; b. move via haptotactic motion to either of the positions
$x_{i-1}$
and
$x_{i+1}$
with probabilities
${P}_{HL_{i,j}}^k$
and
${P}_{HR_{i,j}}^k$
defined via (6) or do not undergo haptotactic movement with probability
${P}_{HS_{i,j}}^k=1-({P}_{HL_{i,j}}^k+{P}_{HR_{i,j}}^k)$
; c. undergo a phenotypic change and thus enter either of the phenotype states
$y_{j-1}$
and
$y_{j+1}$
with probabilities
${P}_{D_{i,j}}^k$
and
${P}_{U_{i,j}}^k$
defined via (7) or remain in the same phenotypic state with probability
${P}_{N_{i,j}}^k=1-({P}_{D_{i,j}}^k+{P}_{U_{i,j}}^k)$
; d. die and divide with probabilities
${P}_{A_{i,j}}^k$
and
${P}_{B_{i,j}}^k$
defined via (8) or remain quiescent with probability
${P}_{Q_{i,j}}^k=1-({P}_{A_{i,j}}^k+{P}_{B_{i,j}}^k)$
. The concentration of MDEs will change over time through: e. diffusion at the rate
$D_M$
; f. secretion by cells at the phenotypic-dependent rate
$p$
; g. natural decay at rate
$\kappa _M$
. The ECM density will change over time through: h. degradation by MDEs at rate
$\kappa _E$
.
2.1. Modelling the dynamics of cancer cells
As summarised in Figure 2a-d, between time-steps
$k$
and
$k+1$
, each cell in phenotypic state
$y_j\in (0,Y)$
at position
$x_i \in \mathbb{R}$
can undergo undirected random movement and haptotactic movement (which are regarded as independent processes), heritable spontaneous phenotypic changes and cell division and death according to the rules provided in the following subsections.
2.1.1. Random cell movement
We model undirected cell movement as a random walk along the spatial dimension, with movement probability
$0\lt \theta \le 1$
. In particular, for a focal cell in the phenotypic state
$y_j$
at spatial position
$x_i$
at time
$t_k$
, we define the probability of moving left or right to spatial positions
$x_{i-1}$
or
$x_{i+1}$
as
${P}_{L_{i,j}}^k$
or
${P}_{R_{i,j}}^k$
, respectively. As we consider this random movement to be undirected and not affected by the cell phenotype, we define
\begin{equation}{P}_{L_{i,j}}^k={P}_{R_{i,j}}^k \;:\!=\;\frac{\theta }{2}. \end{equation}
Note that cells will not undergo random movement with probability
\begin{equation*} P_{S_{i,j}}^k \;:\!=\; 1-\left ({P}_{L_{i,j}}^k+{P}_{R_{i,j}}^k\right ). \end{equation*}
2.1.2. Haptotactic cell movement
We model haptotactic cell movement in response to the ECM as a biased random walk along the spatial dimension. Since cells move up the gradient of the ECM (i.e. they move towards higher ECM densities), we let the haptotactic movement probabilities depend on the difference between the ECM density at the position occupied by the cell and the ECM density at neighbouring positions. Furthermore, we consider the case where larger values of
$y_j$
correlate with a higher haptotaxis sensitivity (cf. Figure 1). Hence, we modulate the probabilities of haptotactic cell movement by the function
$\mu (y_j)$
, which provides a measure of the sensitivity to matrix adhesivity gradients of cells in phenotypic state
$y_j$
and thus satisfies the following assumptions
\begin{equation} \mu (0)=0, \quad \frac{\mathrm{d} \mu (y)}{\mathrm{d} y}\gt 0 \;\; \text{ for} \quad y\in (0,Y]. \end{equation}
We then assume that between time-steps
$k$
and
$k+1$
a cell in phenotypic state
$y_j$
at position
$x_i$
may move to the position
$x_{i-1}$
(i.e. move left) with probability
${P}_{HL_{i,j}}^k$
or move to the position
$x_{i+1}$
(i.e. move right) with probability
${P}_{HR_{i,j}}^k$
, where we define
\begin{equation}{P}_{HL_{i,j}}^k \;:\!=\; \eta \, \mu (y_j) \frac{\left (E^k_{i-1}-E^k_{i}\right )_+}{2 \, E_{\textrm{max}}}, \quad{P}_{HR_{i,j}}^k \;:\!=\; \eta \, \mu (y_j) \frac{\left (E^k_{i+1}-E^k_{i}\right )_+}{2 \, E_{\textrm{max}}}, \quad \text{with} \quad (\!\cdot\!)_+\;:\!=\;\max\!(0,\cdot ). \end{equation}
Here,
$E_{\textrm{max}} \in \mathbb{R}^+_*$
is the maximum value of the ECM density before cell invasion starts (see also Section 2.2.2). Moreover, the parameter
$\eta \in \mathbb{R}^+_*$
is a scaling factor which we consider small enough to ensure
$\eta \, \mu (y_j) \leq 1$
. Hence, the quantities defined via (6) satisfy
$0\lt{P}_{HL_{i,j}}^k+{P}_{HR_{i,j}}^k \le 1$
for all values of
$i$
,
$j$
and
$k$
. Note that cells will not undergo haptotactic movement with probability
\begin{equation*} P_{HS_{i,j}}^k \;:\!=\; 1-\left ({P}_{HL_{i,j}}^k+{P}_{HR_{i,j}}^k\right ). \end{equation*}
2.1.3. Cell phenotypic changes
We model phenotypic changes by allowing cells to update their phenotypic states according to a random walk along the phenotypic dimension. Between time-steps
$k$
and
$k+1$
every cell enters a new phenotypic state with probability
$0\lt \beta \le 1$
, or remains in its current phenotypic state with probability
$1-\beta$
. Since we consider spontaneous phenotypic changes, we assume that a cell originally in phenotypic state
$y_j$
enters state
$y_{j-1}$
with probability
${P}_{D_{i,j}}^k$
or enters state
$y_{j+1}$
with probability
${P}_{U_{i,j}}^k$
, where we define
\begin{equation}{P}_{D_{i,j}}^k={P}_{U_{i,j}}^k \;:\!=\;\frac{\beta }{2}. \end{equation}
Therefore, cells will not undergo phenotypic changes with probability
\begin{equation*} {{P}_{N_{i,j}}^k\;:\!=\; 1-\left ({P}_{D_{i,j}}^k+{P}_{U_{i,j}}^k\right )}. \end{equation*}
Moreover, no-flux boundary conditions are implemented by aborting any attempted phenotypic change of a cell if it requires moving into a phenotypic state outside of the interval
$[0,Y]$
.
2.1.4. Cell division and death
To incorporate the effects of cell proliferation, we assume that a dividing cell is instantly replaced by two identical progeny cells that inherit the spatial position and phenotypic state of the parent cell. Conversely, a cell undergoing cell death is instantly removed from the population. To take into account phenotypic heterogeneity along with density-dependent inhibition of growth, at time-step
$t_k$
we assume that the probabilities of division and death for a cell at spatial position
$x_i$
depend both on the phenotypic state of the cell and the local cell density
$\rho _i^k$
.
In particular, to define the probabilities of cell division and death, we introduce the function
$R(y_j,\rho _i^k)$
, which describes the net growth rate of the cell population density at spatial position
$x_i$
and time
$t_k$
due to division and death of cells in the phenotypic state
$y_j$
, and assume that between time-steps
$k$
and
$k+1$
a cell in phenotypic state
$y_j$
at position
$x_i$
may die with probability
${P}_{A_{i,j}}^k$
, divide with probability
${P}_{B_{i,j}}^k$
, or remain quiescent with probability
${P}_{Q_{i,j}}^k \;:\!=\; 1-\left ({P}_{A_{i,j}}^k+{P}_{B_{i,j}}^k\right )$
, where
\begin{equation} {P}_{A_{i,j}}^k\;:\!=\;\tau \ R(y_j,\rho _i^k)_-, \quad{P}_{B_{i,j}}^k \;:\!=\;\tau \ R(y_j,\rho _i^k)_+, \quad \text{with} \quad (\!\cdot\!)_-\;:\!=\;-\!\min\!(0,\cdot ), \quad (\!\cdot\!)_+\;:\!=\;\max\!(0,\cdot ). \end{equation}
By considering the time-step
$\tau$
sufficiently small, we ensure
${P}_{A_{i,j}}^k+{P}_{B_{i,j}}^k\le 1$
for all values of
$i$
,
$j$
, and
$k$
We consider the scenario where: larger values of
$y_j$
correlate with a lower cell proliferation rate (cf. Figure 1); cells stop dividing if the cell density at their current position becomes larger than a critical value
$\rho _{\textrm{max}} \in \mathbb{R}^+_*$
. Therefore, we make the following assumptions
\begin{equation} R(Y,0)=0,\quad R(0,\rho _{\textrm{max}})=0,\quad \frac{\partial R(y,\rho )}{\partial \rho }\lt 0 \quad \text{and}\quad \frac{\partial R(y,\rho )}{\partial y}\lt 0\quad \text{for } (y,\rho )\in (0,Y)\times \mathbb{R}^+. \end{equation}
In particular, we focus on a similar case to that considered in [Reference Macfarlane, Ruan and Lorenzi46], that is, we assume
\begin{equation} R(y,\rho )\;:\!=\;\alpha \left (r(y)-\frac{\rho }{\rho _{\textrm{max}}}\right ), \quad r(Y)=0,\quad r(0)=1, \quad \text{and} \quad \frac{\mathrm{d}r(y)}{\mathrm{d}y}\lt 0\quad \text{for } y\in (0,Y), \end{equation}
with
$\alpha \in \mathbb{R}^+_*$
. Note that, under assumptions (9), the definitions given by (8) ensure that if
$\rho _i^k\ge \rho _{\textrm{max}}$
then every cell at position
$x_i$
can only die or remain quiescent between time-steps
$k$
and
$k+1$
. Hence, throughout the rest of the paper, we will assume
\begin{equation} \max _{i\in \mathbb{Z}}\ \rho _i^0\le \rho _{\textrm{max}} \end{equation}
so that
\begin{equation} \rho _i^k\le \rho _{\textrm{max}} \quad \text{for all }(k,i)\in \mathbb{N}_0\times \mathbb{Z}. \end{equation}
2.2. Modelling the dynamics of the MDEs and the ECM
The dynamics of the MDE concentration and the ECM density are governed by the rules provided in the following subsections, which are summarised by the schematics in Figure 2e-h and are coupled with the individual-based model for the dynamics of cancer cells that is presented in Section 2.1.
2.2.1. Dynamics of the MDE concentration
We let
$D_M \in \mathbb{R}^+_*$
be the diffusivity of the MDEs and we denote by
$\kappa _M \in \mathbb{R}^+_*$
the rate at which the MDEs undergo natural decay. To incorporate into the model the secretion of MDEs by cells in the phenotypic state
$y_j$
, we introduce the function
$p(y_j)$
. We focus on the scenario where larger values of
$y_j$
correlate with a higher MDE secretion rate (cf. Figure 1), i.e. we make the assumptions
\begin{equation} p(0)=p_{\textrm{min}} \in \mathbb{R}^+_*, \quad \frac{ \mathrm{d}p(y)}{\mathrm{d}y}\gt 0 \quad \text{for }y\in (0,Y). \end{equation}
In this framework, the principle of mass balance gives us the following difference equation for the concentration of MDEs
\begin{equation} M_i^{k+1}=M_{i}^k + \tau \left [ D_M\ (\mathscr{L}\ M^k)_i- \kappa _M \ M_i^k +\Delta _y \sum _{j} \left (p(y_j)\ n_{i,j}^k\right )\right ], \end{equation}
where
$\mathscr{L}$
is the finite-difference Laplacian on the lattice
$\{x_i\}_{i\in \mathbb{Z}}$
, that is,
\begin{equation*} (\mathscr {L}\ M^k)_i\;:\!=\;\frac {M_{i+1}^k+M_{i-1}^k-2\ M_i^k}{\Delta _x^2}.\end{equation*}
2.2.2. Dynamics of the ECM density
We denote by
$\kappa _E \in \mathbb{R}^+_*$
the rate at which the ECM is degraded by MDEs. The principle of mass balance gives us the following difference equation for the density of ECM
\begin{equation} E_i^{k+1}=E_i^k-\tau \ \kappa _E \ M_i^k \, E_i^k. \end{equation}
Recalling that, as mentioned earlier,
$E_{\textrm{max}} \in \mathbb{R}^+_*$
is the maximum value of the ECM density before cell invasion starts, we complement the difference equation (15) with an initial condition such that
\begin{equation} \max _{i\in \mathbb{Z}}\ E_i^0\le E_{\textrm{max}} \end{equation}
so that
\begin{equation} E_i^k \le E_{\textrm{max}} \quad \text{for all }(k,i)\in \mathbb{N}_0\times \mathbb{Z}. \end{equation}
3. The corresponding continuum model
Through an extension of the limiting procedure that we previously employed in [Reference Bubba, Lorenzi and Macfarlane12, Reference Chaplain, Lorenzi and Macfarlane16, Reference Macfarlane, Chaplain and Lorenzi45, Reference Macfarlane, Ruan and Lorenzi46], letting the time-step
$\tau \rightarrow 0$
, the space-step
$\Delta _x\rightarrow 0$
and the phenotype-step
$\Delta _y\rightarrow 0$
in such a way that
\begin{equation} \frac{\Delta _x^2 \theta }{2 \tau }\rightarrow D \in \mathbb{R}^+_{*}, \quad \frac{\Delta _x^2 \eta }{2 E_{\textrm{max}} \tau }\rightarrow \nu \in \mathbb{R}^+_*,\quad \text{and} \quad \frac{\Delta _y^2 \beta }{2\tau }\rightarrow \lambda \in \mathbb{R}^+_*, \end{equation}
and introducing the definition
\begin{equation} \chi (y) \;:\!=\; \nu \, \mu (y), \end{equation}
one can formally show (cf. Appendix A.1) that the deterministic continuum counterpart of the individual-based model presented in Section 2 is given by the PIDE (2)
$_1$
for the local cell population density function,
$n(t,x,y)$
, subject to zero Neumann (i.e. no-flux) boundary conditions at
$y=0$
and
$y=Y$
, complemented with the relation (2)
$_2$
for the cell density,
$\rho (t,x)$
and coupled with the PDE (2)
$_3$
for the MDE concentration,
$M(t,x)$
, along with the infinite-dimensional ODE (2)
$_4$
for the ECM density,
$E(t,x)$
. Consistently with assumptions (11) and (16), the PIDE (2)
$_1$
and the infinite-dimensional ODE (2)
$_4$
are subject to some initial conditions such that
\begin{equation} \max _{x \in \mathbb{R}} \int _0^Y n(0,x,y) \ \mathrm{d}y \le \rho _{\textrm{max}}, \quad \max _{x \in \mathbb{R}} E(0,x) \le E_{\textrm{max}}. \end{equation}
4. Formal asymptotic analysis
In this section, building on the formal asymptotic method that we developed in [Reference Lorenzi and Painter40, Reference Lorenzi, Perthame and Ruan41], which relies on the Hamilton-Jacobi approach developed in [Reference Barles, Mirrahimi and Perthame8, Reference Diekmann, Jabin, Mischler and Perthame20, Reference Lorz, Mirrahimi and Perthame44, Reference Perthame and Barles55], we carry out travelling wave analysis of an appropriately rescaled version of the model (2).
4.1. Rescaled model
We focus on a biological scenario wherein cell proliferation, cell production of MDEs and ECM degradation have a stronger impact on the dynamics of the system than haptotactic cell movement and diffusion of MDEs, which in turn have a stronger impact than random cell movement and cell phenotypic changes [Reference Anderson and Chaplain4, Reference Anderson, Chaplain, Newman, Steele and Thompson5, Reference Huang33, Reference Orsolits, Kovács, Kriston-Vizi, Merkely and Földes49, Reference Textor, Sinn and De Boer67]. To this end, we introduce a small parameter
$\varepsilon \in \mathbb{R}^+_*$
and choose the parameter scaling
\begin{equation} \nu \;:\!=\;{\varepsilon }, \quad D_M \;:\!=\;{\varepsilon }, \quad D \;:\!=\;{\varepsilon }^2, \quad \lambda \;:\!=\;{\varepsilon }^2. \end{equation}
Moreover, in order to explore the long-time behaviour of the system, we use the time scaling
$\displaystyle{t \to \frac{t}{\varepsilon }}$
in the model (2). In so doing, recalling the definition given by (19), we obtain the following rescaled system for the local cell population density function,
$n_{\varepsilon }(t,x,y)$
, the MDE concentration,
$M_{\varepsilon }(t,x)$
and the ECM density,
$E_{\varepsilon }(t,x)$
:
\begin{equation} \begin{cases} \displaystyle{{\varepsilon }\, \partial _t n_{\varepsilon } -{\varepsilon }\, \partial _x\Big ({\varepsilon }\, \partial _x n_{\varepsilon } - n_{\varepsilon } \, \mu (y) \, \partial _x E_{\varepsilon }\Big ) = R(y,\rho _{\varepsilon }) \, n_{\varepsilon } +{\varepsilon }^2 \, \partial _{yy} n_{\varepsilon }},\\[5pt] \displaystyle{\rho _{\varepsilon }(t,x) \;:\!=\; \int _0^Y n_{\varepsilon }(t,x,y)\ \mathrm{d}y},\\[5pt] \displaystyle{{\varepsilon }\, \partial _t M_{\varepsilon } -{\varepsilon }\, \partial _{xx}M_{\varepsilon } = \int _0^Y p(y) \, n_{\varepsilon }(t,x,y) \ \mathrm{d}y- \kappa _M \, M_{\varepsilon }},\\[5pt] \displaystyle{{\varepsilon }\, \partial _t E_{\varepsilon }=-\kappa _E \, E_{\varepsilon } \, M_{\varepsilon }}, \end{cases} \quad (x,y) \in \mathbb{R} \times (0,Y). \end{equation}
4.2. Formal limit for
${\varepsilon }\to 0$
We make the real phase WKB ansatz [Reference Barles, Evans and Souganidis7, Reference Evans and Souganidis22, Reference Fleming and Souganidis25]
\begin{equation} n_{\varepsilon }(t,x,y) = e^{\frac{u_{\varepsilon }(t,x,y)}{\varepsilon }}, \end{equation}
which gives
\begin{equation*} \partial _t n_{\varepsilon } = \frac {\partial _t u_{\varepsilon }}{\varepsilon } n_{\varepsilon }, \quad \partial _x n_{\varepsilon } = \frac {\partial _x u_{\varepsilon }}{\varepsilon } n_{\varepsilon }, \quad \partial ^2_{yy} n_{\varepsilon } = \left (\frac {1}{\varepsilon ^2} \left (\partial _y u_{\varepsilon } \right )^2 + \frac {1}{\varepsilon } \partial ^2_{yy} u_{\varepsilon } \right ) n_{\varepsilon }. \end{equation*}
Substituting the above expressions into the PIDE (22)
$_1$
for
$n_\varepsilon$
and rearranging terms gives the following Hamilton-Jacobi equation for
$u_{\varepsilon } \equiv u_{\varepsilon }(t,x,y)$
\begin{eqnarray*} \partial _t u_{\varepsilon } + \mu (y) \, \partial _{x} E_{\varepsilon } \, \partial _{x} u_{\varepsilon } &=& R(y,\rho _{\varepsilon }) + \left (\partial _x u_{\varepsilon } \right )^2 + \left (\partial _y u_{\varepsilon } \right )^2 + \\[5pt] && + \, \varepsilon \left (\partial ^2_{xx} u_{\varepsilon } - \mu (y) \, \partial ^2_{xx} E_{\varepsilon } + \partial ^2_{yy} u_{\varepsilon } \right ), \quad (x,y) \in \mathbb{R} \times (0,Y). \end{eqnarray*}
Now let
$\rho (t,x)$
be the leading-order term of the asymptotic expansion for
$\rho _{\varepsilon }(t,x)$
as
$\varepsilon \to 0$
. Considering
$x \in \mathbb{R}$
such that
$\rho (t,x) \gt 0$
(i.e.
$x \in \textrm{Supp}(\rho )$
) and letting
$\varepsilon \to 0$
in the above PDE we formally obtain the following equation for the leading-order term
$u \equiv u(t,x,y)$
of the asymptotic expansion for
$u_{\varepsilon }(t,x,y)$
\begin{equation} \partial _t u + \mu (y) \, \partial _{x} E \, \partial _{x} u = R(y,\rho ) + \left (\partial _x u \right )^2 + \left (\partial _y u \right )^2, \quad (x,y) \in{\textrm{Supp}(\rho )} \times (0,Y), \end{equation}
where
$E \equiv E(t,x)$
is the leading-order term of the asymptotic expansion for
$E_{\varepsilon }(t,x)$
.
Constraint on
$u$
. When
$\rho _{\varepsilon } \lt \infty$
for all
${\varepsilon }\in \mathbb{R}^+_*$
, if
$u_{\varepsilon }$
is a strictly concave function of
$y$
and
$u$
is also a strictly concave function of
$y$
whose unique maximum point is
$\bar{y}(t,x)$
then considering
$x \in \textrm{Supp}(\rho )$
and letting
$\varepsilon \to 0$
in (23) formally gives the following constraint on
$u$
\begin{equation} u(t,x,\bar{y}(t,x)) = \max _{y \in [0,Y]} u(t,x,y) = 0, \quad x \in \textrm{Supp}(\rho ), \end{equation}
which implies that
\begin{equation} \partial _y u(t,x,\bar{y}(t,x)) = 0 \quad \text{and} \quad \partial _x u(t,x,\bar{y}(t,x)) = 0, \quad x \in \textrm{Supp}(\rho ). \end{equation}
Relation between
$\bar{y}(t,x)$
and
$\rho (t,x)$
. Evaluating (24) at
$y=\bar{y}(t,x)$
and using (25) and (26), we find
\begin{equation*} R(\bar {y}(t,x),\rho (t,x)) = 0, \quad x \in \textrm {Supp}(\rho ), \end{equation*}
from which, using the fact that the function
$R(y,\rho )$
is defined via (10), we obtain the following formula
\begin{equation} \rho (t,x) = \rho _{\textrm{max}} \, r(\bar{y}(t,x)), \quad x \in \textrm{Supp}(\rho ). \end{equation}
The monotonicity assumption (10) ensures that (27) gives a one-to-one correspondence between
$\bar{y}(t,x)$
and
$\rho (t,x)$
.
Expressions of
$M(t,x)$
and
$E(t,x)$
. When
$n_{\varepsilon }$
is in the form (23), if
$u_{\varepsilon }$
is a strictly concave function of
$y$
and
$u$
is also a strictly concave function of
$y$
that satisfies the constraint (25) then the following asymptotic result formally holds
where
$\delta _{\bar{y}(t,x)}(y)$
is the Dirac delta centred at
$y=\bar{y}(t,x)$
. In this case, focussing on a biological scenario wherein the ECM density is at the maximum level
$E_{\textrm{max}}$
before cell invasion starts at
$t=0$
, letting
$\varepsilon \to 0$
in the PDE (22)
$_3$
for
$M_\varepsilon$
and in the infinite-dimensional ODE (22)
$_4$
for
$E_\varepsilon$
we formally obtain the following expressions of the leading-order terms of the asymptotic expansions for
$M_{\varepsilon }(t,x)$
and
$E_{\varepsilon }(t,x)$
where 
denotes the indicator function of the set
$(\!\cdot\!)$
.
Remark 1. Note that the behaviour of
$E(t,x)$
depicted by (28) shares similarities with the behaviour of the nutrient concentration in the model analysed in [Reference Jabin and Perthame35].
Transport equation for
$\bar{y}$
. When
$R(y,\rho )$
is defined via (10), differentiating (24) with respect to
$y$
, evaluating the resulting equation at
$y=\bar{y}(t,x)$
and using (25) and (26) yields
\begin{equation} \partial ^2_{yt} u(t,x,\bar{y}) + \mu (\bar{y}) \, \partial _{x} E \, \partial ^2_{yx} u(t,x,\bar{y}) = \partial _{y} r(\bar{y}), \quad x\in \textrm{Supp}(\rho ). \end{equation}
Moreover, differentiating (26) with respect to
$t$
and
$x$
we find, respectively,
\begin{equation*} \partial ^2_{ty} u(t,x,\bar {y}) + \partial ^2_{yy} u(t,x,\bar {y}) \, \partial _{t} \bar {y}(t,x) = 0 \; \Rightarrow \; \partial ^2_{yt} u(t,x,\bar {y}) = - \partial ^2_{yy} u(t,x,\bar {y}) \, \partial _{t} \bar {y}(t,x), \quad x\in \textrm {Supp}(\rho ) \end{equation*}
and
\begin{equation*} \partial ^2_{xy} u(t,x,\bar {y}) + \partial ^2_{yy} u(t,x,\bar {y}) \, \partial _{x} \bar {y}(t,x) = 0 \; \Rightarrow \; \partial ^2_{yx} u(t,x,\bar {y}) = - \partial ^2_{yy} u(t,x,\bar {y}) \, \partial _{x} \bar {y}(t,x), \quad x\in \textrm {Supp}(\rho ). \end{equation*}
Substituting the above expressions of
$\partial ^2_{yt} u(t,x,\bar{y})$
and
$\partial ^2_{yx} u(t,x,\bar{y})$
into (29), and using the fact that if
$u$
is a strictly concave function of
$y$
whose unique maximum point is
$\bar{y}(t,x)$
then
$\partial ^2_{yy} u(t,x,\bar{y}) \lt 0$
, gives the following transport equation for
$\bar{y}(t,x)$
\begin{equation} \partial _{t} \bar{y} + \mu (\bar{y}) \, \partial _{x} E \, \partial _{x} \bar{y} = \frac{\partial _{y} r(\bar{y})}{-\partial ^2_{yy} u(t,x,\bar{y})}, \quad x \in \textrm{Supp}(\rho ). \end{equation}
4.3. Travelling wave analysis
4.3.1. Travelling wave problem
Substituting the travelling-wave ansatz
\begin{equation*} u(t,x,y) = u(z,y), \quad \rho (t,x) = \rho (z), \quad \bar {y}(t,x) = \bar {y}(z), \quad M(t,x) = M(z), \quad \text {and} \quad E(t,x) = E(z) \end{equation*}
with
\begin{equation*} z = x - c \, t, \quad c \in \mathbb {R}^+_*, \end{equation*}
into (27), (28), and (30) gives
\begin{equation} \rho (z) = \rho _{\textrm{max}} \, r(\bar{y}(z)), \quad z \in \textrm{Supp}(\rho ), \end{equation}
\begin{equation} \left (c - \mu (\bar{y}) E' \right ) \bar{y}' = \frac{\partial _{y} r(\bar{y})}{\partial ^2_{yy} u(z,\bar{y})}, \quad z \in \textrm{Supp}(\rho ). \end{equation}
The expression (32) of
$E(z)$
makes it possible to simplify the differential equation (33) as follows
\begin{equation} c \, \bar{y}' = \frac{\partial _{y} r(\bar{y})}{\partial ^2_{yy} u(z,\bar{y})}, \quad z \in \textrm{Supp}(\rho ). \end{equation}
We complement the differential equation (34) with the following asymptotic condition
\begin{equation} \lim _{z \to - \infty } \bar{y}(z) =0, \end{equation}
so that, since
$r(0)=1$
(cf. assumptions (10)), the relation (31) gives
\begin{equation} \lim _{z \to - \infty } \rho (z) = \rho _{\textrm{max}}. \end{equation}
4.3.2. Shape of travelling waves
Since
$\partial _y r(y)\lt 0$
for
$y \in (0,Y]$
(cf. assumptions (10)) and given the fact that if
$u$
is a strictly concave function of
$y$
whose unique maximum point is
$\bar{y}(t,x)$
then
$\partial ^2_{yy} u(t,x,\bar{y}) \lt 0$
, the differential equation (34) along with the relation (31) ensure that
\begin{equation} \bar{y}'(z) \gt 0 \quad \text{and} \quad \rho '(z) \lt 0, \quad z \in \textrm{Supp}(\rho ). \end{equation}
The relation (31) and the monotonicity results (37) along with the fact that
$r(Y)=0$
(cf. assumptions (10)) imply that the position of the edge of the travelling front
$\rho (z)$
coincides with the unique point
$\ell \in \mathbb{R}$
such that
$\bar{y}(\ell )=Y$
and
$\bar{y}(z) \lt Y$
on
$(\!-\!\infty, \ell )$
. Hence,
$\textrm{Supp}(\rho ) = (\!-\!\infty, \ell )$
and, since
$p(y)\gt 0$
for all
$y \in [0,Y]$
(cf. assumptions (13)), the expressions (32) of
$M(z)$
and
$E(z)$
yield
5. Numerical simulations
In this section, we report on numerical solutions of the rescaled continuum model (22) and numerical simulations of the corresponding rescaled version of the individual-based model, and we compare them with the results of travelling wave analysis presented in the previous section.
5.1. Set-up of numerical simulations
We start by describing the set-up used to carry out numerical simulations.
5.1.1. Model functions and parameters
To allow the individual-based model to represent the same scenario as the rescaled continuum model (22), we use the same time scaling
$t_k\ \rightarrow \ \frac{t_k}{\varepsilon }=k\frac{\tau }{\varepsilon }$
and reformulate the governing rules for the cell dynamics detailed in Section 2 in terms of
\begin{equation*} n_{\varepsilon _{i,j}}^k\equiv n_\varepsilon (t_k,x_i,y_j)=n\left (\frac {t_k}{\varepsilon },x_i,y_j\right )\;:\!=\;\frac {N_{\varepsilon _{i,j}}^k}{\Delta _x \Delta _y}, \quad \rho _{\varepsilon _{i}}^k\equiv \rho _\varepsilon (t_k,x_i)=\rho \left (\frac {t_k}{\varepsilon },x_i\right )\;:\!=\; \Delta _y \sum _j n_{\varepsilon _{i,j}}^k, \end{equation*}
\begin{equation*} M_{\varepsilon _{i}}^k\equiv M_\varepsilon (t_k,x_i)=M\left (\frac {t_k}{\varepsilon },x_i\right ),\quad E_{\varepsilon _{i}}^k\equiv E_\varepsilon (t_k,x_i)=E\left (\frac {t_k}{\varepsilon },x_i\right ). \end{equation*}
To ensure that conditions (18) and (21) are simultaneously satisfied, we additionally set
\begin{equation*} \theta =\frac {2\tau }{\Delta _x^2}\varepsilon ^2, \quad \eta =\frac {2E_{\textrm {max}}\tau }{\Delta _x^2}\varepsilon, \quad \beta =\frac {2\tau }{\Delta _y^2}\varepsilon ^2. \end{equation*}
In order to carry out numerical simulations, we consider the time interval
$[0,T]$
with
$T=30$
. Furthermore, we restrict the physical domain to the interval
$[0,X]$
, with
$X=100$
, and choose
$Y=1$
. Moreover, we specifically choose
$\Delta _x=5 \times 10^{-2}$
,
$\Delta _y=2 \times 10^{-2}$
and
$\tau = \dfrac{\Delta _x^2}{2}$
.
To satisfy assumptions (5), we use the definition
\begin{equation} \mu (y)\;:\!=\;y^2, \end{equation}
while in order to satisfy assumptions (9), we define
$R(y,\rho )$
via (10) with
$\alpha =0.1$
and, having chosen
$Y=1$
, we further define
\begin{equation} r(y) \;:\!=\;1-y^2. \end{equation}
To satisfy assumptions (13) on
$p(y)$
, we also define
\begin{equation} p(y)\;:\!=\;p_{\textrm{min}}+ \zeta y^2, \end{equation}
where
$\zeta =10^{-5}$
and
$p_{\textrm{min}}=10^{-7}$
. Furthermore, in the simulations we choose
$\kappa _M=1$
,
$\kappa _E=1$
and
$E_{\rm{max}}=1$
.
5.1.2. Initial conditions
We consider a biological scenario in which, initially, the cell population is localised along the
$x=0$
boundary of the spatial domain and most of the cells are in the phenotypic state
$y=\bar{y}^0$
at every position
$x \in [0,X]$
. Specifically, we implement the following initial cell distribution for the IB model
\begin{equation} N_{i,j}^0=\text{int}(F(x_i,y_j)) \; \text{ with } \; F(x,y)\;:\!=\; A_0 \ C\ e^{-x^2} \ e^{-\frac{\left (y-\bar{y}^0\right )^2}{\varepsilon }}, \end{equation}
where
$\text{int}(\!\cdot\!)$
is the integer part of
$(\!\cdot\!)$
and
$C$
is a normalisation constant such that
\begin{equation*} C\int _0^Y e^{-\frac {\left (y-\bar {y}^0\right )^2}{\varepsilon }} \ \mathrm {d}y=1. \end{equation*}
We choose
$\bar{y}^0=0.2$
and
$A_0=100$
. The initial cell density
$\rho _{i}^0$
is then calculated from (42) according to the definition given by (3), and we set
$\displaystyle{\rho _{\textrm{max}}=\max _{i} \rho ^0_i}$
.
Moreover, we assume that there are initially no MDEs and the density of ECM is uniform, that is,
\begin{equation*} M^0_i=0 \quad \text {and} \quad E^0_i=E_{\rm {max}}\quad \text {for all } i. \end{equation*}
Finally, we consider different values of
$\varepsilon$
, that is,
${\varepsilon }\in \left \{10^{-2}, 5 \times 10^{-3}, 10^{-3}\right \}$
, in order to verify whether, for
$\varepsilon$
small enough, there is a good agreement between the results of numerical simulations and the results of formal asymptotic analysis for
${\varepsilon }\to 0$
presented in Section 4.
5.2. Computational implementation of the individual-based model
All simulations of the individual-based model were performed in Matlab, and the numerical scheme used to solve the rescaled system (22) was also implemented in Matlab.
Figure 3. Numerical simulation results of the individual-based model (top row) and numerical solutions of the corresponding continuum model (22) (bottom row) in the case where
$\varepsilon =10^{-2}$
. The plots display, from left to right, the cell population density,
$n_\varepsilon$
, the cell density,
$\rho _\varepsilon$
, the MMP concentration,
$M_\varepsilon$
, and the ECM density,
$E_\varepsilon$
, at progressive times (i.e.
$t = 10$
,
$t = 20$
and
$t = 30$
) for both modelling approaches. Top row. The results from the individual-based model were obtained by averaging over five simulations (solid blue lines), and we additionally plot the corresponding results of each simulation (solid cyan lines). We also include the equivalent numerical solutions of the continuum model (dotted red lines) for comparison. Bottom row. The values of
$\rho _\varepsilon$
,
$M_\varepsilon$
and
$E_\varepsilon$
(solid red lines) are plotted along with the corresponding values obtained by substituting
$\bar{y}(t,x)=\bar{y}_{\varepsilon }(t,x)$
into the formulas given by (27) and (28) (black dotted lines), with
$\bar{y}_{\varepsilon }(t,x)$
being the maximum point of the numerical solution
$n_\varepsilon (t,x,y)$
to the PIDE (22)
$_1$
at position
$x$
at time
$t$
.
In the individual-based model, at each time-step, every single cell undergoes a four-step process: (i) random cell movement, according to the probabilities defined via (4); (ii) haptotactic cell movement, according to the probabilities defined via (6); (iii) phenotypic changes, according to the probabilities defined via (7); and (iv) cell division and death, according to the probabilities defined via (8). For every single cell, during each step of this process, a random number is drawn from the standard uniform distribution on the interval
$(0,1)$
using the built-in Matlab function rand. It is then evaluated whether this number is lower than the probability of the event occurring and if so the event occurs. Since to carry out numerical simulations we have to restrict the spatial domain to the interval
$[0, X]$
, the attempted movement of a cell is aborted if it requires moving out of this interval. Furthermore, the concentration of MDEs and the density of ECM are calculated using the discrete difference equations (14) and (15), respectively.
5.3. Numerical methods for the continuum model
To solve numerically the rescaled system (22) posed on
$(0, T] \times (0,X) \times (0,Y)$
subject to zero-flux boundary conditions and complemented with the continuum analogues of the initial conditions selected for the individual-based model, we employ a finite volume scheme modified from our previous work [Reference Bubba, Lorenzi and Macfarlane12].
5.4. Main results of numerical simulations
Our main results of the numerical simulations of the individual-based model and the corresponding numerical solutions of the continuum model (22) for three distinct values of the scaling parameter
$\varepsilon$
are summarised by the plots in Figures 3-5, which correspond to
$\varepsilon =10^{-2}$
,
$\varepsilon =5 \times 10^{-3}$
and
$\varepsilon =10^{-3}$
, respectively.
Figure 4. Numerical simulation results of the individual-based model (top row) and numerical solutions of the corresponding continuum model (22) (bottom row) in the case where
$\varepsilon =5 \times 10^{-3}$
. The plots display, from left to right, the cell population density,
$n_\varepsilon$
, the cell density,
$\rho _\varepsilon$
, the MMP concentration,
$M_\varepsilon$
and the ECM density,
$E_\varepsilon$
, at progressive times (i.e.
$t = 10$
,
$t = 20$
and
$t = 30$
) for both modelling approaches. Top row. The results from the individual-based model were obtained by averaging over five simulations (solid blue lines), and we additionally plot the corresponding results of each simulation (solid cyan lines). We also include the equivalent numerical solutions of the continuum model (dotted red lines) for comparison. Bottom row. The values of
$\rho _\varepsilon$
,
$M_\varepsilon$
and
$E_\varepsilon$
(solid red lines) are plotted along with the corresponding values obtained by substituting
$\bar{y}(t,x)=\bar{y}_{\varepsilon }(t,x)$
into the formulas given by (27) and (28) (black dotted lines), with
$\bar{y}_{\varepsilon }(t,x)$
being the maximum point of the numerical solution
$n_\varepsilon (t,x,y)$
to the PIDE (22)
$_1$
at position
$x$
at time
$t$
.
Figure 5. Numerical simulation results of the individual-based model (top row) and numerical solutions of the corresponding continuum model (22) (bottom row) in the case where
$\varepsilon =10^{-3}$
. The plots display, from left to right, the cell population density,
$n_\varepsilon$
, the cell density,
$\rho _\varepsilon$
, the MMP concentration,
$M_\varepsilon$
and the ECM density,
$E_\varepsilon$
, at progressive times (i.e.
$t = 5$
,
$t = 10$
and
$t = 15$
) for both modelling approaches. Top row. The results from the individual-based model were obtained by averaging over five simulations (solid blue lines), and we additionally plot the corresponding results of each simulation (solid cyan lines). We also include the equivalent numerical solutions of the continuum model (dotted red lines) for comparison. Bottom row. The values of
$\rho _\varepsilon$
,
$M_\varepsilon$
and
$E_\varepsilon$
(solid red lines) are plotted along with the corresponding values obtained by substituting
$\bar{y}(t,x)=\bar{y}_{\varepsilon }(t,x)$
into the formulas given by (27) and (28) (black dotted lines), with
$\bar{y}_{\varepsilon }(t,x)$
being the maximum point of the numerical solution
$n_\varepsilon (t,x,y)$
to the PIDE (22)
$_1$
at position
$x$
at time
$t$
.
The plots in the top rows of Figures 3
–5 are the results of the individual-based model averaged over 5 simulations. In particular, from left to right, we have the cell population density,
$n_{\varepsilon _{i,j}}^k$
, the cell density,
$\rho _{\varepsilon _{i}}^k$
, the MDE concentration,
$M_{\varepsilon _{i}}^k$
and the ECM density,
$E_{\varepsilon _{i}}^k$
, at progressive times. On the other hand, the plots in the bottom rows of Figures 3
–5 are the numerical solutions of the continuum model, which are plotted along with the corresponding analytical results presented in Section 4.
These plots show that there is a good agreement between the results of numerical simulations of the individual-based model and numerical solutions of the continuum model. This validates the limiting procedure that we employed to formally derive the continuum model. The same plots also demonstrate that the smaller is the value of
$\varepsilon$
, then the better the agreement between numerical solutions of the rescaled continuum model and the analytical results presented in Section 4. This validates the formal asymptotic method that we used to construct, in the limit as
${\varepsilon } \to 0$
, invading fronts with spatial structuring of cell phenotypes. In particular, the plots in Figure 5 demonstrate that, when
$\varepsilon$
is sufficiently small:
-
(i) The local cell population density function
$n_{\varepsilon }(t,x,y)$
becomes concentrated as a sharp Gaussian with maximum at a point
$\bar{y}_{\varepsilon }(t,x)$
for all
$x$
where
$\rho _{\varepsilon }(t,x)\gt 0$
. -
(ii) The maximum point
$\bar{y}_{\varepsilon }(t,x)$
behaves like a compactly supported and monotonically increasing travelling front that connects
$0$
to
$Y$
– recall that here
$Y=1$
. This indicates that cells in phenotypic states
$y\approx Y$
are concentrated towards the leading edge of the invading front, while cells in phenotypic states corresponding to smaller values of
$y$
make up the bulk of the population in the rear. -
(iii) The cell density
$\rho _{\varepsilon }(t,x)$
behaves like a one-sided compactly supported and monotonically decreasing travelling front that connects
$\rho _{\textrm{max}}$
to
$0$
. -
(iv) The values of
$\rho _\varepsilon$
,
$M_\varepsilon$
and
$E_\varepsilon$
are consistent with the values obtained by substituting
$\bar{y}(t,x)=\bar{y}_{\varepsilon }(t,x)$
into the formulas given by (27) and (28).
6. Conclusion
We have formulated a model for cancer invasion in which the infiltrating cancer cells can occupy a spectrum of states in the phenotype space, ranging from ‘fully mesenchymal’ to ‘fully epithelial’. More precisely, the more mesenchymal cells are those that display stronger haptotaxis responses and have greater capacity to modify the ECM through enhanced secretion of MDEs. However, as a trade-off, they have lower proliferative capacity than the more epithelial cells. The framework is multiscale in that we start with an individual-based model that tracks the dynamics of single cells and is based on a branching random walk over a lattice, where cell movements take place through both physical and phenotype space. By applying limiting techniques, we have formally derived the corresponding continuum model, which takes the form of system (2). Despite the intricacy of the model, we showed, through formal asymptotic techniques, that for certain parameter regimes it is possible to carry out a detailed travelling wave analysis and obtain the form for the wave profile. Simulations have been performed of both the individual-based model and the continuum model, which generally show excellent correspondence. Moreover, when parameter values are chosen from the appropriate parameter regime, numerical solutions to the continuum model match closely with the corresponding analytical form. This validates the formal limiting procedure that we employed to derive the continuum model alongside the formal asymptotic method that we used to characterise the wave profile.
Notably, solutions to the model reveal a capacity for self-organisation, in the sense that an initially almost homogeneous population resolves itself into an invading front with spatial structuring of phenotypes. Precisely, the most mesenchymal cells dominate the leading edge of the invasion wave and the most epithelial (and most proliferative) dominate the rear, representing a bulk tumour population. As such, the model recapitulates similar observations into a front to back structuring of invasion waves into leader-type and follower-type cells, witnessed in an increasing number of experimental studies over recent years [Reference Vilchez Mercedes, Bocci, Levine, Onuchic, Jolly and Wong68].
A number of other continuum models have been formulated to study how phenotypic diversity alters cancer invasion processes. These include those intended to describe “go-or-grow” dynamics, a term coined for glioma growth processes where a dichotomy of cells into proliferating or migrating cell types has been suggested [Reference Giese, Loo, Tran, Haskett, Coons and Berens30]. In many of these models, heterogeneity is restricted to binary states (a proliferating class and a migrating class), with functions defining the state to state transitions. Often these models have restricted to relatively simple assumptions for cell migration processes, e.g. a simple diffusion process [Reference Pham, Chauviere and Hatzikirou57, Reference Stepien, Rutter and Kuang63], although more complex movement models have also been considered, e.g. [Reference Conte and Surulescu19]. The model here expands the potential framework for modelling go-or-growth processes, to cover all potential states between fully proliferative and fully migratory. Binary cell state models with distinct proliferative and migratory characteristics have also been developed to describe acid-mediated cancer invasion [Reference Strobl, Krause, Damaghi, Gillies, Anderson and Maini65], where a similar wave structuring of the distinct populations can be found under certain configurations. Invasion models with continuous phenotypes that range from more migratory to more proliferative states have been applied to avascular cancer growth, where further variables are included to describe tissue oxygen levels [Reference Fiandaca, Bernardi, Scianna and Delitala23]. The study here provides a structure for more detailed analytical studies into the travelling wave dynamics observed in these primarily numerically-based investigations. While not specifically focussing on cancer invasion, continuous phenotypic structuring models have also been developed and analysed for chemotaxis-driven wave invasion in [Reference Lorenzi and Painter40] and density- or pressure-driven wave invasion in [Reference Lorenzi, Perthame and Ruan41, Reference Macfarlane, Ruan and Lorenzi46]; the study here expands on the methodologies introduced therein, reinforcing their utility to study a diverse range of models used to explain invasion processes.
There are, clearly, various further extensions that could be considered. The current study has concentrated on invasion processes in 1D for each of physical space (i.e. a transect across the invading front) and phenotype space (from epithelial-like to mesenchymal-like). It is possible, of course, to extend the dimensionality of either or both of these spaces. As a way of illustration, preliminary simulations are presented in Figure 6 for an extension to 2D for the physical space in the individual-based model, where from an initial population concentrated at the origin we observe the emergence of a quasi-symmetric growing tumour with leader-to-follower structuring across the radial transect (consistent with the corresponding 1D model). A natural question to explore in two dimensions would be whether there are conditions under which the symmetric growth breaks, e.g. whether spatial structuring such as ‘tumour fingering’ emerges. Previous models have shown that fingering can occur in various scenarios, such as tumour infiltration into heterogeneous ECM environments [Reference Anderson, Weaver, Cummings and Quaranta6] and tumour growth in the presence of cells with different mobilities [Reference Drasdo and Hoehme21, Reference Lorenzi, Lorz and Perthame43]. Whether it is also possible for such phenomena to develop under certain assumptions for the manner in which phenotypic transitions occur could be a point of focus.
Figure 6. Preliminary 2D results from the individual-based model, obtained by averaging over 15 simulations, in the case where
$\varepsilon =10^{-2}$
. Top and middle rows. The plots display, from left to right, the maximum point of the cell population density,
$\bar{y}_{\varepsilon } = \mathrm{arg\,max}_{y \in [0,Y]} n_\varepsilon$
, the cell density,
$\rho _\varepsilon$
, the MMP concentration,
$M_\varepsilon$
, and the ECM density,
$E_\varepsilon$
, at the start and end of simulations – i.e.
$t = 0$
(top) and
$t = 5$
(middle). Bottom row. The plots display, from left to right,
$\displaystyle{\bar{y}_{\varepsilon }}$
,
$\rho _\varepsilon$
,
$M_\varepsilon$
and
$E_\varepsilon$
across the radial transect at the end of simulations (i.e. at
$t = 5$
). Here
$x\equiv (x_1,x_2) \in [0,10] \times [0,10]$
with
$\Delta _{x_1}=\Delta _{x_2}=0.1$
,
$y \in [0,1]$
with
$\Delta _y=0.02$
, and the initial cell distribution is the 2D analogue of (42) with
$A_0=1$
, while all the other parameters and functions are kept the same as in the 1D simulations of Figures 3
–5.
Extensions in the dimensionality of the phenotype space may also be of interest. Here we have considered a linear pathway from the fully epithelial state to the fully mesenchymal state – i.e. where reduced proliferation is accompanied simultaneously by an upregulation in both haptotaxis and secretion of MDEs. Given that EMTs in cancers can often be ‘partial’ in nature [Reference Vilchez Mercedes, Bocci, Levine, Onuchic, Jolly and Wong68], it is possible that different pathways could be taken from epithelial to mesenchymal: cells could follow separate pathways in which first haptotactic movement is upregulated and then secretion of MDEs, or vice versa. Extensions to a higher dimensionality in the phenotype space would allow exploration into whether this can give rise to additional subtlety in the positioning of different phenotypes.
Another natural extension would be to target the model towards particular experimental studies of leader–follower behaviour, by incorporating system-specific phenomena. For example, studies of collective invasion in non-small cell lung cancer (NSCLC) tumour spheroids have led to an experimental model of invasion with intricate signalling between leader and follower cells [Reference Konen, Summerbell and Dwivedi38]: in terms of movement, leader cells secrete fibronectin and release VEGF that guides follower cell movements through chemotaxis; in terms of proliferation, followers secrete factors to promote leader cell proliferation, while leaders secrete factors that hinder follower growth. Adapting the model to include these additional factors and their impact on follower/leader behaviour would provide a means to test the experimental model and, for example, investigate how perturbations to various aspects of the signalling system would impact on the rate of infiltration.
Further exploration could be made into the processes that lead to phenotypic changes. Here, we have adopted the relatively simple assumption of (unbiased) random phenotype switches, which at the cell-population level leads to diffusion across phenotype space. It is also quite plausible that phenotype changes may be biased in particular directions – e.g. from epithelial towards mesenchymal (or vice versa) – and that the direction and strength of the bias changes with the tumour microenvironment [Reference Aggarwal, Montoya, Donnenberg and Sant1]. Our model has also decoupled proliferation from phenotypic changes, effectively assuming that proliferative events lead to daughter cells of the same phenotype; divisions may also occur in an asymmetric manner – e.g. division of a follower cell leading to a daughter of leader phenotype. With our framework, the impact of such changes can be investigated both at the individual and continuous level.
Summarising, the work here provides the framework for developing and analysing sophisticated haptotaxis models for cancer invasion in which the cell population contains significant phenotypic heterogeneity. While these models present significant challenges at both a numerical and analytical level, we believe the methods developed and described here can allow for further progress in this area.
Acknowledgments
TL and KJP would also like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme Mathematics of Movement, where work on this paper was undertaken.
Financial support
TL gratefully acknowledges support from the Italian Ministry of University and Research (MUR) through the grant PRIN 2020 project (No. 2020JLWP23) “Integrated Mathematical Approaches to Socio-Epidemiological Dynamics” (CUP: E15F21005420006) and the grant PRIN2022-PNRR project (No. P2022Z7ZAJ) “A Unitary Mathematical Framework for Modelling Muscular Dystrophies” (CUP: E53D23018070001), from the CNRS International Research Project “Modelisation de la biomecanique cellulaire et tissulaire” (MOCETIBI), and from the Istituto Nazionale di Alta Matematica (INdAM) and the Gruppo Nazionale per la Fisica Matematica (GNFM).
KJP is a member of INdAM-GNFM and acknowledges “Miur-Dipartimento di Eccellenza” funding to the Dipartimento di Scienze, Progetto e Politiche del Territorio (DIST).
Competing interests
The authors declare that they have no competing interests.
Appendix A
A.1 Formal derivation of the continuum model (2)
We provide the full details of the formal derivation of the continuum model (2) from the individual-based model presented in Section 2, which relies on an extension of the limiting procedure that we previously employed in [Reference Bubba, Lorenzi and Macfarlane12, Reference Chaplain, Lorenzi and Macfarlane16, Reference Macfarlane, Chaplain and Lorenzi45, Reference Macfarlane, Ruan and Lorenzi46].
A.1.1 Equation for cell population density
When cell dynamics are governed by the rules underlying the individual-based model presented in Section 2, the principle of mass balance gives, for the cell population density,
\begin{eqnarray} n^{k+1}_{i,j}&=&n^{k}_{i+1,j+1}\left \{ \frac{\beta }{2}\left [ 1+\tau R(y_j, \rho ^{k}_{i})\right ]\left [\frac{\theta }{2}+\frac{\eta \mu (y_{j})}{2 E_{\textrm{max}}} \left (E^{k}_{i}-E^{k}_{i+1}\right )_+\right ]\right \}\nonumber\\[5pt] &&+n^{k}_{i-1,j+1}\left \{ \frac{\beta }{2}\left [ 1+\tau R(y_j, \rho ^{k}_{i})\right ]\left [\frac{\theta }{2}+\frac{\eta \mu (y_{j})}{2 E_{\textrm{max}}} \left (E^{k}_{i}-E^{k}_{i-1}\right )_+\right ]\right \}\nonumber \\[5pt] &&+n^{k}_{i+1,j-1}\left \{ \frac{\beta }{2}\left [ 1+\tau R(y_j, \rho ^{k}_{i})\right ]\left [\frac{\theta }{2}+\frac{\eta \mu (y_{j})}{2 E_{\textrm{max}}} \left (E^{k}_{i}-E^{k}_{i+1}\right )_+\right ]\right \}\nonumber \\[5pt] &&+n^{k}_{i-1,j-1}\left \{ \frac{\beta }{2}\left [ 1+\tau R(y_j, \rho ^{k}_{i})\right ]\left [\frac{\theta }{2}+\frac{\eta \mu (y_{j})}{2 E_{\textrm{max}}} \left (E^{k}_{i}-E^{k}_{i-1}\right )_+\right ]\right \}\nonumber \\[5pt] &&+n^{k}_{i,j+1}\left \{ \frac{\beta }{2}\left [ 1+\tau R(y_j, \rho ^{k}_{i})\right ]\left [1-\theta -\frac{\eta \mu (y_{j})}{2 E_{\textrm{max}}}\left [ \left (E^{k}_{i+1}-E^{k}_{i}\right )_+ + \left (E^{k}_{i-1}-E^{k}_{i}\right )_+\right ]\right ] \right \}\nonumber \\[5pt] &&+n^{k}_{i,j-1}\left \{ \frac{\beta }{2}\left [ 1+\tau R(y_j, \rho ^{k}_{i})\right ]\left [1-\theta -\frac{\eta \mu (y_{j})}{2 E_{\textrm{max}}}\left [ \left (E^{k}_{i+1}-E^{k}_{i}\right )_+ + \left (E^{k}_{i-1}-E^{k}_{i}\right )_+\right ]\right ] \right \}\\[5pt] &&+n^{k}_{i+1,j}\left \{ (1-\beta )\left [ 1+\tau R(y_j, \rho ^{k}_{i})\right ]\left [\frac{\theta }{2}+\frac{\eta \mu (y_{j})}{2 E_{\textrm{max}}} \left (E^{k}_{i}-E^{k}_{i+1}\right )_+\right ]\right \}\nonumber \\[5pt] &&+n^{k}_{i-1,j}\left \{ (1-\beta )\left [ 1+\tau R(y_j, \rho ^{k}_{i})\right ]\left [\frac{\theta }{2}+\frac{\eta \mu (y_{j})}{2 E_{\textrm{max}}} \left (E^{k}_{i}-E^{k}_{i-1}\right )_+\right ]\right \}\nonumber \\[5pt] &&+n^{k}_{i,j}\left \{ (1-\beta )\left [ 1+\tau R(y_j, \rho ^{k}_{i})\right ] \left [1-\theta -\frac{\eta \mu (y_{j})}{2 E_{\textrm{max}}}\left [\left (E^{k}_{i+1}-E^{k}_{i}\right )_++\left (E^{k}_{i-1}-E^{k}_{i}\right )_+ \right ]\right ]\right \}.\nonumber \end{eqnarray}
Using the fact that for
$\tau$
,
$\Delta _x$
and
$\Delta _y$
sufficiently small the following relations hold
\begin{align*} n^{k}_{i,j}\approx n(t,x,y)\equiv n,\;\;\; n^{k+1}_{i,j}\approx n(t+\tau,x,y),\;\;\; n^{k}_{i\pm 1,j}\approx n(t,x\pm \Delta _x,y),\;\;\; n^{k}_{i,j\pm 1}\approx n(t,x,y\pm \Delta _y)&\\[5pt] E^{k}_{i}\approx E(t,x)\equiv E, \;\;\; E^{k}_{i\pm 1}\approx E(t,x\pm \Delta _x), \;\;\; \mu (y_j)\approx \mu (y) \equiv \mu, \;\;\; \rho ^{k}_{i} \approx \rho (t,x)\equiv \rho, \;\;\; R(y_j,\rho ^{k}_{i})\approx R(y,\rho ) \equiv R,& \end{align*}
Equation (A1) can be rewritten as
\begin{align} &n(t+\tau,x,y)\nonumber \\[5pt] & \quad = n(t,x+\Delta _x,y+\Delta _y)\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\qquad +n(t,x-\Delta _x,y+\Delta _y)\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\qquad +n(t,x+\Delta _x,y-\Delta _y)\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\qquad +n(t,x-\Delta _x,y-\Delta _y)\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \} \\[5pt] &\qquad +n(t,x,y+\Delta _y)\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [1-\theta -\frac{\eta \mu }{2 E_{\textrm{max}}}\left [ \left (E(t,x+\Delta _x)-E\right )_+ + \left (E(t,x-\Delta _x)-E\right )_+\right ]\right ] \right \}\nonumber \\[5pt] &\qquad +n(t,x,y-\Delta _y)\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [1-\theta -\frac{\eta \mu }{2 E_{\textrm{max}}}\left [ \left (E(t,x+\Delta _x)-E\right )_+ + \left (E(t,x-\Delta _x)-E\right )_+\right ]\right ] \right \}\nonumber \end{align}
\begin{align*} &\qquad +n(t,x+\Delta _x,y)\left \{ (1-\beta )\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\qquad +n(t,x-\Delta _x,y)\left \{ (1-\beta )\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\qquad +n\left \{ (1-\beta )\left [ 1+\tau R\right ] \left [1-\theta -\frac{\eta \mu }{2 E_{\textrm{max}}}\left [\left (E(t,x+\Delta _x)-E\right )_++\left (E(t,x-\Delta _x)-E\right )_+ \right ]\right ]\right \}.\nonumber \end{align*}
Assuming the function
$n$
to be sufficiently regular, we now use the following Taylor expansions
\begin{eqnarray*} n(t,x,y\pm \Delta _y)= n\pm \Delta _y\frac{\partial n}{\partial y}+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}+h.o.t., \quad n(t,x\pm \Delta _x,y)= n\pm \Delta _x\frac{\partial n}{\partial x}+\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}+h.o.t.,\\[5pt] n(t,x+\Delta _x,y\pm \Delta _y)= n+\Delta _x\frac{\partial n}{\partial x} \pm \Delta _y\frac{\partial n}{\partial y}+\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2} \pm \Delta _x\Delta _y \frac{\partial ^2 n}{\partial x\partial y}+h.o.t.,\\[5pt] n(t,x-\Delta _x,y\pm \Delta _y)= n-\Delta _x\frac{\partial n}{\partial x} \pm \Delta _y\frac{\partial n}{\partial y}+\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2} \mp \Delta _x \Delta _y \frac{\partial ^2 n}{\partial x\partial y}+h.o.t., \end{eqnarray*}
which allow us to rewrite (A2) as
\begin{align} &n(t+\tau,x,y)=\\[5pt] &n\left \{\frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}}\left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber\\[5pt] &+\Delta _x\frac{\partial n}{\partial x}\left \{\frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}}\left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\Delta _y\frac{\partial n}{\partial y}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\Delta _x \Delta _y \frac{\partial ^2 n}{\partial x\partial y}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+n\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &-\Delta _x\frac{\partial n}{\partial x}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\Delta _y\frac{\partial n}{\partial y}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &-\Delta _x \Delta _y \frac{\partial ^2 n}{\partial x\partial y}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \end{align}
\begin{align*} &+\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+n\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\Delta _x\frac{\partial n}{\partial x}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &-\Delta _y\frac{\partial n}{\partial y}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &-\Delta _x \Delta _y \frac{\partial ^2 n}{\partial x\partial y}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt]&+n\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &-\Delta _x\frac{\partial n}{\partial x}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &-\Delta _y\frac{\partial n}{\partial y}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber\\[5pt] &+\Delta _x \Delta _y \frac{\partial ^2 n}{\partial x\partial y}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+n\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [1-\theta -\frac{\eta \mu }{2 E_{\textrm{max}}}\left [ \left (E(t,x+\Delta _x)-E\right )_+ + \left (E(t,x-\Delta _x)-E\right )_+\right ]\right ] \right \}\nonumber \\[5pt] &+\Delta _y\frac{\partial n}{\partial y}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [1-\theta -\frac{\eta \mu }{2 E_{\textrm{max}}}\left [ \left (E(t,x+\Delta _x)-E\right )_+ + \left (E(t,x-\Delta _x)-E\right )_+\right ]\right ] \right \}\nonumber \\[5pt] &+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [1-\theta -\frac{\eta \mu }{2 E_{\textrm{max}}}\left [ \left (E(t,x+\Delta _x)-E\right )_+ + \left (E(t,x-\Delta _x)-E\right )_+\right ]\right ] \right \}\nonumber \\[5pt] &+n\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [1-\theta -\frac{\eta \mu }{2 E_{\textrm{max}}}\left [ \left (E(t,x+\Delta _x)-E\right )_+ + \left (E(t,x-\Delta _x)-E\right )_+\right ]\right ] \right \}\nonumber \\[5pt] &-\Delta _y\frac{\partial n}{\partial y}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [1-\theta -\frac{\eta \mu }{2 E_{\textrm{max}}}\left [ \left (E(t,x+\Delta _x)-E\right )_+ + \left (E(t,x-\Delta _x)-E\right )_+\right ]\right ] \right \}\nonumber \end{align*}
\begin{align*} &+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [1-\theta -\frac{\eta \mu }{2 E_{\textrm{max}}}\left [ \left (E(t,x+\Delta _x)-E\right )_+ + \left (E(t,x-\Delta _x)-E\right )_+\right ]\right ] \right \}\nonumber \\[5pt] &+n\left \{ (1-\beta )\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\Delta _x\frac{\partial n}{\partial x}\left \{ (1-\beta )\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\left \{ (1-\beta )\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+n\left \{ (1-\beta )\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &-\Delta _x\frac{\partial n}{\partial x}\left \{ (1-\beta )\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\left \{ (1-\beta )\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+n\left \{ (1-\beta )\left [ 1+\tau R\right ] \left [1-\theta -\frac{\eta \mu }{2 E_{\textrm{max}}}\left [\left (E(t,x+\Delta _x)-E\right )_++\left (E(t,x-\Delta _x)-E\right )_+ \right ]\right ]\right \} + h.o.t.\nonumber \end{align*}
Then rearranging and collecting terms of derivatives of
$n$
, we obtain
\begin{align} &n(t+\tau,x,y)\nonumber\\[5pt] &=n\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\\[5pt] &\quad +n\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad +n\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad +n\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad +n\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [1-\theta -\frac{\eta \mu }{2 E_{\textrm{max}}}\left [ \left (E(t,x+\Delta _x)-E\right )_+ + \left (E(t,x-\Delta _x)-E\right )_+\right ]\right ] \right \}\nonumber \\[5pt] &\quad +n\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [1-\theta -\frac{\eta \mu }{2 E_{\textrm{max}}}\left [ \left (E(t,x+\Delta _x)-E\right )_+ + \left (E(t,x-\Delta _x)-E\right )_+\right ]\right ] \right \}\nonumber \\[5pt] &\quad +n\left \{ (1-\beta )\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad +n\left \{ (1-\beta )\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad +n\left \{ (1-\beta )\left [ 1+\tau R\right ] \left [1-\theta -\frac{\eta \mu }{2 E_{\textrm{max}}}\left [\left (E(t,x+\Delta _x)-E\right )_++\left (E(t,x-\Delta _x)-E\right )_+ \right ]\right ]\right \}\nonumber \end{align}
\begin{align*} &\quad +\Delta _x\frac{\partial n}{\partial x}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber\\[5pt] &\quad -\Delta _x\frac{\partial n}{\partial x}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad +\Delta _x\frac{\partial n}{\partial x}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad -\Delta _x\frac{\partial n}{\partial x}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad +\Delta _x\frac{\partial n}{\partial x}\left \{ (1-\beta )\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad -\Delta _x\frac{\partial n}{\partial x}\left \{ (1-\beta )\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad +\Delta _y\frac{\partial n}{\partial y}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad +\Delta _y\frac{\partial n}{\partial y}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad -\Delta _y\frac{\partial n}{\partial y}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber\\[5pt] &\quad -\Delta _y\frac{\partial n}{\partial y}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad +\Delta _y\frac{\partial n}{\partial y}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [1-\theta -\frac{\eta \mu }{2 E_{\textrm{max}}}\left [ \left (E(t,x+\Delta _x)-E\right )_+ + \left (E(t,x-\Delta _x)-E\right )_+\right ]\right ] \right \}\nonumber \\[5pt] &\quad -\Delta _y\frac{\partial n}{\partial y}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [1-\theta -\frac{\eta \mu }{2 E_{\textrm{max}}}\left [ \left (E(t,x+\Delta _x)-E\right )_+ + \left (E(t,x-\Delta _x)-E\right )_+\right ]\right ] \right \}\nonumber \\[5pt] &\quad +\Delta _x \Delta _y \frac{\partial ^2 n}{\partial x\partial y}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad -\Delta _x \Delta _y \frac{\partial ^2 n}{\partial x\partial y}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad -\Delta _x \Delta _y \frac{\partial ^2 n}{\partial x\partial y}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad +\Delta _x \Delta _y \frac{\partial ^2 n}{\partial x\partial y}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad +\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad +\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad +\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \end{align*}
\begin{align*} &\quad +\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad +\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\left \{ (1-\beta )\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad +\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\left \{ (1-\beta )\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad +\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad +\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad +\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad +\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &\quad +\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [1-\theta -\frac{\eta \mu }{2 E_{\textrm{max}}}\left [ \left (E(t,x+\Delta _x)-E\right )_+ + \left (E(t,x-\Delta _x)-E\right )_+\right ]\right ] \right \}\nonumber \\[5pt] &\quad +\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \frac{\beta }{2}\left [ 1+\tau R\right ]\left [1-\theta -\frac{\eta \mu }{2 E_{\textrm{max}}}\left [ \left (E(t,x+\Delta _x)-E\right )_+ + \left (E(t,x-\Delta _x)-E\right )_+\right ]\right ] \right \} +h.o.t.\nonumber \end{align*}
Further simplifying yields
\begin{align} &n(t+\tau,x,y)=\\[5pt] &n\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}+n\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+n\left \{ \beta \left [ 1+\tau R\right ]\left [1-\theta -\frac{\eta \mu }{2 E_{\textrm{max}}}\left [ \left (E(t,x+\Delta _x)-E\right )_+ + \left (E(t,x-\Delta _x)-E\right )_+\right ]\right ] \right \}\nonumber \\[5pt] &+n\left \{ (1-\beta )\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+n\left \{ (1-\beta )\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+n\left \{ (1-\beta )\left [ 1+\tau R\right ] \left [1-\theta -\frac{\eta \mu }{2 E_{\textrm{max}}}\left [\left (E(t,x+\Delta _x)-E\right )_++\left (E(t,x-\Delta _x)-E\right )_+ \right ]\right ]\right \}\nonumber \\[5pt] &+\Delta _x\frac{\partial n}{\partial x}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber\\[5pt] &-\Delta _x\frac{\partial n}{\partial x}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\Delta _x\frac{\partial n}{\partial x}\left \{ (1-\beta )\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &-\Delta _x\frac{\partial n}{\partial x}\left \{ (1-\beta )\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \end{align}
\begin{align*} &+\Delta _y\frac{\partial n}{\partial y}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber\\[5pt] &-\Delta _y\frac{\partial n}{\partial y}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\left \{ (1-\beta )\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\left \{ (1-\beta )\left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \beta \left [ 1+\tau R\right ]\left [1-\theta -\frac{\eta \mu }{2 E_{\textrm{max}}}\left [ \left (E(t,x+\Delta _x)-E\right )_+ + \left (E(t,x-\Delta _x)-E\right )_+\right ]\right ] \right \} +h.o.t.\nonumber \end{align*}
Cancelling out
$\beta$
terms, where possible, we find
\begin{align} &n(t+\tau,x,y)= n\left \{ \left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\\[5pt] &+n\left \{ \left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+n\left \{ \left [ 1+\tau R\right ]\left [1-\theta -\frac{\eta \mu }{2 E_{\textrm{max}}}\left [ \left (E(t,x+\Delta _x)-E\right )_+ + \left (E(t,x-\Delta _x)-E\right )_+\right ]\right ] \right \}\nonumber \\[5pt] &+\Delta _x\frac{\partial n}{\partial x}\left \{ \left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber\\[5pt] &-\Delta _x\frac{\partial n}{\partial x}\left \{ \left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\Delta _y\frac{\partial n}{\partial y}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber\\[5pt] &-\Delta _y\frac{\partial n}{\partial y}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\left \{ \left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber\\[5pt] &+\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\left \{ \left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \end{align}
\begin{align*} &+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber\\[5pt] &+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\theta }{2}+\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \beta \left [ 1+\tau R\right ]\left [1-\theta -\frac{\eta \mu }{2 E_{\textrm{max}}}\left [ \left (E(t,x+\Delta _x)-E\right )_+ + \left (E(t,x-\Delta _x)-E\right )_+\right ]\right ] \right \} +h.o.t.\nonumber \end{align*}
Cancelling out
$\theta$
terms, where possible, we obtain
\begin{align} &n(t+\tau,x,y)=n\left \{ \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber\\[5pt] &+n\left \{ \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\\[5pt] &+n\left \{ \left [ 1+\tau R\right ]\left [1-\frac{\eta \mu }{2 E_{\textrm{max}}}\left [ \left (E(t,x+\Delta _x)-E\right )_+ + \left (E(t,x-\Delta _x)-E\right )_+\right ]\right ] \right \}\nonumber \\[5pt] &+\Delta _x\frac{\partial n}{\partial x}\left \{ \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &-\Delta _x\frac{\partial n}{\partial x}\left \{ \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\Delta _y\frac{\partial n}{\partial y}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &-\Delta _y\frac{\partial n}{\partial y}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\left \{ \left [ 1+\tau R\right ]\left [\theta +\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\left \{ \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \beta \left [ 1+\tau R\right ]\left [1-\frac{\eta \mu }{2 E_{\textrm{max}}}\left [ \left (E(t,x+\Delta _x)-E\right )_+ + \left (E(t,x-\Delta _x)-E\right )_+\right ]\right ] \right \} +h.o.t.\nonumber \end{align}
We can then rearrange the equations to obtain
\begin{align} &n(t+\tau,x,y)=n\left \{ \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber\\[5pt]&+n\left \{ \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \} \end{align}
\begin{align*} &+n\left [ 1+\tau R\right ]-n\left \{ \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}}\left [ \left (E(t,x+\Delta _x)-E\right )_+\right ]\right ]\right \}\nonumber\\[5pt] &-n\left \{ \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}}\left [ \left (E(t,x-\Delta _x)-E\right )_+\right ]\right ]\right \}\nonumber \\[5pt] &+\Delta _x\frac{\partial n}{\partial x}\left \{ \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber\\[5pt] &-\Delta _x\frac{\partial n}{\partial x}\left \{ \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\Delta _y\frac{\partial n}{\partial y}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber\\[5pt] &-\Delta _y\frac{\partial n}{\partial y}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\theta \left [ 1+\tau R\right ] +\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\left \{ \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\left \{ \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber\\[5pt] &+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\beta \left [ 1+\tau R\right ]-\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}}\left [ \left (E(t,x+\Delta _x)-E\right )_+\right ]\right ]\right \}\nonumber \\[5pt] &-\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}}\left [ \left (E(t,x-\Delta _x)-E\right )_+\right ]\right ]\right \} +h.o.t.\nonumber \end{align*}
Now, using the fact that, for real functions
$f$
, the relation
$(f)_+-(-\!f)_+=f$
holds, we have
\begin{align} &n(t+\tau,x,y)=n\left \{ \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )\right ]\right \}\\[5pt]&+n\left \{ \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )\right ]\right \}+n\left [ 1+\tau R\right ]\nonumber\\[5pt] &+\Delta _x\frac{\partial n}{\partial x}\left \{ \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &-\Delta _x\frac{\partial n}{\partial x}\left \{ \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\Delta _y\frac{\partial n}{\partial y}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \end{align}
\begin{align*} &-\Delta _y\frac{\partial n}{\partial y}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\theta \left [ 1+\tau R\right ] +\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\left \{ \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\left \{ \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )_+\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x+\Delta _x)\right )\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (E-E(t,x-\Delta _x)\right )\right ]\right \}\nonumber \\[5pt] &+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\beta \left [ 1+\tau R\right ]+h.o.t.\nonumber \end{align*}
Next, assuming the function
$E$
to be sufficiently regular, substituting the following Taylor expansion
\begin{equation*} E(t,x\pm \Delta _x)= E\pm \Delta _x\frac {\partial E}{\partial x}+\frac {\Delta _x^2}{2}\frac {\partial ^2 E}{\partial x^2}+h.o.t. \end{equation*}
into (A9) and removing higher-order terms, we obtain
\begin{eqnarray} &&n(t+\tau,x,y)=n\left \{ \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (-\Delta _x^2 \frac{\partial ^2 E}{\partial x^2}\right )\right ]\right \}+n\left [ 1+\tau R\right ]\\[5pt] &&+\Delta _x\frac{\partial n}{\partial x}\left \{ \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (-\Delta _x \frac{\partial E}{\partial x}\right )_+\right ]\right \}-\Delta _x\frac{\partial n}{\partial x}\left \{ \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (\Delta _x \frac{\partial E}{\partial x}\right )_+\right ]\right \}\nonumber \\[5pt] &&+\Delta _y\frac{\partial n}{\partial y}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (-\Delta _x \frac{\partial E}{\partial x}\right )_+\right ]\right \}-\Delta _y\frac{\partial n}{\partial y}\left \{ \beta \left [ 1+\tau R\right ]\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (\Delta _x \frac{\partial E}{\partial x}\right )_+\right ]\right \}\nonumber \\[5pt] &&+\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\theta \left [ 1+\tau R\right ]+\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\beta \left [ 1+\tau R\right ]+h.o.t.\nonumber \end{eqnarray}
Absorbing terms of order
$\mathscr{O}(\tau \Delta _x^2)$
and
$\mathscr{O}(\tau \Delta _y^2)$
into h.o.t. yields
\begin{eqnarray} &&n(t+\tau,x,y)=n\left \{ \left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (-\Delta _x^2 \frac{\partial ^2 E}{\partial x^2}\right )\right ]\right \}+n\left [ 1+\tau R\right ]\\[5pt] &&+\Delta _x\frac{\partial n}{\partial x}\left \{\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (-\Delta _x \frac{\partial E}{\partial x}\right )_+\right ]\right \}-\Delta _x\frac{\partial n}{\partial x}\left \{ \left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (\Delta _x \frac{\partial E}{\partial x}\right )_+\right ]\right \}\nonumber \\[5pt] &&+\Delta _y\frac{\partial n}{\partial y}\left \{ \beta \left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (-\Delta _x \frac{\partial E}{\partial x}\right )_+\right ]\right \}-\Delta _y\frac{\partial n}{\partial y}\left \{ \beta \left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (\Delta _x \frac{\partial E}{\partial x}\right )_+\right ]\right \}\nonumber \\[5pt] &&+\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\theta +\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\beta +h.o.t.\nonumber \end{eqnarray}
Once again, using the relation
$(f)_+-(\!-\!f)_+=f$
for real functions
$f$
, we find
\begin{eqnarray} &&n(t+\tau,x,y)=n\left \{ \left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (-\Delta _x^2 \frac{\partial ^2 E}{\partial x^2}\right )\right ]\right \}+n\left [ 1+\tau R\right ]+\Delta _x\frac{\partial n}{\partial x}\left \{\left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (-\Delta _x \frac{\partial E}{\partial x}\right )\right ]\right \}\\[5pt] &&+\Delta _y\frac{\partial n}{\partial y}\left \{ \beta \left [\frac{\eta \mu }{2 E_{\textrm{max}}} \left (-\Delta _x \frac{\partial E}{\partial x}\right )\right ]\right \}+\frac{\Delta _x^2}{2}\frac{\partial ^2 n}{\partial x^2}\theta +\frac{\Delta _y^2}{2}\frac{\partial ^2 n}{\partial y^2}\beta +h.o.t.\nonumber \end{eqnarray}
Further rearranging gives
\begin{align} n(t+\tau,x,y)&=\frac{\theta \Delta _x^2 }{2}\frac{\partial ^2 n}{\partial x^2}-\frac{\Delta _x^2 \eta \mu }{2 E_{\textrm{max}}}\left [n\frac{\partial ^2 E}{\partial x^2}+\frac{\partial n}{\partial x}\frac{\partial E}{\partial x}\right ] +n+\tau R n +\frac{\Delta _y^2 \beta }{2}\frac{\partial ^2 n}{\partial y^2} \\[5pt] &- \beta \Delta _y \Delta _x\left [\frac{\eta \mu }{2 E_{\textrm{max}}}\frac{\partial E}{\partial x}\right ]\frac{\partial n}{\partial y}+h.o.t.\nonumber \end{align}
Dividing both sides of (A13) by
$\tau$
yields
\begin{align*} \frac{n(t+\tau,x,y)-n}{\tau }&=\frac{\theta \Delta _x^2 }{2\tau }\frac{\partial ^2 n}{\partial x^2}-\frac{\Delta _x^2 \eta \mu }{2 E_{\textrm{max}}\tau }\left [n\frac{\partial ^2 E}{\partial x^2}+\frac{\partial n}{\partial x}\frac{\partial E}{\partial x}\right ] + R n\nonumber \\[5pt] &+\frac{\Delta _y^2 \beta }{2\tau }\frac{\partial ^2 n}{\partial y^2}- \frac{\beta \Delta _y \Delta _x}{\tau }\left [\frac{\eta \mu }{2 E_{\textrm{max}}}\frac{\partial E}{\partial x}\right ]\frac{\partial n}{\partial y}+h.o.t. \end{align*}
Now letting the time-step
$\tau \rightarrow 0$
, the space-step
$\Delta _x\rightarrow 0$
and the phenotype-step
$\Delta _y\rightarrow 0$
in such a way that
\begin{equation*} \frac {\Delta _x^2 \theta }{2 \tau }\rightarrow D \in \mathbb {R}^+_{*}, \quad \frac {\Delta _x^2 \eta }{2 E_{\textrm {max}}\tau }\rightarrow \nu \in \mathbb {R}^+_*,\quad \text {and} \quad \frac {\Delta _y^2 \beta }{2\tau }\rightarrow \lambda \in \mathbb {R}^+_* \end{equation*}
and using the definition for
$\chi (y)$
given by (19), we formally obtain
\begin{eqnarray} &&\frac{\partial n}{\partial t}=D\frac{\partial ^2 n}{\partial x^2}-\chi (y)\left [n\frac{\partial ^2 E}{\partial x^2}+\frac{\partial n}{\partial x}\frac{\partial E}{\partial x}\right ] + R n+\lambda \frac{\partial ^2 n}{\partial y^2}, \end{eqnarray}
which further simplifies to
\begin{eqnarray} &&\frac{\partial n(t,x,y)}{\partial t}=D\frac{\partial ^2 n(t,x,y)}{\partial x^2}-\chi (y)\frac{\partial }{\partial x}\left (n(t,x,y)\frac{\partial E(t,x)}{\partial x}\right )+R(y,\rho ) n(t,x,y)+\lambda \frac{\partial ^2 n(t,x,y)}{\partial y^2},\quad \end{eqnarray}
from which, rearranging terms, we recover the PIDE (2)
$_1$
for
$n(t,x,y)$
.
Remark 2. Under the initial conditions and the assumptions on the model functions and in the asymptotic regime considered here, the cell density
$\rho (t,x)$
behaves like a one-side compactly supported and monotonically decreasing travelling front, while the ECM density
$E(t,x)$
is identically equal to
$0$
on the support of
$\rho (t,x)$
and identically equal to
$E_{\rm max}=1$
outside the support of
$\rho (t,x)$
(i.e. ahead of the cell travelling front). Hence, since for each
$t$
the relation
$\textrm{Supp}\left (n(t,x,y)\right ) \subseteq \textrm{Supp}\left (\rho (t,x)\right )$
holds
$y$
by
$y$
, despite the fact that the ECM density
$E(t,x)$
jumps from
$0$
to
$E_{\rm max}=1$
moving from inside to outside
$\textrm{Supp}\left (\rho (t,x)\right )$
, we expect the formal method that we have employed to derive the PIDE (A15) from the underlying IB model to apply for
$x \in \textrm{Supp}\left (\rho (t,x)\right )$
. This is confirmed by the good agreement between the cell dynamics predicted by numerical simulations of the IB model and numerical solutions of the corresponding continuum model.
A.1.2 Equations for the MDE concentration and the ECM density
Using first the fact that for
$\tau$
,
$\Delta _x$
and
$\Delta _y$
sufficiently small the following relations hold
\begin{eqnarray*} n^{k}_{i,j}\approx n(t,x,y),\quad E^{k}_{i}\approx E(t,x),\quad M^{k}_{i}\approx M(t,x),\quad M^{k}_{i\pm 1}\approx M(t,x\pm \Delta _x),\\[5pt] n^{k+1}_{i,j}\approx n(t+\tau,x,y),\quad E^{k+1}_{i}\approx E(t+\tau,x),\quad M^{k+1}_{i}\approx M(t+\tau,x),\quad p(y_j)\approx p(y), \end{eqnarray*}
and then letting
$\Delta _{\tau } \to 0$
,
$\Delta _{x} \to 0$
and
$\Delta _y \to 0$
in the difference equations (14) and (15), one formally obtains the PDE (2)
$_3$
for
$M(t,x)$
and the infinite-dimensional ODE (2)
$_4$
for
$E(t,x)$
.
Open access
