Phenotype structuring in collective cell migration:a tutorial of mathematical models and methods

TL;DR

This work frames phenotype-structured populations within collective cell migration by extending classical reaction-advection-diffusion models to phenotype-structured PDEs (PS-PDEs) that track density $n(t,\bm x,\bm y)$ across physical and phenotypic spaces. It provides a tutorial-style guide: deriving PS-PDEs from agent-based models, analyzing travelling waves and concentration phenomena, and implementing simulations via the method of lines, while illustrating with diffusion-based (Fisher–KPP), pressure-based, and taxis-based movement forms. Key contributions include formal derivations linking micro-scale rules to PS-PDEs, analytical insights into how phenotypic structure sorts across travelling waves and drives concentration, and practical numerical schemes for high-dimensional, non-local PDEs. The framework enables mechanistic exploration of how phenotypic heterogeneity shapes invasion, adaptation, and therapy responses, with broad implications for cancer invasion, bacterial chemotaxis, and tissue dynamics; it also highlights methodological challenges in connecting PS-PDEs with data and in developing efficient numerical tools.

Abstract

Populations are heterogeneous, deviating in numerous ways. Phenotypic diversity refers to the range of traits or characteristics across a population, where for cells this could be the levels of signalling, movement and growth activity, etc. Clearly, the phenotypic distribution -- and how this changes over time and space -- could be a major determinant of population-level dynamics. For instance, across a cancerous population, variations in movement, growth, and ability to evade death may determine its growth trajectory and response to therapy. In this review, we discuss how classical partial differential equation (PDE) approaches for modelling cellular systems and collective cell migration can be extended to include phenotypic structuring. The resulting non-local models -- which we refer to as phenotype-structured partial integro-differential equations (PS-PIDEs) -- form a sophisticated class of models with rich dynamics. We set the scene through a brief history of structured population modelling, and then review the extension of several classic movement models -- including the Fisher-KPP and Keller-Segel equations -- into a PS-PIDE form. We proceed with a tutorial-style section on derivation, analysis, and simulation techniques. First, we show a method to formally derive these models from underlying agent-based models. Second, we recount travelling waves in PDE models of spatial spread dynamics and concentration phenomena in non-local PDE models of evolutionary dynamics, and combine the two to deduce phenotypic structuring across travelling waves in PS-PIDE models. Third, we discuss numerical methods to simulate PS-PIDEs, illustrating with a simple scheme based on the method of lines and noting the finer points of consideration. We conclude with a discussion of future modelling and mathematical challenges.

Figures (7)

