A hybrid discrete-continuum approach to model Turing pattern formation

Since its introduction in 1952, Turing's (pre-)pattern theory ("the chemical basis of morphogenesis") has been widely applied to a number of areas in developmental biology. The related pattern formation models normally comprise a system of reaction-diffusion equations for interacting chemical species ("morphogens"), whose heterogeneous distribution in some spatial domain acts as a template for cells to form some kind of pattern or structure through, for example, differentiation or proliferation induced by the chemical pre-pattern. Here we develop a hybrid discrete-continuum modelling framework for the formation of cellular patterns via the Turing mechanism. In this framework, a stochastic individual-based model of cell movement and proliferation is combined with a reaction-diffusion system for the concentrations of some morphogens. As an illustrative example, we focus on a 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 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, and then turn to the case of growing domains. In both cases, we formally derive the corresponding deterministic continuum limit and show that that there is an excellent quantitative match between the spatial patterns produced by the stochastic individual-based model and its deterministic continuum counterpart, when sufficiently large numbers of cells are considered. This paper is intended to present 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 in the future.


Introduction
Turing's (pre-)pattern theory In 1952, A.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 [49]. 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 A. Gierer and H. 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" [13]. 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' [33], shortly after the specific application of the theory to animal coat markings by J.D. Murray [36]. Turing (pre-)patterns and resulting cellular patterns have now been discussed widely since their introduction and investigated through further mathematical modelling approaches.
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 [27], 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 [1,4,6,7,19,21,22,28,26]. 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 [25,44]. 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 [23] 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 [37]. The different scenarios leading to these three distinct classes of pre-patterns have been explored mathematically by various authors [10,48]. For example, in [34] 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 [23,27]. 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 [28].
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 [9,17,18,35,41,51]. 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 [14]. In particular, softwares like CompuCell [16] and CompuCell3D [5] 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.
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 [15].
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, the numerical simulation of which requires computational times that are typically much smaller than those required by the numerical exploration of the corresponding IB models for large cell numbers. Hence, having both types of models available allows one to use IB models in the regime of low cells numbersi.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 consideredi.e., when stochastic effects associated with small cell population levels are negligible.
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 [29,42,43] and references therein, in the future.
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). In both cases, we formally derive the deterministic continuum limit of the model, using methods similar to those we previously employed in [2,3], and show that that there is an excellent quantitative match between the spatial patterns produced by the stochastic IB model and its deterministic continuum counterpart, when sufficiently large numbers of cells are considered. In the case of static domains, we also present the results of numerical simulations showing that possible differences between the spatial patterns produced by the two modelling approaches can emerge in the regime of sufficiently low cell numbers. In fact, lower cell numbers correlate 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. Section 4 concludes the paper providing a brief overview of possible research perspectives.

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.

A hybrid discrete-continuum model
We let cells and morphogens be distributed across a d-dimensional static domain. In particular, we consider the case where the spatial domain is represented by the interval [0, ] when d = 1 and the square [0, ] × [0, ] when d = 2, with ∈ R * + . The position of the cells and the molecules of morphogens at time t ∈ R + is modelled by the variable x ∈ [0, ] when d = 1 and by the vector x = (x, y) ∈ [0, ] 2 when d = 2.
We discretise the time variable t as t k = kτ with k ∈ N 0 and the space variables x and y as x i = i χ and y j = j χ with (i, j) ∈ [0, I] 2 ⊂ N 2 0 , where τ ∈ R * + and χ ∈ R * + are the time-and space-step, respectively, and I := 1 + χ .
Throughout this section we use the notation i ≡ i and x i ≡ x i when d = 1, and i ≡ (i, j) and x i ≡ (x i , y j ) when d = 2. The concentrations of the morphogens at position x i and at time t k are modelled by the discrete, non-negative functions u k i and v k i , and we denote by n k i the local cell density, which is defined as where the dependent variable N k i ∈ N 0 models the number of cells at position x i and at time t k . We present here the model when d = 2 but analogous considerations hold for d = 1.

