A hybrid discrete-continuum approach to model Turing pattern formation

1.   Introduction

1.1.   Turing's (pre-)pattern theory

In 1952, Alan M. Turing's seminal work 'The chemical basis of morphogenesis' introduced the theory according to which the heterogeneous spatial distribution of some chemicals that regulate cellular differentiation, called "morphogens", acts as a template (i.e., a pre-pattern) for cells to form different kinds of patterns or structures during the embryonic development of an organism ^[1]. In his work, Turing proposed a system of reaction-diffusion equations, with linear reaction terms, modelling the space-time dynamics of the concentrations of two morphogens as the basis for the formation of such pre-patterns. The system had stable homogenous steady states which were driven unstable by diffusion, resulting in spatially heterogeneous morphogen distributions, as long as the diffusion rate of one of the chemical was much larger (order 10) than the other. Twenty years after the publication of Turing's paper, in 1972 Alfred Gierer and Hans Meinhardt further developed this theory through the introduction of activator-inhibitor systems (i.e., reaction-diffusion systems with nonlinear reaction terms) and the notion of "short-range activation and long-range inhibition" ^[2]. The application of this theory to many problems in developmental biology was discussed a further ten years later in 1982, in Meinhardt's book 'Models of Biological Pattern Formation' ^[3], shortly after the specific application of the theory to animal coat markings by James D. Murray ^[4]. Turing (pre-)patterns and resulting cellular patterns have now been discussed widely since their introduction and investigated through further mathematical modelling approaches.

1.2.   Mathematical exploration of cell pattern formation via the Turing mechanism

Several continuum models formulated as systems of partial differential equations (PDEs) have been used to study cell pattern formation via the Turing mechanism, in different spatial dimensions and on domains of various shapes and sizes. Generally, spatial domains can be static or growing to represent the growth of an organism over time. In ^[5], the authors highlighted that there can be a minimum domain size required for pattern formation to occur, and that when considering a growing domain Turing patterns generally become more complex. Multiple authors have analytically and numerically studied pattern formation on growing domains ^{[6,7,8,9,10,11,12,13,14]}. Specifically, in the case where chemotaxis of cells is included (i.e., when cells move up the concentration gradient of the activator), various authors have considered pattern formation mediated by the Turing mechanism on exponentially growing domains ^[15,16]. Along with spatial aspects of cellular patterning, temporal aspects can be considered, such as the role of time-delays in pattern formation. For instance, in ^[17] the authors investigated Turing pattern formation on a morphogen-regulated growing domain where there was a delay in the signalling, gene expression and domain-growth processes.

Turing patterns can arise as stripes, spots (peaks of high density) or reverse-spots (troughs of low density) depending on the particular choice of parameter values and initial distributions of morphogens ^[18]. The different scenarios leading to these three distinct classes of pre-patterns have been explored mathematically by various authors ^[19,20]. For example, in ^[21] the authors showed that, if there is a low level of production of the morphogens, striped patterns are formed by a wider range of parameter settings than spotted patterns. However, Turing patterns can be sensitive to small perturbations in the parameter values and the initial distributions of the morphogens, often leading to a discussion on the robustness of such patterns, or lack thereof ^[5,17]. In regard to a lack of robustness of the Turing mechanism to perturbations in the initial morphogen distributions, it has been found that incorporating stochastic aspects can improve robustness of pattern formation ^[13].

Discrete models and hybrid discrete-continuum models have also been used to describe cell pattern formation via the Turing mechanism in a range of biological and theoretical scenarios ^{[22,23,24,25,26,27,28,29,30]}. In contrast to continuum models formulated as PDEs, such models permit the representation of biological processes at the level of single cells and account for possible stochastic variability in cell dynamics. This leads to greater adaptability and higher accuracy in the mathematical modelling of morphogenesis and pattern formation in biological systems ^[31]. In particular, softwares like CompuCell ^[32] and CompuCell3D ^[33] have been employed to implement hybrid discrete-continuum models to investigate the interplay between gene regulatory networks and cellular processes (e.g., haptotaxis, chemotaxis, cell adhesion and division) during morphogenesis. The three main components of models for cell pattern formation implemented using these softwares are: A Cellular Potts model for the dynamics of the cells and the extracellular matrix; a reaction-diffusion model for the dynamics of the diffusible morphogens; a combination of a state automaton and a set of ordinary differential equations to model the dynamics of gene regulatory networks.

1.3.   A hybrid discrete-continuum approach to model cell pattern formation via the Turing mechanism

Here we develop a discrete-continuum modelling framework for the formation of cellular patterns via the Turing mechanism. In this framework, a reaction-diffusion system for the concentrations of some morphogens is combined with a stochastic individual-based (IB) model that tracks the dynamics of single cells. Such an IB model consists of a set of rules that describe cell movement and proliferation as a discrete-time branching random walk ^[34].

A key advantage of this modelling framework is that it can be easily adapted to both static and growing spatial domains, thus covering a broad spectrum of applications. Moreover, the deterministic continuum limits of the IB models defined in this framework can be formally derived. Such continuum models are formulated as PDEs, which cannot capture phenomena that are driven by stochastic effects in the dynamics of single cells but are more amenable to analytical and numerical approaches. For instance, the numerical simulation of these models requires computational times that are typically much smaller than those required by the numerical exploration of the corresponding IB models for large cell numbers. Continuum models for spatial dynamics of living organisms have been derived from underlying discrete models through different mathematical methods in several previous works. Possible examples include the derivation of continuum models of chemotaxis from velocity-jump process ^{[35,36,37,38,39]} or from different types of random walks ^{[40,44,41,42,43]}; the derivation of diffusion and nonlinear diffusion equations from random walks ^{[45,46,47,48,49,50,51]}, from systems of discrete equations of motion ^{[52,53,54,55,56,57,58]}, from discrete lattice-based exclusion processes ^{[59,60,61,62,63,64,65,66]} and from cellular automata ^[67,68,69]; and the derivation of nonlocal models of cell-cell adhesion from position-jump processes ^[70].

