Derivation and travelling wave analysis of phenotype-structured haptotaxis models of cancer invasion

We formulate haptotaxis models of cancer invasion wherein the infiltrating cancer cells can occupy a spectrum of states in phenotype space, ranging from `fully mesenchymal' to `fully epithelial'. The more mesenchymal cells are those that display stronger haptotaxis responses and have greater capacity to modify the extracellular matrix (ECM) through enhanced secretion of matrix-degrading enzymes (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, which is based on a branching random walk over a lattice representing both physical and phenotype space. We formally derive the corresponding continuum model, which takes the form of a coupled system comprising a partial integro-differential equation for the local cell population density function, a partial differential equation for the MDE concentration and an infinite-dimensional ordinary differential equation for the ECM density. Despite the intricacy of the model, we show, through formal asymptotic techniques, that for certain parameter regimes it is possible to carry out a detailed travelling wave analysis and obtain invading fronts 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.


Introduction 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 [27].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 [68].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 [68].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 [11].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 [38,60,68]; the directed movement that results from migration up a matrix adhesivity gradient is referred to as haptotaxis [14].Cells may also realign fibres, leading to an oriented matrix that directs invasion along certain paths [58].Leader cells may also start to secrete (or increase the secretion of) matrix-degrading enzymes (MDEs) [17,36,60,68,73].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 [17,32,38,68].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 [38].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 [38].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' [11,17,39,48,59,68].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 [52,68,72,74].
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 [9,68]; however, if leader cells are removed then the invasion of the tumour stops [38].Furthermore, through inhibiting MDE secretion by leader cells then further invasion can be slowed or prevented [13].

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 [2,26,61].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.[5,56].Focussing on a 1D spatial scenario, as a prototypical example of these haptotaxis models of cancer invasion we can consider the following system where the functions ρ ≡ ρ(t, x), M ≡ M(t, x), and E ≡ E(t, x) model, respectively, the cancer cell density, the MDE concentration, and the ECM density at time t ∈ R + and position x ∈ R. In the system (1.1), the parameters D and µ 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(ρ) 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 ρ 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 γ is the rate of MDE production by cancer cells, and the parameter κ M is the rate of natural decay of MDEs.Finally, the parameter κ E is the rate at which the ECM is degraded by MDEs upon contact.Haptotaxis models of cancer invasion of the form of (1.1) have been studied analytically through proof of local existence, global existence, boundedness, and uniqueness of solutions [66], investigations into blow-up of solutions [62], and analysis of 2D radially symmetric solutions [10].Numerical simulations have also been used to investigate models of the form of (1.1) [28].Further models extend systems of this form, for example to include more detailed mechanisms of ECM remodelling and enzyme activities [3,15,47], or through the inclusion of non-local terms to incorporate the effects of cell-cell adhesion [29,51].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 [18].While model (1.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 [37,57,63].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 [65].
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.1) to include a continuous structuring variable representing intercellular variability in certain phenotypic characteristics [54].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 [40], and how trade-offs between mobility and proliferation may impact on density-or pressure-driven growth processes [41,46]; directly relevant to cancer, structured phenotype models of this type have also been developed to explore avascular tumour growth [23] and the evolutionary dynamics underpinning the emergence of intra-tumour phenotypic heterogeneity [24,42,69].In recent work by Guilberteau et al. [31], 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 populationlevel 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) [70,71].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 [5].Here, a model of the type of (1.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 [6], 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 [26].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 [5].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.[50,53,64].

Outline
In this paper, we consider the following generalisation of the haptotaxis model of cancer invasion (1.1), where the continuous structuring variable y ∈ [0, Y] ⊂ R + , with Y > 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: Compared to model (1.1), here the PDE for the cell density, ρ(t, x), is replaced by the partial integro-differential equation (PIDE) (1.2) 1 for the local cell population density function, n ≡ n(t, x, y), which is linked to the cell density through the relation (1.2) 2 .Moreover, the functions χ(y), p(y), and R(y, ρ) 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 (1.2) 1 models the effect of phenotypic changes, which occur at rate λ.
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 [34], 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 [12,16,45,46], we formally derive the model (1.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 [40,41], we carry out travelling wave analysis of an appropriately rescaled version of the model (1.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).

The individual-based model
In this section, integrating the modelling approaches that we developed in [12,46], 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 R. Furthermore, we describe the phenotypic state of each cell by means of a structuring variable y ∈ [0, Y] ⊂ 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 [38].Hence, cells in phenotypic states characterised by values of y closer to 0 display a more epitheliallike 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 ∈ R + and the space variable x ∈ R, respectively, as t k = kτ and and ∆ x ∈ R + * , where R + * denotes the set of positive real numbers.Moreover, we discretise the phenotype variable via y j = j∆ y ∈ [0, Y] with j ∈ N 0 and ∆ y ∈ R + * .Here, τ, ∆ x , and ∆ 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∈Z ×{y j } j∈N 0 and we introduce the dependent variable N k i, j ∈ 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 The concentration of MDEs and the density of ECM at position x i at time t k are denoted by 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.

Modelling the dynamics of cancer cells
As summarised in Figures 2a.-d., between time-steps k and k + 1, each cell in phenotypic state y j ∈ (0, Y) at position x i ∈ 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.

Random cell movement
We model undirected cell movement as a random walk along the spatial dimension, with movement probability 0 < θ ≤ 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 k L i, j or P k R i, j , respectively.As we consider this random movement to be undirected and not affected by the cell phenotype, we define Note that cells will not undergo random movement with probability

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 in phenotypic state y j ∈ (0, Y) at spatial position x i ∈ R may: a. move via random motion to either of the positions x i−1 and x i+1 with probabilities P k L i, j and P k R i, j defined via (2.2) or do not undergo random movement with probability P k S i, j = 1 − P k L i, j + P k R i, j ; b. move via haptotactic motion to either of the positions x i−1 and x i+1 with probabilities P k HL i, j and P k HR i, j defined via (2.4) or do not undergo haptotactic movement with probability P k HS i, j = 1 − P k HL i, j + P k HR i, j ; c. undergo a phenotypic change and thus enter either of the phenotype states y j−1 and y j+1 with probabilities P k D i, j and P k U i, j defined via (2.5) or remain in the same phenotypic state with probability P k N i, j = 1 − P k D i, j + P k U i, j ; d. die and divide with probabilities P k A i, j and P k B i, j defined via (2.6) or remain quiescent with probability P k Q i, j = 1 − P k A i, j + P k B i, j .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 κ M .The ECM density will change over time through: h.degradation by MDEs at rate κ E .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 µ(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 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 k HL i, j or move to the position x i+1 (i.e.move right) with probability P k HR i, j , where we define Here, E max ∈ R + * is the maximum value of the ECM density before cell invasion starts (see also Section 2.2.2).Moreover, the parameter η ∈ R + * is a scaling factor which we consider small enough to ensure η µ(y j ) ≤ 1.Hence, the quantities defined via (2.4) satisfy 0 < P k HL i, j + P k HR i, j ≤ 1 for all values of i, j, and k.Note that cells will not undergo haptotactic movement with probability P k HS i, j := 1 − P k HL i, j + P k HR i, j .

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 < β ≤ 1, or remains in its current phenotypic state with probability 1 − β.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 k D i, j or enters state y j+1 with probability P k U i, j , where we define Therefore cells will not undergo phenotypic changes with probability 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].

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 timestep 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 ρ k i .In particular, to define the probabilities of cell division and death, we introduce the function R(y j , ρ k i ), 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 k A i, j , divide with probability P k B i, j , or remain quiescent with probability P k Q i, j := 1 − P k A i, j + P k B i, j , where By considering the time-step τ sufficiently small, we ensure P k A i, j + P k B i, j ≤ 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 ρ max ∈ R + * .Therefore, we make the following assumptions R(Y, 0) = 0, R(0, ρ max ) = 0, ∂R(y, ρ) ∂ρ < 0 and ∂R(y, ρ) ∂y < 0 for (y, ρ) ∈ (0, Y)×R + .(2.7) In particular, we focus on a similar case to that considered in [46], that is, we assume with α ∈ R + * .Note that, under assumptions (2.7), the definitions given by (2.6) ensure that if ρ k i ≥ ρ 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

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 the following subsections, which are summarised by the schematics in Figures 2e.-h.and are coupled with the individual-based model for the dynamics of cancer cells that is presented in Section 2.1.

Dynamics of the MDE concentration
We let D M ∈ R + * be the diffusivity of the MDEs and we denote by κ M ∈ 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 In this framework, the principle of mass balance gives us the following difference equation for the concentration of MDEs where L is the finite-difference Laplacian on the lattice {x i } i∈Z , that is,

Dynamics of the ECM density
We denote by κ E ∈ 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 Recalling that, as mentioned earlier, E max ∈ R + * is the maximum value of the ECM density before cell invasion starts, we complement the difference equation (2.13) with an initial condition such that max so that 3 The corresponding continuum model Through an extension of the limiting procedure that we previously employed in [12,16,45,46], letting the time-step τ → 0, the space-step ∆ x → 0, and the phenotype-step ∆ y → 0 in such a way that and introducing the definition χ(y) := ν µ(y),

Formal asymptotic analysis
In this section, building on the formal asymptotic method that we developed in [40,41], which relies on the Hamilton-Jacobi approach developed in [8,20,44,55], we carry out travelling wave analysis of an appropriately rescaled version of the model (1.2).

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 [4,5,33,49,67].To this end, we introduce a small parameter ε ∈ R + * and choose the parameter scaling Moreover, in order to explore the long-time behaviour of the system, we use the time scaling t → t ε in the model (1.2).In so doing, recalling the definition given by (3.2), we obtain the following rescaled system for the local cell population density function, n ε (t, x, y), the MDE concentration, M ε (t, x), and the ECM density, E ε (t, x): 4.2 Formal limit for ε → 0 We make the real phase WKB ansatz [7,22,25] which gives Substituting the above expressions into the PIDE (4.2) 1 for n ε and rearranging terms gives the following Hamilton-Jacobi equation for u ε ≡ u ε (t, x, y) x) be the leading-order term of the asymptotic expansion for ρ ε (t, x) as ε → 0. Considering x ∈ R such that ρ(t, x) > 0 (i.e.x ∈ Supp(ρ)) and letting ε → 0 in the above PDE we formally obtain the following equation for the leading-order term u ≡ u(t, x, y) of the asymptotic expansion for u ε (t, x, y) where E ≡ E(t, x) is the leading-order term of the asymptotic expansion for E ε (t, x).
Expressions of M(t, x) and E(t, x) When n ε is in the form (4.3), if u ε is a strictly concave function of y and u is also a strictly concave function of y that satisfies the constraint (4.5) then the following asymptotic result formally holds where δ ȳ(t,x) (y) is the Dirac delta centred at y = ȳ(t, x).In this case, focussing on a biological scenario wherein the ECM density is at the maximum level E max before cell invasion starts at t = 0, letting ε → 0 in the PDE (4.2) 3 for M ε and in the infinite-dimensional ODE (4.2) 4 for E ε we formally obtain the following expressions of the leading-order terms of the asymptotic expansions for M ε (t, x) and E ε (t, x) where 1 (•) denotes the indicator function of the set (•).
Remark 1.Note that the behaviour of E(t, x) depicted by (4.8) shares similarities with the behaviour of the nutrient concentration in the model analysed in [35].
Transport equation for ȳ When R(y, ρ) is defined via (2.8), differentiating (4.4) with respect to y, evaluating the resulting equation at y = ȳ(t, x), and using (4.5) and (4.6) yields Moreover, differentiating (4.6) with respect to t and x we find, respectively, and x ∈ Supp(ρ).Substituting the above expressions of ∂ 2 yt u(t, x, ȳ) and ∂ 2 yx u(t, x, ȳ) into (4.9), and using the fact that if u is a strictly concave function of y whose unique maximum point is ȳ(t, x) then ∂ 2 yy u(t, x, ȳ) < 0, gives the following transport equation for ȳ(t, x) , x ∈ Supp(ρ).(4.10)

Numerical simulations
In this section, we report on numerical solutions of the rescaled continuum model (4.2) 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.

Set-up of numerical simulations
We start by describing the set-up used to carry out numerical simulations.

Model functions and parameters
To allow the individual-based model to represent the same scenario as the rescaled continuum model (4.2), we use the same time scaling t k → t k ε = k τ ε and reformulate the governing rules for the cell dynamics detailed in Section 2 in terms of To ensure that conditions (3.1) and (4.1) are simultaneously satisfied, we additionally set 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

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 = ȳ0 at every position x ∈ [0, X].Specifically, we implement the following initial cell distribution for the IB model where int(•) is the integer part of (•) and C is a normalisation constant such that We choose ȳ0 = 0.2 and A 0 = 100.The initial cell density ρ 0 i is then calculated from (5.4) according to the definition given by (2.1), and we set ρ max = max i ρ 0 i .Moreover, we assume that there are initially no MDEs and the density of ECM is uniform, that is, Finally, we consider different values of ε, that is, ε ∈ 10 −2 , 5 × 10 −3 , 10 −3 , in order to verify whether, for ε small enough, there is a good agreement between the results of numerical simulations and the results of formal asymptotic analysis for ε → 0 presented in Section 4.

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 (4.2) was also implemented in Matlab.
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 (2.2); (ii) haptotactic cell movement, according to the probabilities defined via (2.4); (iii) phenotypic changes, according to the probabilities defined via (2.5); (iv) cell division and death, according to the probabilities defined via (2.6).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 (2.12) and (2.13), respectively.

Numerical methods for the continuum model
To solve numerically the rescaled system (4.2) posed on (0, T ] × (0, X) × (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 [12].

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 (4.2) for three distinct values of the scaling parameter ε are summarised by the plots in Figures 3-5, which correspond to ε = 10 −2 , ε = 5 × 10 −3 , and ε = 10 −3 , respectively.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 k ε i, j , the cell density, ρ k ε i , the MDE concentration, M k ε i , and the ECM density, E k ε i , 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 ε, 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 ε → 0, invading fronts with spatial structuring of cell phenotypes.In particular, the plots in Figure 5 demonstrate that, when ε is sufficiently small: (i) The local cell population density function n ε (t, x, y) becomes concentrated as a sharp Gaussian with maximum at a point ȳε (t, x) for all x where ρ ε (t, x) > 0.
(ii) The maximum point ȳε (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 ≈ 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 ρ ε (t, x) behaves like a one-sided compactly supported and monotonically decreasing travelling front that connects ρ max to 0.
(iv) The values of ρ ε , M ε , and E ε are consistent with the values obtained by substituting ȳ(t, x) = ȳε (t, x) into the formulas given by (4.7) and (4.8).

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 (1.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 [68].
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 [30].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 [57,63], although more complex movement models have also been considered, e.g.[19].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 [65], 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 [23].
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 [40] and density-or pressure-driven wave invasion in [41,46]; 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 [6] and tumour growth in the presence of cells with different mobilities [21,43].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.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 [68], 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 [38]: 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 [1].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. .
Using the fact that for τ, ∆ x , and ∆ y sufficiently small the following relations hold Eq (A.1) can be rewritten as Assuming the function n to be sufficiently regular, we now use the following Taylor expansions which allow us to rewrite (A.2) as Then rearranging and collecting terms of derivatives of n, we obtain Further simplifying yields Cancelling out β terms, where possible, we find Remark 2. Under the initial conditions and the assumptions on the model functions and in the asymptotic regime considered here, the cell density ρ(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 ρ(t, x) and identically equal to 1 outside the support of ρ(t, x) (i.e.ahead of the cell travelling front).Hence, since for each t the relation Supp (n(t, x, y)) ⊆ Supp (ρ(t, x)) holds y by y, despite the fact that the ECM density E(t, x) jumps from 0 to 1 moving from inside to outside Supp (ρ(t, x)), we expect the formal method that we have employed to derive the PIDE (A.15) from the underlying IB model to apply for x ∈ Supp (ρ(t, x)).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 τ, ∆ x , and ∆ y sufficiently small the following relations hold

Figure 1 :
Figure 1: Schematic overview of the trade-offs between haptotactic and MDE production ability and proliferative potential incorporated into the individual-based model.

Figure 2 :
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 ∈ (0, Y) at spatial position x i ∈ R may: a. move via random motion to either of the positions x i−1 and x i+1 with probabilities P kL i, j and P k R i, j defined via (2.2) or do not undergo random movement with probability P k S i, j = 1 − P k L i, j + P k R i, j ; b. move via haptotactic motion to either of the positions x i−1 and x i+1 with probabilities P k HL i, j and P k HR i, j defined via (2.4) or do not undergo haptotactic movement with probability P k HS i, j = 1 − P k HL i, j + P k HR i, j ; c. undergo a phenotypic change and thus enter either of the phenotype states y j−1 and y j+1 with probabilities P k D i, j and P k U i, j defined via (2.5) or remain in the same phenotypic state with probability P k N i, j = 1 − P k D i, j + P k U i, j ; d. die and divide with probabilities P k A i, j and P k B i, j defined via (2.6) or remain quiescent with probability P k Q i, j = 1 − P k A i, j + P k B i, j .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 κ M .The ECM density will change over time through: h.degradation by MDEs at rate κ E .

Figure 3 :
Figure 3: Numerical simulation results of the individual-based model (top row) and numerical solutions of the corresponding continuum model (4.2) (bottom row) in the case where ε = 10 −2 .The plots display, from left to right, the cell population density, n ε , the cell density, ρ ε , the MMP concentration, M ε , and the ECM density, E ε , 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 5 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 ρ ε , M ε , and E ε (solid red lines) are plotted along with the corresponding values obtained by substituting ȳ(t, x) = ȳε (t, x) into the formulas given by (4.7) and (4.8) (black dotted lines), with ȳε (t, x) being the maximum point of the numerical solution n ε (t, x, y) to the PIDE (4.2) 1 at position x at time t.

Figure 4 :
Figure 4: Numerical simulation results of the individual-based model (top row) and numerical solutions of the corresponding continuum model (4.2) (bottom row) in the case where ε = 5 × 10 −3 .The plots display, from left to right, the cell population density, n ε , the cell density, ρ ε , the MMP concentration, M ε , and the ECM density, E ε , 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 5 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 ρ ε , M ε , and E ε (solid red lines) are plotted along with the corresponding values obtained by substituting ȳ(t, x) = ȳε (t, x) into the formulas given by (4.7) and (4.8) (black dotted lines), with ȳε (t, x) being the maximum point of the numerical solution n ε (t, x, y) to the PIDE (4.2) 1 at position x at time t.

Figure 5 :
Figure 5: Numerical simulation results of the individual-based model (top row) and numerical solutions of the corresponding continuum model (4.2) (bottom row) in the case where ε = 10 −3 .The plots display, from left to right, the cell population density, n ε , the cell density, ρ ε , the MMP concentration, M ε , and the ECM density, E ε , 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 5 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 ρ ε , M ε , and E ε (solid red lines) are plotted along with the corresponding values obtained by substituting ȳ(t, x) = ȳε (t, x) into the formulas given by (4.7) and (4.8) (black dotted lines), with ȳε (t, x) being the maximum point of the numerical solution n ε (t, x, y) to the PIDE (4.2) 1 at position x at time t.

Figure 6 :
Figure 6: Preliminary 2D results from the individual-based model, obtained by averaging over 15 simulations, in the case where ε = 10 −2 .Top and middle rows.The plots display, from left to right, the maximum point of the cell population density, ȳε = arg max y∈[0,Y] n ε , the cell density, ρ ε , the MMP concentration, M ε , and the ECM density, E ε , 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, ȳε , ρ ε , M ε , and E ε across the radial transect at the end of simulations (i.e. at t = 5).Here x ≡ (x 1 , x 2 ) ∈ [0, 10]×[0, 10] with ∆ x 1 = ∆ x 2 = 0.1, y ∈ [0, 1] with ∆ y = 0.02, and the initial cell distribution is the 2D analogue of (5.4) with A 0 = 1, while all the other parameters and functions are kept the same as in the 1D simulations of Figures 3-5.