Dynamics of the morphogens
The dynamics of u k i and v k i are governed by the following coupled system of difference equations subject to zero-flux boundary conditions. Here, δ 2 i is the second-order central difference operator on the lattice {x i } i and δ 2 j is the second-order central difference operator on the lattice {y j } j , that is, and Moreover, D u ∈ R * + and D v ∈ R * + represent the diffusion coefficients of the morphogens and the functions P(u k i , v k i ) and Q(u k i , v k i ) are the rates of change of u k i and v k i due to local reactions. We consider an activator-inhibitor system whereby u k i models the concentration of the activator while v k i models the concentration of the inhibitor. Hence, we let the functions P and Q satisfy the following assumptions ∂P ∂v < 0 and ∂Q ∂u > 0 (2.4) and be such that for some maximal concentrations u max ∈ R * + and v max ∈ R * + .

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 θ ∈ R * + , where 0 < θ ≤ 1. In particular, for a cell on the lattice site i, we define the probability of moving to one of the lattice sites in the von Neumann neighbourhood of i via undirected, random movement as θ 4 , while the probability of not undergoing undirected, random movement is defined as 1 − θ. 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 Mathematical modelling of chemotactic cell movement In line with [25,44], 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 [2], for a cell on the lattice site 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 (2.10) Here, (·) + denotes the positive part of (·) and the parameter η ∈ R + , with 0 ≤ η ≤ 1, is directly proportional to the chemotactic sensitivity of the cells. Notice that since (2.5) holds the quantities defined via (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 [39], 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 [47], 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 [3], between time-steps k and k + 1, we let a cell on the lattice site i divide (i.e., be replaced by two identical daughter cells that are placed on the lattice site i) with probability 12) or remain quiescent (i.e., do not divide nor die) with probability Here, (·) + and (·) − denote the positive part and the negative part of (·). The parameters α n ∈ R * + and β n ∈ R * + are, respectively, the intrinsic rates of cell division and cell death. Moreover, the function ψ model the effects of saturating growth and, therefore, it is assumed to be such that ψ (·) < 0 and ψ(n max ) = 0, (2.14) where n max ∈ R * + is the local carrying capacity of the cell population. Finally, the functions φ u and φ v satisfy the following assumptions Notice that we are implicitly assuming that the time-step τ is sufficiently small that 0 < P h < 1 for all h ∈ {b, d, q}.