This paper is intended to be a proof of concept for the ideas underlying the modelling framework, with the aim to then apply the related methods to the study of specific patterning and morphogenetic processes, such as those discussed in ^[71,72,73] and references therein, in the future.

1.4.   Contents of the paper

As an illustrative example, we focus on a hybrid discrete-continuum model in which the dynamics of the morphogens are governed by an activator-inhibitor system that gives rise to Turing pre-patterns. The cells then interact with the morphogens in their local area through either of two forms of chemically-dependent cell action: Chemotaxis and chemically-controlled proliferation. We begin by considering such a hybrid model posed on static spatial domains (see Section 2) and then turn to growing domains (see Section 3). Using methods similar to those we have previously employed in ^[44,46], we formally derive the deterministic continuum limit of the model. When sufficiently large numbers of cells are considered, the results of numerical simulations demonstrate an excellent quantitative match between the spatial patterns produced by the stochastic IB model and its deterministic continuum counterpart. Section 4 concludes the paper and provides a brief overview of possible research perspectives.

2.   Mathematical modelling of cell pattern formation on static domains

In this section, we illustrate our hybrid discrete-continuum modelling framework by developing a model for the formation of cellular patterns via the Turing mechanism on static spatial domains (see Section 2.1). The corresponding deterministic continuum model is provided (see Section 2.2) and results of numerical simulations of both models are presented (see Section 2.3). We report on numerical results demonstrating a good match between cellular patterns produced by the stochastic IB model and its deterministic continuum counterpart, in different spatial dimensions and biological scenarios, as well as on results showing the emergence of possible differences between the cell patterns produced by the two models for relatively low cell numbers.

2.1.   A hybrid discrete-continuum model

We let cells and morphogens be distributed across a $d$-dimensional static domain, with $d = 1$ or $d = 2$. In particular, we consider the case where the spatial domain is represented by the interval $[0, \ell]$ when $d = 1$ or the square $[0, \ell] \times [0, \ell]$ when $d = 2$, with $\ell \in \mathbb{R}^*_+$, where $\mathbb{R}^*_+$ is the set of positive real numbers not including zero. The position of the cells and the molecules of morphogens at time $t \in \mathbb{R}_+$ is modelled by the variable $x \in [0, \ell]$ when $d = 1$ and by the vector ${\bf x} = (x, y) \in [0, \ell]^2$ when $d = 2$.

We discretise the time variable $t$ as $t_{k} = k \tau$ with $k\in\mathbb{N}_{0}$ and the space variables $x$ and $y$ as $x_{i} = i \, \chi$ and $y_{j} = j \, \chi$ with $(i, j) \in [0, I]^2 \subset \mathbb{N}^2_0$, where $\tau \in \mathbb{R}^*_+$ and $\chi \in \mathbb{R}^*_+$ are the time- and space-step, respectively, and $I: = 1 + \left\lceil {\frac{\ell }{\chi }} \right\rceil$, where $\left\lceil \cdot \right\rceil$ denotes the ceiling function. Here, $\mathbb{N}_{0}$ is the set of natural numbers including zero. Throughout this section we use the notation ${\bf{i}} \equiv i$ and ${\bf x}_{\bf{i}} \equiv x_i$ when $d = 1$, and ${\bf{i}} \equiv (i, j)$ and ${\bf x}_{\bf{i}} \equiv (x_i, y_j)$ when $d = 2$. The concentrations of the morphogens at position ${\bf x}_{\bf i}$ and at time $t_k$ are modelled by the discrete, non-negative functions $u^{k}_{{\bf i}}$ and $v^{k}_{{\bf i}}$, and we denote by $n^{k}_{{\bf i}}$ the local cell density, which is defined as the number of cells at position ${\bf x}_{\bf{i}}$ and at time $t_k$, which is modelled by the dependent variable $N^{k}_{\bf i}\in\mathbb{N}_0$, divided by the size of the ${\bf i}^{th}$ site of the spatial grid, that is

We present here the model when $d = 2$ but analogous considerations hold for $d = 1$.

2.1.1.   Dynamics of the morphogens

The dynamics of $u^{k}_{{\bf i}}$ and $v^{k}_{{\bf i}}$ are governed by the following coupled system of difference equations

subject to zero-flux boundary conditions. Here, $\delta^2_{i}$ is the second-order central difference operator on the lattice $\left\{x_i\right\}_i$ and $\delta^2_{j}$ is the second-order central difference operator on the lattice $\{y_j\}_j$, that is,

Moreover, $D_{u} \in \mathbb{R}^*_+$ and $D_{v} \in \mathbb{R}^*_+$ represent the diffusion coefficients of the morphogens and the functions $P(u^{k}_{{\bf i}}, v^{k}_{{\bf i}})$ and $Q(u^{k}_{{\bf i}}, v^{k}_{{\bf i}})$ are the rates of change of $u^{k}_{{\bf i}}$ and $v^{k}_{{\bf i}}$ due to local reactions.

The system of difference equations (2.2) is a standard discretisation of a generic reaction-diffusion system of the type that is commonly used to describe morphogen dynamics—see ^[5,74] and references therein. The specific forms of the functions $P$ and $Q$ depend on the reaction kinetics involved in the biological problem at stake—we refer the interested reader to ^[4,5,13] and references therein. We consider an activator-inhibitor system whereby $u^{k}_{{\bf i}}$ models the concentration of the activator while $v^{k}_{{\bf i}}$ models the concentration of the inhibitor. Hence, we let the functions $P$ and $Q$ satisfy the following standard assumptions for activator-inhibitor kinetics