• Figure 1: Examples of phenotype heterogeneity and its experimental exposition. (a) Heterogeneity has been intensively studied for its role during cancer invasion, one convenient concept being a division into 'followers' and 'leaders'. Leader cells have proteolytic and migratory capabilities that allow them to modify and infiltrate the surrounding extracellular matrix; follower cells expand into the space behind the migrating leaders. For further information, e.g. see vilchez2021decoding. (b) An example of how experiments can reveal phenotypic structuring in a population. A T-maze sorts E. coli cells according to their chemotactic sensitivity: those with stronger responses penetrate deeper into the maze. Figure adapted from the study of salek2019bacterial. (c) A range of proteomic methods can provide quantitative data on phenotypic states, including: (i) population-wide average levels of protein expression (e.g. western blot); (ii) protein expression distributions across the population, also from bulk measurements (e.g. flow cytometry); (iii) spatial distributions of protein expression levels, at single-cell and tissue levels (e.g. immunofluorescence). Adapted from celia2018hysteresis and licensed under CC BY 4.0 (https://creativecommons.org/licenses/by/4.0/), see the original manuscript for precise information on the depicted data. (d) Potential ways of classifying phenotypic variation: (Top) Simple binary variation, where a population is divided into two principal types (e.g. the 'follower' or 'leader' types in (a)); (Middle) Phenotypes spanning a 1D continuum (e.g. chemotactic sensitivities revealed in T-maze experiments); (Bottom) Phenotypes spanning a higher dimension continuum, e.g. migratory capacity, proliferative capacity, and resistance (i.e. ability to evade death).
• Figure 2: Illustrative output of the McKendrick-von Foerster equation. Top panels show the evolving rescaled age distribution, $n(t,a)/\rho(t)$; bottom panels plot the age distribution $n(t,a)$ at select times. Here we assume in (\ref{['vonfoerster']}) that: (i) the population is initially concentrated towards newborns -- i.e. $n^0(a):=\rho^0 e^{-a}$, where $\rho^0$ is proportional to the initial size of the population (for convenience, we set $\rho_0 = 1$); (ii) the death rate increases linearly with age (i.e. $\mu(t,a,n,\rho) \equiv \mu(a) := 0.001a$); and (iii) births result from a Gaussian-like distribution centred about age $a=25$ (i.e. $\beta(t,a,n,\rho) \equiv \beta(a):=(4\exp^{-(a-25)^2/(2\sigma^2)})/(\sqrt{2\pi \sigma^2})$), for (a) dispersed births ($\sigma =4$), and (b) concentrated births ($\sigma =0.25$). In the initial phase, we see a generation that ages with time and new birth "spikes" yielding a new generation as the previous generation reaches $a=25$. However, over longer times, the rescaled age distribution reaches a steady-state profile (though this takes longer when births are more concentrated). Note that a steady-state rescaled age distribution does not imply the population size itself is stabilising -- e.g. for these simulations $\rho(t)$ increases exponentially in time.
• Figure 3: Illustration of travelling wave like solutions from the Fisher-KPP model. Illustrative outputs of the Fisher-KPP model \ref{['FKPP']}-\ref{['def:RFKPP']}. In particular, evolution towards a travelling wave form is shown for the rescaled PDE \ref{['FisherKPP:rescaled']}, complemented with an initial condition of form \ref{['FisherKPP:rescaledIC']}, with (a) $\varepsilon = 1$, (b) $\varepsilon = 0.1$, (c) $\varepsilon=0.01$. In each case the density $\rho(t,x)$ is plotted at times $t=5$ (blue lines), $t=10$ (green lines), $t=15$ (red lines), and $t=20$ (red lines). We see an evolution of the solutions towards a travelling wave profile with wave speed $c=2$ (numerically calculated wave speeds $v$ reported for each simulation at times $t=5$ and $t=20$). As $\varepsilon \rightarrow 0$, this occurs on a fast timescale and the profile becomes increasingly step-like, from $\rho \equiv 1$ across the wave to $\rho \equiv 0$ ahead of the wave, in agreement with analytical predictions. While travelling waves are usually considered on an unbounded domain, the numerical method requires a bounded domain. For these simulations we have set ${\bm x} \equiv x \in[0,50]$, subject to zero-flux boundary conditions.
• Figure 4: Illustration of travelling wave like solutions from a PS-PDE model of diffusion-based movement. Illustrative output of the PS-PDE model of diffusion-based movement \ref{['structuredFKPP']}. Here, we consider one dimension in each of physical and phenotype space, such that ${\bm x} \equiv x \in[0,L]$ with $L \in \mathbb{R}^+_*$ and ${\bm y} \equiv y \in[0,1]$. Zero-flux conditions are imposed at the boundaries and the population is initially set such that $n(0,x,y)=0.1 e^{-100x}$. We take $D(y):=y$ and define $R(y,\rho) := r(y) - \kappa \rho$, with $r(y) := 10 (1-y)$ and $\kappa=10$. These definitions describe a trade-off in which the most motile population members ($y \approx 1$) have the lowest proliferation rate, while the least motile members ($y \approx 0$) have the highest proliferation rate. Simulations are for (a) $\bar{D}=10^{-1}$, (b) $\bar{D}=10^{-3}$, and (c) $\bar{D}=10^{-5}$. Note that $L=100$ in (a-b) and $L=200$ in (c), due to the longer time it takes to evolve to the travelling wave profile in (c). Top row displays the phenotype density function $n(t,x,y)$ with the different colour-themes corresponding to different times (colour bars in (c)). Central row displays the density $\rho(t,x)$ at the same times, showing the evolution to travelling wave solutions. Bottom row displays the density $\rho$ at (a-b) $t=25$ and (c) $t=50$, where the colour code indicates the locally prevailing phenotype, $\bar{\bm y}(t,{\bm x}) \equiv \bar{y}(t,x)$, across the wave for (a-b) $t=25$ and (c) $t=50$. Broadly speaking, lowering the diffusivity in phenotype space (i.e. reducing the value of the rate of phenotypic changes $\bar{D}$) leads to an increasingly structured profile, where more mobile members dominate the front of the wave and more proliferative members dominate the rear.
• Figure 5: Illustration of Gaussian-like solutions and concentration phenomena from a non-local PDE model of evolutionary dynamics. Illustrative output of the non-local PDE model of evolutionary dynamics \ref{['nonlocalFKPP']}-\ref{['def:pEx1']}, with ${\bm y} \equiv y \in [-5,5]$, subject to zero-flux boundary conditions. In particular, Gaussian-like solutions and concentration phenomena are shown for the rescaled non-local PDE \ref{['nonlocal-FisherKPP:rescaled']}, under definitions \ref{['def:REx1']},\ref{['def:pEx1']} (with $\kappa=0.1$, $\gamma=10$, and $\varphi=1$), subject to initial condition \ref{['eq:ICNLFKPPrev']} (with $\rho^0=60$ and $\bar{y}^0=-1$). Simulations are for (a) $\varepsilon=1$, (b) $\varepsilon=0.1$, and (c) $\varepsilon=0.01$. Top panels display the dynamics of the rescaled phenotype density $n(t,y)/\rho(t)$; bottom panels display the phenotype density $n(t,y)$ (red solid line) at select times. We compare the numerical solution with the analytically predicted one (black dashed lines) -- i.e. the Gaussian-like solution \ref{['eq:SOLNLFKPPreveps']} with mean $\bar{y}_\varepsilon(t)$ and weight $\rho_\varepsilon(t)$ obtained by solving numerically the Cauchy problem \ref{['eq:NLFKPPodesreveps']}. We see the solution to be of a Gaussian form, which becomes concentrated as a weighted increasingly sharp Gaussian as $\varepsilon$ gets smaller, corroborating the emergence of concentration phenomena of type \ref{['eq:concphen']}. We see also that the mean or prevailing phenotype $\bar{y}(t)$ converges to $y=\varphi$ (black dashed lines, top panels) as $t \to \infty$, corresponding to selection of the fittest trait.