Corresponding continuum model
Letting the time-step τ → 0 and the space-step χ → 0 in such a way that using the formal method employed in [2,3] it is possible to formally show that the deterministic continuum counterpart of the stochastic IB model presented in Section 2.1 is given by the following PDE for the cell density n(t, x) subject to zero-flux boundary conditions. Here, D n ∈ R * + defined via (2.16) is the diffusion coefficient (i.e., the motility) of the cells, while C n ∈ R + defined via (2.16) represents the chemotactic sensitivity of the cells to the activator. In (2.17), the concentration of the activator u(t, x) and the concentration of the inhibitor v(t, 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 τ → 0 and χ → 0 in (2.2).

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.

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 [45], that is, where α u , α v , β, γ ∈ R * + . The system of difference equations (2.2) and the system of PDEs (2.18) complemented with (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 [27] and references therein. We choose parameter values such that Turing prepatterns 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 ψ in (2.11)-(2.13) and (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 chemotaxisi.e., we assume η = 0 in (2.8) and (2.9), which implies that C n = 0 in (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 φ u and φ v in (2.11)-(2.13) and (2.17) In the second scenario, chemotaxis up the concentration gradient of the activator occursi.e., we assume η > 0, which implies that C n > 0 -but cell division and death are not regulated by the morphogensi.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   Figure D1 summarise the dynamics of the continuum concentrations of morphogens u(t, x) and v(t, 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 k i and v k i 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 defined via (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 2 summarise the dynamics of the cell density in the case where there is no chemotaxis and chemically-controlled cell proliferation occursi.e., when η = 0, C n = 0, and the functions φ u and φ v are defined via (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 i given by (2.1), with N k i obtained through computational simulations of the IB model, and the continuum cell density n(t, 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 3 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 occursi.e., when the functions φ u and φ v are defined via (2.22), η > 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 i given by (2.1), with N k i obtained through computational simulations of the IB model, and the continuum cell density n(t, 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 max or higher rates of cell death β 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 i given by (2.1), with N k i obtained through computational simulations of the IB model, and the continuum cell density n(t, 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 β n are considered (results not shown).

Mathematical modelling of cell pattern formation on growing domains
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).  Figure D1. Here, the parameter setting is the same as that of Figure 2 but with a smaller initial cell density and a smaller local carrying capacity of the cell population n max . A complete description of the set-up of numerical simulations is given in Appendix B. d = 2. The real, positive function L(t), with L(·) ≥ 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 ∈ [0, L(t)] and x = (x, y) ∈ [0, L(t)] 2 . Moreover, we make the change of variables [6,24] x →x and x →x = (x,ŷ) withx := x L(t)

A hybrid-discrete continuum model
,ŷ := y L(t) We discretise the time variable t as t k = kτ with k ∈ N 0 and the space variablesx andŷ asx i = i χ andŷ j = j χ with (i, j) ∈ [0, I] 2 ⊂ N 2 0 , where τ ∈ R * + and χ ∈ R * + are the time-and space-step, respectively, and I := 1 + 1 χ . Throughout this section we use the notation i ≡ i andx i ≡x i when d = 1, and i ≡ (i, j) andx i ≡ (x i ,ŷ j ) when d = 2. We also use the notation L k = L(t k ). The concentrations of the morphogens at positionx i and at time t k are modelled by the discrete, non-negative functions u k i and v k i , and we denote by n k i the local cell density, which is defined as where the dependent variable N k i ∈ N 0 models the number of cells at positionx i and at time t k . As in the case of static domains, we present the model for d = 2 but analogous considerations hold for d = 1.

Dynamics of the morphogens
The dynamics of u k i and v k i are governed by the following coupled system of difference equations subject to zero-flux boundary conditions. Here, δ 2 i is the second-order central difference operator on the lattice {x i } i and δ 2 j is the second-order central difference operator on the lattice {ŷ j } j , which are defined via (2.3). Moreover, D u ∈ R * + and D v ∈ R * + represent the diffusion coefficients of the morphogens, which are rescaled by L 2 k for consistency with the change of variables (3.1), and the functions P(u k i , v k i ) and Q(u k i , v k i ) are the rates of change of u k i and v k i due to local reactions, which satisfy assumptions (2.4) and (2.5), as in the case of static domains. Finally, the last terms on the right-hand sides of (3.3) 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 [6,7] Definiton (3.4) captures the effects of dilution of the concentrations of the morphogens due to local volume changes of the spatial domain [6,21]. On the other hand, when apical growth of the domain occurs one has [7,24]  Moreover, definitions (2.8) and (2.9) are modified as

Dynamics of the cells
Finally, definitions (2.11)-(2.13) are modified as and Here, the function g i (n k i , L k ) is defined via (3.4) in the case of uniform domain growth and via (3.5) in the case of apical growth. The functions ψ, φ u and φ v satisfy assumptions (2.14) and (2.15), and we assume τ and L k to be such that that 0 < P h < 1 for all h ∈ {b, d, q}.

Corresponding continuum model
Similarly to the case of static domains, letting the time-step τ → 0 and the space-step χ → 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,x) subject to zero-flux boundary conditions, with either (3.14) in the case of uniform domain growth, or (3.15) in the case of apical growth. Here, D n ∈ R * + defined via (2.16) is the rescaled diffusion coefficient (i.e., the rescaled motility) of the cells, while C n ∈ R + defined via (2.16) represents the chemotactic sensitivity of cells to the activator. In (3.13), the concentration of the activator u(t,x) and the concentration of the inhibitor v(t,x) are governed by the continuum counterpart of the difference equations (3.3) 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 τ → 0 and χ → 0 in (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.

Summary of the set-up of numerical simulations
We define the functions P, Q, ψ, φ u and φ v as in the case of static domains. In more detail, P and Q are defined via (2.19), ψ is defined via (2.20), and φ u and φ v are defined via either (2.21) or (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 := 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 define D n and C n via (2.23) so that conditions (2.16) are met. A complete description of the set-up of numerical simulations is given in Appendix C.

Main results of numerical simulations
Dynamics of the morphogens The plots in the top rows of Figures 7 and 9 and in the Supplementary Figure D2 summarise the dynamics of the continuum concentrations of morphogens u(t,x) and v(t,x) obtained by solving numerically the system of PDEs (3.16) subject to zero-flux boundary conditions and with G(x, u, L) and G(x, v, L) defined via (3.14), while the plots in the top rows of Figures 11 and 13 and in the Supplementary Figure D3 refer to the case where G(x, u, L) and G(x, v, L) are defined via (3.15). Identical results hold for the discrete morphogen concentrations u k i and v k i obtained by solving the system of difference equations (3.3) (results not shown). These plots demonstrate that, when the spatial domain grows over time, a dynamical rescaling and repetition of the Turing prepatterns observed in the case of static domains occursi.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. Figure 7 and the plots in Figure 8 summarise the dynamics of the cell density in the case where there is no chemotaxis, chemically-controlled cell proliferation occursi.e., when η = 0, C n = 0, and the functions φ u and φ v are defined via (2.21) -and the functions g i (n k i , L k ) and G(x, n, L) are defined via (3.4) and (3.14), respectively. On the other hand, the plots in the bottom row of Figure 11 and the plots in Figure 12 refer to the case where the functions g i (n k i , L k ) and G(x, n, L) are defined via (3.5) and (3.15). 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 i given by (3.2), with N k i obtained through computational simulations of the IB model, and the continuum cell density n(t,x) obtained by solving numerically the PDE (3.13) subject to zero-flux boundary conditions and complemented with either (3.14) or (3.15), 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 occursi.e., when the functions φ u and φ v are defined via (2.22), η > 0 and C n > 0 -see the plots in the bottom row of Figure 9 along with the plots in Figure 10 and the plots in the bottom row of Figure 13 along with the plots in Figure 14.

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 [30,31,32], and mechanical interactions between cells and the extra cellular matrix can control cell aggregation [38]. Moreover, cellular patterning leading to the emergence of spatial structures often requires the interplay between non-diffusible species, transcription factors and cell signallingviz. the process underlying digit formation in tetrapods [46]. 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 [50], and by incorporating more complex cellular processes such as anoikis [12,11] and cell deformation [8,40].       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 [18,20], it would also be interesting to explore possible ways of integrating such alternative modelling strategies into our framework.
Dynamics of the morphogens For the dynamics of the morphogens, we consider the parameter setting used in [28], that is, Moreover, we assume the initial distributions to be small perturbations of the homogeneous steady state (u * , v * ) ≡ (1, 0.9), that is, where ρ := 0.001 and 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 u * − ρ ≤ u 0 i ≤ u * + ρ and v * − ρ ≤ v 0 i ≤ v * + ρ for all i, that is, the parameter ρ 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 (2.6)-(2.9) and the probabilities of cell division and death (2.11)-(2.13). The parameter u max in (2.8) and (2.9) is defined as max 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 defined via (2.6) and (2.7); Phase 2) chemotaxis according to the probabilities defined via (2.8) and (2.9); Phase 3) division and death according to the probabilities defined via (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 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 defined via (2.23).
Dynamics of the morphogens For the dynamics of the morphogens, we use the parameter setting given by (B.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.3) 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 (3.6)-(3.9) and the probabilities of cell division and death (3.10)- (3.12). The parameter u max in (3.8) and (3.9) is defined as max k,i u k i .
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 defined via (3.6) and (3.7); Phase 2) chemotaxis according to the probabilities defined via (3.8) and (3.9); Phase 3) division and death according to the probabilities defined via (3.10)-(3.12). 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.13) and the system of PDEs (3.16) subject to zero-flux boundary conditions are computed thorugh 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.3). In particular, the values of the parameters D n and C n in the PDE (3.13) are defined via (2.23).