In particular, in this paper we will focus on Schnakenberg kinetics ^[77] and, therefore, the functions $P$ and $Q$ are given via the definitions (2.19). Moreover, we will assume

for some maximal concentrations $u_{\rm max} \in \mathbb{R}^*_+$ and $v_{\rm max} \in \mathbb{R}^*_+$.

2.1.2.   Dynamics of the cells

We consider a scenario where the cells proliferate (i.e., divide and die) and can change their position according to a combination of undirected, random movement and chemotactic movement, which are seen as independent processes. This results in the following rules for the dynamic of the cells.

Mathematical modelling of undirected, random cell movement We model undirected, random cell movement as a random walk with movement probability $\theta \in \mathbb{R}^*_+$, where $0 < \theta\leq1$. In particular, for a cell on the lattice site ${\bf i}$, we define the probability of moving to one of the four neighbouring lattice sites in the von Neumann neighbourhood of ${\bf i}$ via undirected, random movement as $\dfrac{\theta}{4}$ (i.e., there is an equal chance of moving to any of the four neighbouring sites). Therefore, the probability of not undergoing undirected, random movement is defined as $1 - \theta$ (i.e., one minus the probability of movement). Furthermore, the spatial domain is assumed to be closed and, therefore, cell moves that require moving out of the spatial domain are not allowed. Under these assumptions, the probabilities of moving to the left and right sites via undirected, random movement are defined as

while the probabilities of moving to the lower and upper sites via undirected, random movement are defined as

Mathematical modelling of chemotactic cell movement In line with ^[15,16], we let cells move up the concentration gradient of the activator via chemotaxis. Chemotactic cell movement is modelled as a biased random walk whereby the movement probabilities depend on the difference between the concentration of the activator at the site occupied by a cell and the concentration of the activator in the von Neumann neighbourhood of the cell's site. Moreover, as similarly done in the case of undirected, random cell movement, cell moves that require moving out of the spatial domain are not allowed. In particular, building on the modelling strategy presented in ^[44], for a cell on the lattice site ${\bf i}$ and at the time-step $k$, we define the probability of moving to the left and right sites via chemotaxis as

while the probabilities of moving to the lower and upper sites via chemotaxis are defined as

Hence, the probability of not undergoing chemotactic movement is

Here, $(\cdot)_{+}$ denotes the positive part of $(\cdot)$ and the parameter $\eta \in \mathbb{R}_+$ with $0 \leq \eta \leq 1$, where $\mathbb{R}_+$ is the set of non-negative real numbers, is directly proportional to the chemotactic sensitivity of the cells. Notice that since relations (2.5) hold the quantities given via the definitions (2.8)–(2.10) are all between $0$ and $1$.

Mathematical modelling of cell proliferation We consider a scenario in which the cell population undergoes saturating growth. Hence, in line with ^[75], we allow every cell to divide or die with probabilities that depend on a monotonically decreasing function of the local cell density. Moreover, building on the ideas presented in ^[76], we model chemically-controlled cell proliferation by letting the probabilities of cell division and death depend on the local concentrations of the activator and of the inhibitor. In particular, building upon the modelling strategy used in ^[46], between time-steps $k$ and $k+1$, we let a cell on the lattice site ${\bf i}$ divide (i.e., be replaced by two identical daughter cells that are placed on the lattice site ${\bf i}$) with probability

die with probability

or remain quiescent (i.e., do not divide nor die) with probability

Here, $(\cdot)_{+}$ and $(\cdot)_{-}$ denote the positive part and the negative part of $(\cdot)$. The parameters $\alpha_n \in \mathbb{R}^*_+$ and $\beta_n \in \mathbb{R}^*_+$ are, respectively, the intrinsic rates of cell division and cell death. Moreover, the function $\psi$ model the effects of saturating growth and, therefore, it is assumed to be such that

where $n_{\rm max} \in \mathbb{R}^*_+$ is the local carrying capacity of the cell population. Finally, the functions $\phi_u$ and $\phi_v$ satisfy the following assumptions

Notice that we are implicitly assuming that the time-step $\tau$ is sufficiently small that $0 < \mathcal{P}_h < 1$ for all $h \in \{b, d, q \}$.

2.2.   Corresponding continuum model

Letting the time-step $\tau \to 0$ and the spdescribed in the definitions (2.6) and (2.7)ace-step $\chi \to 0$ in such a way that

using the formal method employed in ^[44,46] it is possible to formally show (see Remark A.1 in Appendix A) that the deterministic continuum counterpart of the stochastic IB model presented in Section 2.1, which is described via the systems (2.6)–(2.15), is given by the following PDE for the cell density $n(t, {\bf x})$

subject to zero-flux boundary conditions. Here, $D_{n} \in \mathbb{R}^*_+$ defined via conditions (2.16) is the diffusion coefficient (i.e., the motility) of the cells, while $C_{n} \in \mathbb{R}_+$ defined via conditions (2.16) represents the chemotactic sensitivity of the cells to the activator. In Eq (2.17), the concentration of the activator $u(t, {\bf x})$ and the concentration of the inhibitor $v(t, {\bf x})$ are governed by the continuum counterpart of the system of difference equations (2.2) subject to zero-flux boundary conditions, that is, the following system of PDEs complemented with zero-flux boundary conditions

which can be formally obtained by letting $\tau \to 0$ and $\chi \to 0$ in system (2.2).

