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

TL;DR

Despite the intricacy of the model, it is shown 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.

Abstract

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.

Figures (6)

• 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: 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 \ref{['eq:define_randmov']} or do not undergo random movement with probability ${P}_{S_{i,j}}^k=1-\left({P}_{L_{i,j}}^k+{P}_{R_{i,j}}^k\right)$; 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 \ref{['eq:define_haptomov']} or do not undergo haptotactic movement with probability ${P}_{HS_{i,j}}^k=1-\left({P}_{HL_{i,j}}^k+{P}_{HR_{i,j}}^k\right)$; 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 \ref{['eq:define_phenochange']} or remain in the same phenotypic state with probability ${P}_{N_{i,j}}^k=1-\left({P}_{D_{i,j}}^k+{P}_{U_{i,j}}^k\right)$; d. die and divide with probabilities ${P}_{A_{i,j}}^k$ and ${P}_{B_{i,j}}^k$ defined via \ref{['eq:define_proliferation']} or remain quiescent with probability ${P}_{Q_{i,j}}^k=1-\left({P}_{A_{i,j}}^k+{P}_{B_{i,j}}^k\right)$. 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$.
• Figure 3: Numerical simulation results of the individual-based model (top row) and numerical solutions of the corresponding continuum model \ref{['eq:PDE_scaled']} (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 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 $\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 \ref{['eq:Riszero']} and \ref{['eq:lotME']} (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 \ref{['eq:PDE_scaled']}$_1$ at position $x$ at time $t$.
• Figure 4: Numerical simulation results of the individual-based model (top row) and numerical solutions of the corresponding continuum model \ref{['eq:PDE_scaled']} (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 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 $\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 \ref{['eq:Riszero']} and \ref{['eq:lotME']} (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 \ref{['eq:PDE_scaled']}$_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 \ref{['eq:PDE_scaled']} (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 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 $\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 \ref{['eq:Riszero']} and \ref{['eq:lotME']} (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 \ref{['eq:PDE_scaled']}$_1$ at position $x$ at time $t$.

Theorems & Definitions (2)

• Remark 1
• Remark 2