Coarse-graining and hybrid methods for efficient simulation of stochastic multi-scale models of tumour growth

Abstract

The development of hybrid methodologies is of current interest in both multi-scale modelling and stochastic reaction-diffusion systems regarding their applications to biology. We formulate a hybrid method for stochastic multi-scale models of cells populations that extends the remit of existing hybrid methods for reaction-diffusion systems. Such method is developed for a stochastic multi-scale model of tumour growth, i.e. population-dynamical models which account for the effects of intrinsic noise affecting both the number of cells and the intracellular dynamics. In order to formulate this method, we develop a coarse-grained approximation for both the full stochastic model and its mean-field limit. Such approximation involves averaging out the age-structure (which accounts for the multi-scale nature of the model) by assuming that the age distribution of the population settles onto equilibrium very fast. We then couple the coarse-grained mean-field model to the full stochastic multi-scale model. By doing so, within the mean-field region, we are neglecting noise in both cell numbers (population) and their birth rates (structure). This implies that, in addition to the issues that arise in stochastic-reaction diffusion systems, we need to account for the age-structure of the population when attempting to couple both descriptions. We exploit our coarse-graining model so that, within the mean-field region, the age-distribution is in equilibrium and we know its explicit form. This allows us to couple both domains consistently, as upon transference of cells from the mean-field to the stochastic region, we sample the equilibrium age distribution. Furthermore, our method allows us to investigate the effects of intracellular noise, i.e. fluctuations of the birth rate, on collective properties such as travelling wave velocity. We show that the combination of population and birth-rate noise gives rise to large fluctuations of the birth rate in the region at the leading edge of front, which cannot be accounted for by the coarse-grained model. Such fluctuations have non-trivial effects on the wave velocity. Beyond the development of a new hybrid method, we thus conclude that birth-rate fluctuations are central to a quantitatively accurate description of invasive phenomena such as tumour growth.

Keywords: Age-structured model; Hybrid methods; Multi-scale modelling; Reaction–diffusion systems; Tumour growth.

Fig. 1

Schematic representation of the different elements that compose our multi-scale model. We show the different levels of biological organisation as well as associated characteristic time scales , associated to each of these layers: resource scale, i.e. oxygen which is supplied at a constant rate and consumed by the cell population, cellular scale, i.e. oxygen-regulated cell cycle progression which determines the age-dependent birth rate into the cellular layer, and, finally, the cellular scale, which is associated to the stochastic population dynamics.

Fig. 2

Representation of the setting of our model. Diffusible substances (e.g. oxygen), c(t,x), is modelled as a continuous field described by a reaction–diffusion PDE, it is represented by green solid line. The birth-and-death dynamics with diffusion of the cell population is modelled by means of a RDME on a lattice $\mathcal{L}$. Each vertex of the lattice, $x_i\in \mathcal{L}$, is associated to a compartment or voxel within which the population is assumed to be well-mixed and its stochastic dynamics ruled by a local law of mass action. L is the total length of the system and h is the lattice spacing, so that $L=N_{\mathcal{L}}h$ where $N_{\mathcal{L}}=\text{card}(\mathcal{L})$. Here N(t,x_i) depict the number of cells in compartment i and it is calculated as $N(t,x_i)={\int }_0^{\infty }n(t,a,x_i)da$, with n(t,a,x_i) being the number of cells of age a at time t and compartment x_i. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3

This plot shows a schematic representation of the characteristic curves, a = t + a₀, corresponding to our age-structured stochastic dynamics and the emergence of new genealogies (red line) when a birth occurs (indicated by the red dashed line) within a previously existing one. Genealogies terminate when the corresponding population becomes extinct (indicated by the red point). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4

Plot showing the time evolution of the average front position for the three models: the coarse-grained mean-field model (green line), hybrid model (blue line), and stochastic model (red line). Results shown for the hybrid and stochastic models correspond to an average over 100 and 40 realisations respectively. Θ = 2000. For other parameter values see Section 5. Position calculation: average of the positions x where the population is greater than 0 and smaller than (K-100). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5

These plots show the absolute value of the relative difference between the velocity of the front predicted by hybrid simulations and full stochastic age-dependent SSA simulations (plot (a)) and the coarse-grained mean-field system (plot (b)). The green solid line represents the mean value over time. Each point corresponds to an average over 40 realisations of the age-structured SSA and 100 realisations of the hybrid method. The velocity of the front is calculated using the data corresponding to the position of the (average) front shown in Fig. 4. The threshold, Θ, in the hybrid simulations is taken to be Θ = 2000. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6

Comparison of the velocity average of 100 realisations of hybrid simulation with different threshold Θ = 1000,1500,2000,2500,3000 and 3500. Error bars correspond to standard error of the mean.

Fig. 7

Plots showing the time evolution of probability distribution function of the birth rate as obtained from simulation of the full stochastic multi-scale model. We show three snapshots (time increasing from left to right) for the region behind the interface (upper row) and ahead of the interface (lower row). These results show that behind the interface the distribution of the empirical birth rate, calculated using Eq. (38), is centred around the equilibrium birth rate, Eq. (40), vertical red line. By contrast, ahead of the interface the birth rate distribution is much broader. Therefore, whereas behind the interface the equilibrium birth rate is a good approximation, this is not the case ahead of the interface. See also supplementary movies A.11 and A.12. Θ = 2000. For other parameter values see Section 5. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)