2.3.   Numerical simulations

In this section, we carry out a systematic quantitative comparison between the results of numerical simulations of the hybrid model presented in Section 2.1 and numerical solutions of the corresponding continuum model given in Section 2.2, both in one and in two spatial dimensions. All simulations are performed in MATLAB and the final time of simulations is chosen such that the concentrations of morphogens and the cell density are at numerical equilibrium at the end of simulations.

2.3.1.   Summary of the set-up of numerical simulations

Dynamics of the morphogens We consider the case where the functions $P$ and $Q$ that model the rates of change of the concentrations of the morphogens are defined according to Schnakenberg kinetics ^[77], that is,

where $\alpha_{u}, \, \alpha_{v}, \, \beta, \, \gamma \in \mathbb{R}^*_+$. The system of difference equations (2.2) and the system of PDEs (2.18) complemented with definitions (2.19) and subject to zero-flux boundary conditions are known to exhibit Turing pre-patterns. The conditions required for such patterns to emerge have been extensively studied in previous works and, therefore, are omitted here—the interested reader is referred to ^[5] and references therein. We choose parameter values such that Turing pre-patterns arise from the perturbation of homogeneous initial distributions of the morphogens. A complete description of the set-up of numerical simulations is given in Appendix B.

Dynamics of the cells We focus on the case where the cell population undergoes logistic growth and, therefore, we define the function $\psi$ in Eqs (2.11)–(2.13) and Eq (2.17) as

Moreover, we consider two scenarios corresponding to two alternative forms of chemically-dependent cell action. In the first scenario, there is no chemotaxis—i.e., we assume $\eta = 0$ in definitions (2.8) and (2.9), which implies that $C_{n} = 0$ in Eq (2.17)—and the cells interact with the morphogens in their local area through chemically-controlled proliferation. In particular, we use the following definitions of the functions $\phi_u$ and $\phi_v$ in Eqs (2.11)–(2.13) and Eq (2.17)

In the second scenario, chemotaxis up the concentration gradient of the activator occurs—i.e., we assume $\eta > 0$, which implies that $C_{n} > 0$—but cell division and death are not regulated by the morphogens—i.e., we assume

In both scenarios, we let the initial cell distribution be homogeneous and, given the values of the parameters chosen to carry out numerical simulations of the IB model, we use the following definitions

so that conditions (2.16) are met. A complete description of the set-up of numerical simulations is given in Appendix B.

2.3.2.   Main results of numerical simulations

Dynamics of the morphogens The plots in the top rows of Figures 1 and 2 and in the Supplementary Figure D2 summarise the dynamics of the continuum concentrations of morphogens $u(t, {\bf x})$ and $v(t, {\bf x})$ obtained by solving numerically the system of PDEs (2.18) subject to zero-flux boundary conditions. Identical results hold for the discrete morphogen concentrations $u_{\bf i}^k$ and $v_{\bf i}^k$ obtained by solving the system of difference equations (2.2) (results not shown). These plots demonstrate that in the case where the reaction terms $P$ and $Q$ are described via the definitions (2.19), under the parameter setting considered here, Turing pre-patterns arise from perturbation of homogeneous initial distributions of the morphogens. Such pre-patterns consist of spots of activator and inhibitor, whereby maximum points of the concentration of activator coincide with minimum points of the concentration of inhibitor, and vice versa.

Dynamics of the cells in the presence of chemically-controlled cell proliferation The plots in the bottom row of Figure 1 and the plots in Figure 3 summarise the dynamics of the cell density in the case where there is no chemotaxis and chemically-controlled cell proliferation occurs—i.e., when $\eta = 0$, $C_n = 0$, and the functions $\phi_u$ and $\phi_v$ are given via the definitions (2.21). Since a larger concentration of the activator corresponds to a higher cell division rate and a smaller concentration of the inhibitor corresponds to a lower cell death rate, we observe the formation of cellular patterns consisting of spots of cells centred at the same points as the spots of activator. These plots demonstrate that there is an excellent quantitative match between the discrete cell density $n^k_{\bf i}$ given by Eq (2.1), with $N^k_{\bf i}$ obtained through computational simulations of the IB model, and the continuum cell density $n(t, {\bf x})$ obtained by solving numerically the PDE (2.17) subject to zero-flux boundary conditions, both in one and in two spatial dimensions.

Dynamics of the cells in the presence of chemotaxis The plots in the bottom row of Figure 2 and the plots in Figure 4 summarise the dynamics of the cell density in the case where cell proliferation is not regulated by the morphogens and chemotactic movement of the cells up the concentration gradient of the activator occurs—i.e., when the functions $\phi_u$ and $\phi_v$ are given via the definitions (2.22), $\eta > 0$ and $C_n > 0$. Since the cells sense the concentration of the activator and move up its gradient, cellular patterns consisting of spots of cells centred at the same points as the spots of the activator are formed. Compared to the case of chemically-controlled cell proliferation, in this case the spots of cells are smaller and characterised by a larger cell density (i.e., cells are more densely packed). There is again an excellent quantitative match between the discrete cell density $n^k_{\bf i}$ given by Eq (2.1), with $N^k_{\bf i}$ obtained through computational simulations of the IB model, and the continuum cell density $n(t, {\bf x})$ obtained by solving numerically the PDE (2.17) subject to zero-flux boundary conditions, both in one and in two spatial dimensions.

Emergence of possible differences between cell patterns produced by the IB model and the continuum model In all cases discussed so far we have observed excellent agreement between the dynamics of the discrete cell density obtained through computational simulations of the stochastic IB model and the continuum cell density obtained by solving numerically the corresponding deterministic continuum model. However, we expect possible differences between the two models to emerge in the presence of low cell numbers. In order to investigate this, we carried out numerical simulations of the IB model and the PDE model for the case where cells undergo chemically-controlled cell proliferation, considering either lower initial cell densities along with lower values of the local carrying capacity of the cell population $n_{\rm max}$ or higher rates of cell death $\beta_n$, which correspond to lower saturation values of the local cell density. The plots in the bottom row of Figure 5 and in Figure 6 summarise the dynamics of the cell density for relatively small initial cell numbers and local carrying capacities. These plots show that differences between the discrete cell density $n^k_{\bf i}$ given by Eq (2.1), with $N^k_{\bf i}$ obtained through computational simulations of the IB model, and the continuum cell density $n(t, {\bf x})$, obtained by solving numerically the PDE (2.17) subject to zero-flux boundary conditions, can emerge both in one and in two spatial dimensions. Analogous considerations hold for the case in which higher rates of cell death $\beta_n$ are considered (results not shown).

3.   Mathematical modelling of cell pattern formation on growing domains

So far, we have considered the case of static spatial domains. However, in many biological problems the formation of cellular patterns occurs on growing spatial domains, for example, in morphogenesis and embryogenesis ^{[10,11,12,13]}. Therefore, in this section, we extend the hybrid model developed in the previous section to the case of growing spatial domains (see Section 3.1). We consider both the case of uniform domain growth and the case of apical growth (i.e., the case where domain growth is restricted to a region located toward the boundary). Similarly to the previous section, the deterministic continuum limit of the model is provided (see Section 3.2) and the results of numerical simulations demonstrating a good match between the cellular patterns produced by the stochastic IB model and its deterministic continuum counterpart are presented (see Section 3.3).

3.1.   A hybrid-discrete continuum model

Building upon the modelling framework presented in the previous section, we let cells and morphogens be distributed across a $d$-dimensional growing domain represented by the interval $[0, \mathcal{L}(t)]$ when $d = 1$ and the square $[0, \mathcal{L}(t)]^2$ when $d = 2$. The real, positive function $\mathcal{L}(t)$, with $\mathcal{L}(\cdot) \geq 1$, determines the growth of the right-hand and upper boundary of the spatial domain (i.e., we consider the case where the domain grows equally in both the $x$ and $y$ directions). In analogy with the previous section, we use the notation $x \in [0, \mathcal{L}(t)]$ and ${\bf x} = (x, y) \in [0, \mathcal{L}(t)]^2$. Moreover, we make the change of variables

which allows one to describe the spatial position of the cells and the molecules of morphogens by means of the variable $\hat{x}\in[0, 1]$ when $d = 1$ and the vector $\hat{{\bf x}} = (\hat{x}, \hat{y}) \in [0, 1]^2$ when $d = 2$ ^[8,78]. We then discretise the time variable $t$ and space variables $\hat x$ and $\hat y$ in a similar way to the static domain case, as described in Section 2.1. Throughout this section we use the notation ${\bf{i}} \equiv i$ and $\hat {\bf x}_{\bf{i}} \equiv \hat x_i$ when $d = 1$, and ${\bf{i}} \equiv (i, j)$ and $\hat {\bf x}_{\bf{i}} \equiv (\hat x_i, \hat y_j)$ when $d = 2$. We also use the notation $\mathcal{L}_k = \mathcal{L}(t_k)$. The concentrations of the morphogens at position $\hat {\bf x}_{\bf i}$ and at time $t_k$ are modelled by the discrete, non-negative functions $u^{k}_{{\bf i}}$ and $v^{k}_{{\bf i}}$, and we denote by $n^{k}_{{\bf i}}$ the local cell density, which is defined via Eq (2.1). As in the case of static domains, we present the model for $d = 2$ but analogous considerations hold for $d = 1$.

3.1.1.   Dynamics of the morphogens

The dynamics of $u^{k}_{{\bf i}}$ and $v^{k}_{{\bf i}}$ are governed by the following coupled system of difference equations

subject to zero-flux boundary conditions. Here, $\delta^2_{i}$ is the second-order central difference operator on the lattice $\left\{\hat x_{i}\right\}_{i}$ and $\delta^2_{j}$ is the second-order central difference operator on the lattice $\{\hat y_{j}\}_{j}$, which are provided in definitions (2.3). Moreover, $D_{u} \in \mathbb{R}^*_+$ and $D_{v} \in \mathbb{R}^*_+$ represent the diffusion coefficients of the morphogens, which are rescaled by $\mathcal{L}^2_k$ for consistency with the change of variables (3.1), and the functions $P(u^{k}_{{\bf i}}, v^{k}_{{\bf i}})$ and $Q(u^{k}_{{\bf i}}, v^{k}_{{\bf i}})$ are the rates of change of $u^{k}_{{\bf i}}$ and $v^{k}_{{\bf i}}$ due to local reactions, which satisfy conditions (2.4) and (2.5), as in the case of static domains. Finally, the last terms on the right-hand sides of system (3.2) represent the rates of change of the concentrations of morphogens due to variation in the size of the domain. In the case of uniform domain growth, the following definition holds ^[8,9]

Equation (3.3) captures the effects of dilution of the concentrations of the morphogens due to local volume changes of the spatial domain ^[8,11]. On the other hand, when apical growth of the domain occurs one has ^[9,78]

3.1.2.   Dynamics of the cells

Under the change of variables (3.1), the dynamics of the cells in the IB model is governed by rules analogous to those described in Section 2.1.2 for the case of static domains. In summary, definitions (2.6) and (2.7) are modified as

Moreover, definitions (2.8) and (2.9) are modified as

Finally, definitions (2.11)–(2.13) are modified as

and

Here, the function $g_{{\bf i}}(n^{k}_{{\bf i}}, \mathcal{L}_{k})$ is defined via Eq (3.3) in the case of uniform domain growth and via Eq (3.4) in the case of apical growth. The functions $\psi$, $\phi_u$ and $\phi_v$ satisfy assumptions (2.14) and (2.15), and we assume $\tau$ and $\mathcal{L}_{k}$ to be such that that $0 < \mathcal{P}_h < 1$ for all $h \in \{b, d, q \}$.

3.2.   Corresponding continuum model

Similarly to the case of static domains, letting the time-step $\tau \to 0$ and the space-step $\chi \to 0$ in such a way that conditions (2.16) are met, it is possible to formally show (see Appendix A) that the deterministic continuum counterpart of the stochastic IB model on growing domains is given by the following PDE for the cell density $n(t, \hat{{\bf x}})$

subject to zero-flux boundary conditions, with either

in the case of uniform domain growth, or

in the case of apical growth. Here, $D_{n} \in \mathbb{R}^*_+$ in definitions (2.16) is the rescaled diffusion coefficient (i.e., the rescaled motility) of the cells, while $C_{n} \in \mathbb{R}_+$ in definitions (2.16) represents the chemotactic sensitivity of cells to the activator. In Eq (3.12), the concentration of the activator $u(t, \hat{{\bf x}})$ and the concentration of the inhibitor $v(t, \hat{{\bf x}})$ are governed by the continuum counterpart of the difference equations (3.2) complemented with zero-flux boundary conditions, that is, the following system of PDEs subject to zero-flux boundary conditions

which can be formally obtained by letting $\tau \to 0$ and $\chi \to 0$ in system (3.2).

3.3.   Numerical simulations

In this section, we carry out a systematic quantitative comparison between the results of numerical simulations of the hybrid model presented in Section 3.1 and numerical solutions of the corresponding continuum model given in Section 3.2, both in one and in two spatial dimensions. All simulations are performed in MATLAB and the final time of simulations is chosen such that the essential features of the pattern formation process are evident.

3.3.1.   Summary of the set-up of numerical simulations

We define the functions $P$, $Q$, $\psi$, $\phi_u$ and $\phi_v$ as in the case of static domains. In more detail, $P$ and $Q$ are provided in definitions (2.19), $\psi$ is defined via Eq (2.20), and $\phi_u$ and $\phi_v$ are given via either definitions (2.21) or definitions (2.22).

In all simulations, we let the initial spatial distributions of morphogens and cells be the numerical steady state distributions obtained in the case of static domains with $\ell: = 1$, and we assume the domain to slowly grow linearly over time, that is,

Given the values of the parameters chosen to carry out numerical simulations of the IB model, we describe $D_{n}$ and $C_{n}$ via the definitions (2.23) so that conditions (2.16) are met. A complete description of the set-up of numerical simulations is given in Appendix C.

3.3.2.   Main results of numerical simulations

Dynamics of the morphogens The plots in the top rows of Figures 7 and 8 and in the Supplementary Figure D2 summarise the dynamics of the continuum concentrations of morphogens $u(t, \hat {\bf x})$ and $v(t, \hat {\bf x})$ obtained by solving numerically the system of PDEs (3.15) subject to zero-flux boundary conditions and with $G(\hat {\bf x}, u, \mathcal{L})$ and $G(\hat {\bf x}, v, \mathcal{L})$ defined via Eq (3.13), while the plots in the top rows of Figures 9 and 10 and in the Supplementary Figure D3 refer to the case where $G(\hat {\bf x}, u, \mathcal{L})$ and $G(\hat {\bf x}, v, \mathcal{L})$ are defined via Eq (3.14). Identical results hold for the discrete morphogen concentrations $u_{\bf i}^k$ and $v_{\bf i}^k$ obtained by solving the system of difference equations (3.2) (results not shown). These plots demonstrate that, when the spatial domain grows over time, a dynamical rescaling and repetition of the Turing pre-patterns observed in the case of static domains occurs—i.e., spots of high concentration repeatedly split symmetrically. In the case of uniform domain growth, such a self-similar process occurs throughout the whole domain, while in the case of apical growth the process takes place toward the growing edge of the domain.

Dynamics of the cells The plots in the bottom row of Figure 7 and the plots in Figure 11 summarise the dynamics of the cell density in the case where there is no chemotaxis, chemically-controlled cell proliferation occurs—i.e., when $\eta = 0$, $C_n = 0$, and the functions $\phi_u$ and $\phi_v$ are given by definitions (2.21)—and the functions $g_{{\bf i}}(n^{k}_{{\bf i}}, \mathcal{L}_{k})$ and $G(\hat {\bf x}, n, \mathcal{L})$ are defined via Eqs (3.3) and (3.13), respectively. On the other hand, the plots in the bottom row of Figure 9 and the plots in Figure 12 refer to the case where the functions $g_{{\bf i}}(n^{k}_{{\bf i}}, \mathcal{L}_{k})$ and $G(\hat {\bf x}, n, \mathcal{L})$ are defined via Eqs (3.4) and (3.14). These plots indicate that, when the spatial domain grows over time, spots of high cell density stretch either throughout the domain (uniform growth) or at the growing edge (apical growth) before splitting. This process causes cell patterns to rescale and repeat across the domain at a smaller scale. These plots demonstrate also that there is a good quantitative match between the discrete cell density $n^k_{\bf i}$ given by Eq (2.1), with $N^k_{\bf i}$ obtained through computational simulations of the IB model, and the continuum cell density $n(t, \hat {\bf x})$ obtained by solving numerically the PDE (3.12) subject to zero-flux boundary conditions and complemented with either Eq (3.13) or Eq (3.14), both in one and in two spatial dimensions.

Analogous considerations apply to the case where cell proliferation is not regulated by the morphogens and chemotactic movement of the cells up the concentration gradient of the activator occurs—i.e., when the functions $\phi_u$ and $\phi_v$ are described through the definitions (2.22), $\eta > 0$ and $C_n > 0$—see the plots in the bottom row of Figure 8 along with the plots in Figure 13 and the plots in the bottom row of Figure 10 along with the plots in Figure 14.

4.   Conclusions and research perspectives

4.1.   Conclusions

We have developed a hybrid discrete-continuum modelling framework that can be used to describe the formation of cellular patterns, specifically focusing on the Turing mechanism as the driving force behind the patterns. We used reaction-diffusion systems to describe the evolution of morphogens, which dictate the action of cells, while cell dynamics were described by stochastic IB models. We formally derived the deterministic continuum counterparts of the IB models, which were formulated as PDEs for the cell density, and compared the two modelling approaches through numerical simulations both in the case of stationary spatial domains and in the case of two types of growing domains, corresponding to uniform and apical growth. Numerical simulations demonstrated that in the case of sufficiently large cell numbers there was an excellent quantitative match between the spatial patterns produced by the stochastic IB model and its deterministic continuum counterpart. Moreover, in the case of static domains, we also presented the results of numerical simulations showing that possible differences between the spatial patterns produced by the two modelling approaches could emerge in the regime of sufficiently low cell numbers. In fact, lower cell numbers correlated with both lower regularity of the cell density and demographic stochasticity, which may cause a reduction in the quality of the approximations employed in the formal derivation of the deterministic continuum model from the stochastic IB model. Hence, having both types of models available allows one to use IB models in the regime of low cells numbers—i.e., when stochastic effects associated with small cell population levels, which cannot be captured by PDE models, are particularly relevant—and then turn to their less computationally expensive PDE counterparts when large cell numbers need to be considered—i.e., when stochastic effects associated with small cell population levels are negligible.

4.2.   Research perspectives

There are a number of additional elements that would be relevant to incorporate into the modelling framework presented here in order to further broaden its spectrum of applications.

For instance, as was recognised by Turing himself, exogenous diffusing chemicals are not the only vehicle of coordination between cells. In particular, it is known that long range cell-cell interactions can be mediated by signal proteins produced by the cells themselves and also by mechanical forces between cells and components of the cellular microenvironment. For example, vascular endothelial growth factor signalling has been shown to control neural crest cell migration ^[79,80,81], and mechanical interactions between cells and the extra cellular matrix can control cell agammaegation ^[82]. Moreover, cellular patterning leading to the emergence of spatial structures often requires the interplay between non-diffusible species, transcription factors and cell signalling—viz. the process underlying digit formation in tetrapods ^[83]. In this regard, it would be interesting to extend the modelling framework by allowing the cells to consume and/or produce chemicals required for successful coordination of their actions ^[84], and by incorporating more complex cellular processes such as anoikis ^[85,86] and cell deformation ^[87,88]. In the situation where local production and/or consumption of the chemicals by the cells occurs, for particular cases such as those considered in ^[44], we would still expect it to be possible to derive an effective deterministic continuum limit of the IB model for the dynamics of the cells through formal procedures analogous to the one used here. However, there could also be cases in which PDE models derived using similar formal procedures might not be able to faithfully reproduce the dynamics of the branching random walk underlying the IB model, due to the interplay between stochastic effects and nonlinear dynamical interactions between the cells and the chemicals.

To date, only few biological systems have been demonstrated to satisfy the necessary conditions required for the formation of Turing pre-patterns via reaction-diffusion systems. Since mathematical models formulated as scalar integro-differential equations, whereby the formation of Turing-like patterns is governed by suitable integral kernels, have proven capable of faithfully reproduce a variety of pigmentation patterns in fish ^[27,89], it would also be interesting to explore possible ways of integrating such alternative modelling strategies into our framework.

Acknowledgments

MAJC gratefully acknowledges support of EPSRC Grant No. EP/N014642/1 (EPSRC Centre for Multiscale Soft Tissue Mechanics–With Application to Heart & Cancer).

Conflict of interest

All authors declare no conflicts of interest in this paper.

Appendix

A.   Formal derivation of the deterministic continuum model on growing domains

We carry out a formal derivation of the deterministic continuum model given by the PDE (2.16) for $d = 2$. Similar methods can be used in the case where $d = 1$.

When cell dynamics are governed by the rules described in Section 2.1.2 and Section 3.1.2, considering $(i, j) \in [1, I-1] \times [1, I-1]$, the mass balance principle gives

Using the fact that the following relations hold for $\tau$ and $\chi$ sufficiently small

the balance equation (A.1) can be formally rewritten in the approximate form

with

where $n \equiv n(t, \hat x, \hat y)$, $u \equiv u(t, \hat x, \hat y)$, $v \equiv v(t, \hat x, \hat y)$ and $\mathcal{L} \equiv \mathcal{L}(t)$. Dividing both sides of Eq (A.2) by $\tau$ gives

If $n(t, \hat x, \hat y)$ is a twice continuously differentiable function of $\hat x$ and $\hat y$ and a continuously differentiable function of $t$, $u(t, \hat x, \hat y)$ is a twice continuously differentiable function of $\hat x$ and $\hat y$, and the function $\mathcal{L}(t)$ is continuously differentiable, for $\chi$ and $\tau$ sufficiently small we can use the Taylor expansions

Substituting into Eq (A.3), using the elementary property $(a)_+ - (-a)_+ = a$ for $a \in \mathbb{R}$ and letting $\tau \to 0$ and $\chi \to 0$ in such a way that conditions (2.16) are met, after a little algebra, as similarly done in ^[44], we find

where $G(\hat x, \hat y, n, \mathcal{L})$ is given by Eq (3.13) in the case where $g_{(i, j)}(n^{k}_{(i, j)})$ is defined via Eq (3.3) and by equation (3.14) in the case where $g_{(i, j)}(n^{k}_{(i, j)})$ is defined via Eq (3.4). The PDE (A.4) can be easily rewritten in the form of Eq (3.12). Moreover, zero-flux boundary conditions easily follow from the fact that [cf. definitions (3.5)–(3.8)]

and

Remark A.1. The derivation of the continuum limit for the static domain case can be carried out in a similar way by assuming $\mathcal{L}_{k}$ to be constant, which implies that $g_{{\bf i}} \equiv 0$ and results in $G \equiv 0$.

B.   Set-up of numerical simulations on static domains

We let $x\in[0, 1]$, $y\in[0, 1]$ and $\chi: = 0.005$ (i.e., $I = 201$). Moreover, we define $\tau : = 1\times 10^{-3}$.

Dynamics of the morphogens For the dynamics of the morphogens, we consider the parameter setting used in ^[13], that is,

Moreover, we assume the initial distributions to be small perturbations of the homogeneous steady state $(u^*, v^*)\equiv(1, 0.9)$, that is,

where $\rho : = 0.001$ and ${\bf R}$ is either a vector for $d = 1$ or a matrix for $d = 2$ whose components are random numbers drawn from the standard uniform distribution on the interval $(0, 1)$, using the built-in MATLAB function RAND. These choices of the initial distributions of morphogens are such that

that is, the parameter $\rho$ determines the level of perturbation from the homogeneous steady state. Since the difference equations (2.2) governing the dynamics of the morphogens are independent from the dynamics of the cells, such equations are solved first for all time-steps and the solutions obtained are then used to evaluate both the probabilities of cell movement [cf. definitions (2.6)–(2.9)] and the probabilities of cell division and death [cf. definitions (2.11)–(2.13)]. The parameter $u_{\rm max}$ in definitions (2.8) and (2.9) is defined as ${\max_{k, {\bf i}} u_{{\bf i}}^k}$.

Computational implementation of the rules underlying the dynamics of the cells At each time-step, each cell undergoes a three-phase process: Phase 1) undirected, random movement according to the probabilities described in the definitions (2.6) and (2.7); Phase 2) chemotaxis according to the probabilities described via the definitions (2.8) and (2.9); Phase 3) division and death according to the probabilities defined via Eqs (2.11)–(2.13). For each cell, during each phase, 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.

Dynamics of the cells Unless stated otherwise, we assume the initial cell distributions to be homogeneous with

In the case where chemically-controlled cell proliferation occurs and there is no chemotaxis, unless stated otherwise, we use the following parameter values when $d = 1$

and the following ones when $d = 2$

The results shown in Figures 5 and 6 refer to the same settings with the modification that when $d = 1$

and when $d = 2$

In the case where cells undergo chemotaxis and cell proliferation is not chemically-controlled, unless stated otherwise, we use the following parameter values when $d = 1$

and the following ones when $d = 2$

Numerical solutions of the corresponding continuum models Numerical solutions of the PDE (2.17) and the system of PDEs (2.18) subject to zero-flux boundary conditions are computed through standard finite-difference schemes using initial conditions and parameter values that are compatible with those used for the IB model and the system of difference equations (2.2). In particular, the values of the parameters $D_n$ and $C_n$ in the PDE (2.17) are described via the definitions (2.23).

C.   Set-up of numerical simulations on growing domains

We let $x\in[0, 1]$, $y\in[0, 1]$ and $\chi: = 0.005$ (i.e., $I = 201$). Moreover, we assume $\tau : = 1\times 10^{-3}$ and we define $\mathcal{L}$ according to equation (3.16) (i.e., the domain grows linearly over time).

Dynamics of the morphogens For the dynamics of the morphogens, we use the parameter setting given by the definitions (13.1). Moreover, we define the initial distributions as the numerical equilibrium distributions obtained in the case of static domains. Similarly to the case of static domains, since the difference equations (3.2) governing the dynamics of the morphogens are independent from the dynamics of the cells, such equations are solved first for all time-steps and the solutions obtained are then used to evaluate both the probabilities of cell movement given by definitions (3.5)–(3.8) and the probabilities of cell division and death given by Eqs (3.9)–(3.11). The parameter $u_{\rm max}$ in the definitions (3.7) and (3.8) is defined as ${\max_{k, {\bf i}} u_{{\bf i}}^k}$.

Computational implementation of the rules underlying the dynamics of the cells Similarly to the case of static domains, at each time-step, each cell undergoes a three-phase process: Phase 1) undirected, random movement according to the probabilities described through the definitions (3.5) and (3.6); Phase 2) chemotaxis according to the probabilities described through the definitions (3.7) and (3.8); Phase 3) division and death according to the probabilities defined via Eqs (3.9)–(3.11). For each cell, during each phase, 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.

Dynamics of the cells We assume the initial cell distributions and all parameter values to be the same as those used in the static domain case.

Numerical solutions of the corresponding continuum models Numerical solutions of the PDE (3.12) and the system of PDEs (3.15) subject to zero-flux boundary conditions are computed through standard finite-difference schemes using initial conditions and parameter values that are compatible with those used for the IB model and the system of difference equations (3.2). In particular, the values of the parameters $D_n$ and $C_n$ in the PDE (3.12) are described through the definitions (2.23).

D.   Supplementary figures