A hybrid discrete-continuum approach to model Turing pattern formation
AA hybrid discrete-continuum approach to model Turing patternformation
Dedicated to the memory of Federica Bubba
Fiona R Macfarlane , Mark AJ Chaplain and Tommaso Lorenzi School of Mathematics and Statistics, University of St Andrews, UK; Department of Mathematical Sciences “G. L. Lagrange”, Dipartimento di Eccellenza2018-2022, Politecnico di Torino, 10129 Torino, Italy;
Abstract
Since its introduction in 1952, with a further refinement in 1972 by Gierer and Meinhardt, Turing’s (pre-)patterntheory (“the chemical basis of morphogenesis”) has been widely applied to a number of areas in developmentalbiology, where evolving cell and tissue structures are naturally observed. The related pattern formation modelsnormally comprise a system of reaction-di ff usion equations for interacting chemical species (“morphogens”), whoseheterogeneous distribution in some spatial domain acts as a template for cells to form some kind of pattern or structurethrough, for example, di ff erentiation or proliferation induced by the chemical pre-pattern. Here we develop a hybriddiscrete-continuum modelling framework for the formation of cellular patterns via the Turing mechanism. In thisframework, a stochastic individual-based model of cell movement and proliferation is combined with a reaction-di ff usion system for the concentrations of some morphogens. As an illustrative example, we focus on a model inwhich 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 hybridmodel posed on static spatial domains, and then turn to the case of growing domains. In both cases, we formally derivethe corresponding deterministic continuum limit and show that that there is an excellent quantitative match betweenthe spatial patterns produced by the stochastic individual-based model and its deterministic continuum counterpart,when su ffi ciently large numbers of cells are considered. This paper is intended to present a proof of concept for theideas underlying the modelling framework, with the aim to then apply the related methods to the study of specificpatterning and morphogenetic processes in the future. Turing’s (pre-)pattern theory
In 1952, A.M. Turing’s seminal work ‘
The chemical basis of morphogenesis ’ introducedthe theory according to which the heterogeneous spatial distribution of some chemicals that regulate cellular di ff erentiation,called “morphogens”, acts as a template ( i.e., a pre-pattern) for cells to form di ff erent kinds of patterns or structuresduring the embryonic development of an organism [87]. In his work, Turing proposed a system of reaction-di ff usionequations, with linear reaction terms, modelling the space-time dynamics of the concentrations of two morphogensas the basis for the formation of such pre-patterns. The system had stable homogenous steady states which weredriven unstable by di ff usion, resulting in spatially heterogeneous morphogen distributions, as long as the di ff usionrate of one of the chemical was much larger (order 10) than the other. Twenty years after the publication of Turing’spaper, in 1972 A. Gierer and H. Meinhardt further developed this theory through the introduction of activator-inhibitorsystems ( i.e. , reaction-di ff usion systems with nonlinear reaction terms) and the notion of “short-range activationand long-range inhibition” [26]. The application of this theory to many problems in developmental biology wasdiscussed a further ten years later in 1982, in Meinhardt’s book ‘ Models of Biological Pattern Formation ’ [56], shortlyafter the specific application of the theory to animal coat markings by J.D. Murray [60]. Turing (pre-)patterns andresulting cellular patterns have now been discussed widely since their introduction and investigated through furthermathematical modelling approaches. 1 a r X i v : . [ q - b i o . CB ] S e p athematical exploration of cell pattern formation via the Turing mechanism Several continuum models formulatedas systems of partial di ff erential equations (PDEs) have been used to study cell pattern formation via the Turingmechanism, in di ff erent spatial dimensions and on domains of various shapes and sizes. Generally, spatial domainscan be static or growing to represent the growth of an organism over time. In [49], the authors highlighted that therecan be a minimum domain size required for pattern formation to occur, and that when considering a growing domainTuring patterns generally become more complex. Multiple authors have analytically and numerically studied patternformation on growing domains [2, 11, 15, 16, 38, 40, 41, 50, 48]. Specifically, in the case where chemotaxis of cells isincluded ( i.e. , when cells move up the concentration gradient of the activator), various authors have considered patternformation mediated by the Turing mechanism on exponentially growing domains [47, 74]. Along with spatial aspectsof cellular patterning, temporal aspects can be considered, such as the role of time-delays in pattern formation. Forinstance, in [43] the authors investigated Turing pattern formation on a morphogen-regulated growing domain wherethere 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) dependingon the particular choice of parameter values and initial distributions of morphogens [61]. The di ff erent scenariosleading to these three distinct classes of pre-patterns have been explored mathematically by various authors [22, 82].For example, in [57] the authors showed that, if there is a low level of production of the morphogens, striped patternsare formed by a wider range of parameter settings than spotted patterns. However, Turing patterns can be sensitive tosmall perturbations in the parameter values and the initial distributions of the morphogens, often leading to a discussionon the robustness of such patterns, or lack thereof [43, 49]. In regard to a lack of robustness of the Turing mechanismto perturbations in the initial morphogen distributions, it has been found that incorporating stochastic aspects canimprove robustness of pattern formation [50].Discrete models and hybrid discrete-continuum models have also been used to describe cell pattern formationvia the Turing mechanism in a range of biological and theoretical scenarios [12, 13, 20, 35, 36, 37, 58, 69, 89]. Incontrast to continuum models formulated as PDEs, such models permit the representation of biological processesat the level of single cells and account for possible stochastic variability in cell dynamics. This leads to greateradaptability and higher accuracy in the mathematical modelling of morphogenesis and pattern formation in biologicalsystems [27]. In particular, softwares like CompuCell [32] and CompuCell3D [14] have been employed to implementhybrid 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 modelsfor cell pattern formation implemented using these softwares are: a Cellular Potts model for the dynamics of the cellsand the extracellular matrix; a reaction-di ff usion model for the dynamics of the di ff usible morphogens; a combinationof a state automaton and a set of ordinary di ff erential equations to model the dynamics of gene regulatory networks. A hybrid discrete-continuum approach to model cell pattern formation via the Turing mechanism
Here wedevelop a discrete-continuum modelling framework for the formation of cellular patterns via the Turing mechanism. Inthis framework, a reaction-di ff usion system for the concentrations of some morphogens is combined with a stochasticindividual-based (IB) model that tracks the dynamics of single cells. Such an IB model consists of a set of rules thatdescribe cell movement and proliferation as a discrete-time branching random walk [30].A key advantage of this modelling framework is that it can be easily adapted to both static and growing spatialdomains, thus covering a broad spectrum of applications. Moreover, the deterministic continuum limits of the IBmodels defined in this framework can be formally derived. Such continuum models are formulated as PDEs, whichcannot capture phenomena that are driven by stochastic e ff ects in the dynamics of single cells but are more amenableto analytical and numerical approaches. For instance, the numerical simulation of these models requires computationaltimes that are typically much smaller than those required by the numerical exploration of the corresponding IB modelsfor large cell numbers. Continuum models for spatial dynamics of living organisms have been derived from underlyingdiscrete models thorugh di ff erent mathematical methods in several previous works. Possible examples include thederivation of continuum models of chemotaxis from velocity-jump process [28, 29, 71, 75, 76] or from di ff erenttypes of random walks [1, 5, 6, 85, 86]; the derivation of di ff usion and nonlinear di ff usion equations from randomwalks [9, 10, 31, 67, 70, 77, 78], from systems of discrete equations of motion [3, 8, 45, 59, 63, 64, 68], from discretelattice-based exclusion processes [4, 21, 23, 33, 34, 42, 46, 83] and from cellular automata [17, 19, 84]; and thederivation of nonlocal models of cell-cell adhesion from position-jump processes [7].This paper is intended to be a proof of concept for the ideas underlying the modelling framework, with the aim tothen apply the related methods to the study of specific patterning and morphogenetic processes, such as those discussedin [52, 72, 73] and references therein, in the future. 2 ontents of the paper As an illustrative example, we focus on a hybrid discrete-continuum model in which thedynamics of the morphogens are governed by an activator-inhibitor system that gives rise to Turing pre-patterns. Thecells then interact with the morphogens in their local area through either of two forms of chemically-dependent cellaction: chemotaxis and chemically-controlled proliferation. We begin by considering such a hybrid model posed onstatic spatial domains (see Section 2) and then turn to growing domains (see Section 3). Using methods similar tothose we have previously employed in [5, 10], we formally derive the deterministic continuum limit of the model.When su ffi ciently large numbers of cells are considered, the results of numerical simulations demonstrate an excellentquantitative match between the spatial patterns produced by the stochastic IB model and its deterministic continuumcounterpart. Section 4 concludes the paper and provides a brief overview of possible research perspectives. In this section, we illustrate our hybrid discrete-continuum modelling framework by developing a model for theformation of cellular patterns via the Turing mechanism on static spatial domains (see Section 2.1). The correspondingdeterministic continuum model is provided (see Section 2.2) and results of numerical simulations of both models arepresented (see Section 2.3). We report on numerical results demonstrating a good match between cellular patternsproduced by the stochastic IB model and its deterministic continuum counterpart, in di ff erent spatial dimensions andbiological scenarios, as well as on results showing the emergence of possible di ff erences between the cell patternsproduced by the two models for relatively low cell numbers. We let cells and morphogens be distributed across a d -dimensional static domain, with d = d =
2. In particular, weconsider the case where the spatial domain is represented by the interval [0 , (cid:96) ] when d = , (cid:96) ] × [0 , (cid:96) ]when d =
2, with (cid:96) ∈ R ∗ + , where R ∗ + is the set of positive real numbers not including zero. The position of the cellsand the molecules of morphogens at time t ∈ R + is modelled by the variable x ∈ [0 , (cid:96) ] when d = x = ( x , y ) ∈ [0 , (cid:96) ] when d = t as t k = k τ with k ∈ N and the space variables x and y as x i = i χ and y j = j χ with ( i , j ) ∈ [0 , I ] ⊂ N , where τ ∈ R ∗ + and χ ∈ R ∗ + are the time- and space-step, respectively, and I : = + (cid:38) (cid:96)χ (cid:39) , where (cid:100)·(cid:101) denotes the ceiling function. Here, N is the set of natural numbers including zero. Throughout this section weuse 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 themorphogens 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 denoteby n k i the local cell density, which is defined as the number of cells at position x i and at time t k , which is modelled bythe dependent variable N k i ∈ N , divided by the size of the i th site of the spatial grid, that is n k i : = N k i χ d . (2.1)We present here the model when d = d = The dynamics of u k i and v k i are governed by the following coupled system of di ff erence equations u k + i = u k i + τ D u χ (cid:16) δ i u k i + δ j u k i (cid:17) + τ P ( u k i , v k i ) , v k + i = v k i + τ D v χ (cid:16) δ i v k i + δ j v k i (cid:17) + τ Q ( u k i , v k i ) , ( k , i ) ∈ N × (0 , I ) , (2.2)subject to zero-flux boundary conditions. Here, δ i is the second-order central di ff erence operator on the lattice { x i } i and δ j is the second-order central di ff erence operator on the lattice { y j } j , that is, δ i u k i : = u k ( i + , j ) + u k ( i − , j ) − u k ( i , j ) and δ j u k i : = u k ( i , j + + u k ( i , j − − u k ( i , j ) . (2.3)3oreover, D u ∈ R ∗ + and D v ∈ R ∗ + represent the di ff usion coe ffi cients 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.The system of di ff erence equations (2.2) is a standard discretisation of a generic reaction-di ff usion system of thetype that is commonly used to describe morphogen dynamics – see, for instance, [49, 51] and references therein. Thespecific forms of the functions P and Q depend on the reaction kinetics involved in the biological problem at stake – werefer the interested reader to [49, 50, 60] and references therein. 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 thefunctions P and Q satisfy the following standard assumptions for activator-inhibitor kinetics ∂ P ∂ v < ∂ Q ∂ u > . (2.4)In particular, in this paper we will focus on Schnakenberg kinetics [79] and, therefore, we will define the functions P and Q via (2.19). Moreover, we will assume0 < u k i ≤ u max and 0 < v k i ≤ v max ∀ ( k , i ) ∈ N × [0 , I ] (2.5)for some maximal concentrations u max ∈ R ∗ + and v max ∈ R ∗ + . We consider a scenario where the cells proliferate ( i.e. , divide and die) and can change their position according to acombination 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 asa 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 four neighbouring lattice sites in the von Neumann neighbourhoodof i via undirected, random movement as θ i.e. there is an equal chance of moving to any of the four neighbouringsites). Therefore, the probability of not undergoing undirected, random movement is defined as 1 − θ ( i.e. one minusthe probability of movement). Furthermore, the spatial domain is assumed to be closed and, therefore, cell moves thatrequire moving out of the spatial domain are not allowed. Under these assumptions, the probabilities of moving to theleft and right sites via undirected, random movement are defined as T k L( i , j ) : = θ , T k R( i , j ) : = θ i , j ) ∈ [1 , I − × [0 , I ] , (2.6) T k L(0 , j ) : = , T k R( I , j ) : = j ∈ [0 , I ] , while the probabilities of moving to the lower and upper sites via undirected, random movement are defined as T k D( i , j ) : = θ , T k U( i , j ) : = θ i , j ) ∈ [0 , I ] × [1 , I − , (2.7) T k D( i , : = , T k U( i , I ) : = i ∈ [0 , I ] . Mathematical modelling of chemotactic cell movement
In line with [47, 74], we let cells move up the concentrationgradient of the activator via chemotaxis. Chemotactic cell movement is modelled as a biased random walk wherebythe movement probabilities depend on the di ff erence between the concentration of the activator at the site occupiedby a cell and the concentration of the activator in the von Neumann neighbourhood of the cell’s site. Moreover, assimilarly done in the case of undirected, random cell movement, cell moves that require moving out of the spatialdomain are not allowed. In particular, building on the modelling strategy presented in [5], for a cell on the lattice site4 and at the time-step k , we define the probability of moving to the left and right sites via chemotaxis as J k L( i , j ) : = η (cid:16) u k ( i − , j ) − u k ( i , j ) (cid:17) + u max , J k R( i , j ) : = η (cid:16) u k ( i + , j ) − u k ( i , j ) (cid:17) + u max for ( i , j ) ∈ [1 , I − × [0 , I ] , (2.8) J k L(0 , j ) : = , J k R( I , j ) : = j ∈ [0 , I ] , while the probabilities of moving to the lower and upper sites via chemotaxis are defined as J k D( i , j ) : = η (cid:16) u k ( i , j − − u k ( i , j ) (cid:17) + u max , J k U( i , j ) : = η (cid:16) u k ( i , j + − u k ( i , j ) (cid:17) + u max for ( i , j ) ∈ [0 , I ] × [1 , I − , (2.9) J k D( i , : = , J k U( i , I ) : = i ∈ [0 , I ] . Hence, the probability of not undergoing chemotactic movement is1 − (cid:16) J k L i + J k R i + J k D i + J k U i (cid:17) for i ∈ [0 , I ] . (2.10)Here, ( · ) + denotes the positive part of ( · ) and the parameter η ∈ R + with 0 ≤ η ≤
1, where R + is the set of non-negativereal numbers, is directly proportional to the chemotactic sensitivity of the cells. Notice that since (2.5) holds thequantities 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 undergoessaturating growth. Hence, in line with [65], we allow every cell to divide or die with probabilities that depend ona monotonically decreasing function of the local cell density. Moreover, building on the ideas presented in [81], wemodel chemically-controlled cell proliferation by letting the probabilities of cell division and death depend on the localconcentrations of the activator and of the inhibitor. In particular, building upon the modelling strategy used in [10],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 cellsthat are placed on the lattice site i ) with probability P b (cid:16) n k i , u k i (cid:17) : = τ α n (cid:16) ψ ( n k i ) (cid:17) + φ u ( u k i ) , (2.11)die with probability P d (cid:16) n k i , u k i v k i (cid:17) : = τ (cid:16) α n (cid:16) ψ ( n k i ) (cid:17) − φ u ( u k i ) + β n φ v ( v k i ) (cid:17) , (2.12)or remain quiescent ( i.e. , do not divide nor die) with probability P q (cid:16) n k i , u k i , v k i (cid:17) : = − τ (cid:16) α n (cid:12)(cid:12)(cid:12) ψ ( n k i ) (cid:12)(cid:12)(cid:12) φ u ( u k i ) + β n φ v ( v k i ) (cid:17) . (2.13)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 e ff ects of saturatinggrowth and, therefore, it is assumed to be such that ψ (cid:48) ( · ) < ψ ( n max ) = , (2.14)where n max ∈ R ∗ + is the local carrying capacity of the cell population. Finally, the functions φ u and φ v satisfy thefollowing assumptions φ u (0) = , φ (cid:48) u ( · ) > φ v (0) = , φ (cid:48) v ( · ) > . (2.15)Notice that we are implicitly assuming that the time-step τ is su ffi ciently small that 0 < P h < h ∈ { b , d , q } . Letting the time-step τ → χ → θ χ d τ → D n and η χ d τ u max → C n as τ → , χ → , (2.16)5sing the formal method employed in [5, 10] it is possible to formally show (see Remark A in Appendix A) that thedeterministic continuum counterpart of the stochastic IB model presented in Section 2.1, which is defined via (2.6) -(2.15), is given by the following PDE for the cell density n ( t , x ) ∂ t n − ∇ x · ( D n ∇ x n − C n n ∇ x u ) = (cid:16) α n ψ ( n ) φ u ( u ) − β n φ v ( v ) (cid:17) n , ( t , x ) ∈ R ∗ + × (0 , (cid:96) ) d (2.17)subject to zero-flux boundary conditions. Here, D n ∈ R ∗ + defined via (2.16) is the di ff usion coe ffi cient ( i.e. , themotility) of the cells, while C n ∈ R + defined via (2.16) represents the chemotactic sensitivity of the cells to theactivator. In (2.17), the concentration of the activator u ( t , x ) and the concentration of the inhibitor v ( t , x ) are governedby the continuum counterpart of the system of di ff erence equations (2.2) subject to zero-flux boundary conditions, thatis, the following system of PDEs complemented with zero-flux boundary conditions ∂ t u − D u ∆ x u = P ( u , v ) ,∂ t v − D v ∆ x u = Q ( u , v ) , ( t , x ) ∈ R ∗ + × (0 , (cid:96) ) d , (2.18)which can be formally obtained by letting τ → χ → In this section, we carry out a systematic quantitative comparison between the results of numerical simulations ofthe hybrid model presented in Section 2.1 and numerical solutions of the corresponding continuum model given inSection 2.2, both in one and in two spatial dimensions. All simulations are performed in M atlab and the final time ofsimulations is chosen such that the concentrations of morphogens and the cell density are at numerical equilibrium atthe end of simulations.
We consider the case where the functions P and Q that model the rates of change ofthe concentrations of the morphogens are defined according to Schnakenberg kinetics [79], that is, P ( u , v ) : = α u − β u + γ u v , Q ( u , v ) : = α v − γ u v (2.19)where α u , α v , β, γ ∈ R ∗ + . The system of di ff erence equations (2.2) and the system of PDEs (2.18) complementedwith (2.19) and subject to zero-flux boundary conditions are known to exhibit Turing pre-patterns. The conditionsrequired for such patterns to emerge have been extensively studied in previous works and, therefore, are omitted here– the interested reader is referred to [49] 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 descriptionof 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, wedefine the function ψ in (2.11)-(2.13) and (2.17) as ψ ( n ) : = (cid:32) − nn max (cid:33) . (2.20)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 η = C n = φ u and φ v in (2.11)-(2.13) and (2.17) φ u ( u ) : = + uu max and φ v ( v ) : = + vv max . (2.21)In the second scenario, chemotaxis up the concentration gradient of the activator occurs – i.e. , we assume η >
0, whichimplies that C n > i.e. , we assume φ u ( u ) ≡ φ v ( v ) ≡ . (2.22)6n both scenarios, we let the initial cell distribution be homogeneous and, given the values of the parameters chosento carry out numerical simulations of the IB model, we use the following definitions D n : = θ χ d τ and C n : = η χ d τ u max (2.23)so that conditions (2.16) are met. A complete description of the set-up of numerical simulations is given in Appendix B. The plots in the top rows of Figures 1 and 3 and in the Supplementary Figure D1summarise the dynamics of the continuum concentrations of morphogens u ( t , x ) and v ( t , x ) obtained by solvingnumerically the system of PDEs (2.18) subject to zero-flux boundary conditions. Identical results hold for thediscrete morphogen concentrations u k i and v k i obtained by solving the system of di ff erence equations (2.2) (resultsnot shown). These plots demonstrate that in the case where the reaction terms P and Q are defined via (2.19), underthe parameter setting considered here, Turing pre-patterns arise from perturbation of homogeneous initial distributionsof the morphogens. Such pre-patterns consist of spots of activator and inhibitor, whereby maximum points of theconcentration 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 ofFigure 1 and the plots in Figure 2 summarise the dynamics of the cell density in the case where there is no chemotaxisand chemically-controlled cell proliferation occurs – i.e. , when η = C n =
0, and the functions φ u and φ v aredefined via (2.21). Since a larger concentration of the activator corresponds to a higher cell division rate and asmaller concentration of the inhibitor corresponds to a lower cell death rate, we observe the formation of cellularpatterns consisting of spots of cells centred at the same points as the spots of activator. These plots demonstrate thatthere is an excellent quantitative match between the discrete cell density n k i given by (2.1), with N k i obtained throughcomputational simulations of the IB model, and the continuum cell density n ( t , x ) obtained by solving numerically thePDE (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 inFigure 4 summarise the dynamics of the cell density in the case where cell proliferation is not regulated by themorphogens and chemotactic movement of the cells up the concentration gradient of the activator occurs – i.e. , whenthe functions φ u and φ v are defined via (2.22), η > C n >
0. Since the cells sense the concentration of theactivator and move up its gradient, cellular patterns consisting of spots of cells centred at the same points as thespots of the activator are formed. Compared to the case of chemically-controlled cell proliferation, in this case thespots of cells are smaller and characterised by a larger cell density ( i.e. , cells are more densely packed). There isagain an excellent quantitative match between the discrete cell density n k i given by (2.1), with N k i obtained throughcomputational simulations of the IB model, and the continuum cell density n ( t , x ) obtained by solving numerically thePDE (2.17) subject to zero-flux boundary conditions, both in one and in two spatial dimensions. Emergence of possible di ff erences 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 densityobtained through computational simulations of the stochastic IB model and the continuum cell density obtained bysolving numerically the corresponding deterministic continuum model. However, we expect possible di ff erencesbetween the two models to emerge in the presence of low cell numbers. In order to investigate this, we carriedout numerical simulations of the IB model and the PDE model for the case where cells undergo chemically-controlledcell proliferation, considering either lower initial cell densities along with lower values of the local carrying capacityof the cell population n max or higher rates of cell death β n , which correspond to lower saturation values of the localcell density. The plots in the bottom row of Figure 5 and in Figure 6 summarise the dynamics of the cell densityfor relatively small initial cell numbers and local carrying capacities. These plots show that di ff erences between thediscrete cell density n k i given by (2.1), with N k i obtained through computational simulations of the IB model, andthe continuum cell density n ( t , x ), obtained by solving numerically the PDE (2.17) subject to zero-flux boundaryconditions, can emerge both in one and in two spatial dimensions. Analogous considerations hold for the case inwhich higher rates of cell death β n are considered (results not shown).7igure 1: Results of numerical simulations on a one-dimensional static domain in the presence of chemically-controlled cell proliferation . Top row.
Plots of the concentrations of morphogens at four consecutive time instants.The green lines highlight the concentration of activator u ( t , x ) and the red lines highlight the concentration of inhibitor v ( t , x ) obtained by solving numerically the system of PDEs (2.18) for d = Bottom row.
Comparison between the discrete cell density n ki obtained by averagingthe results of computational simulations of the IB model (solid dark blue lines) and the continuum cell density n ( t , x )obtained by solving numerically the PDE (2.17) for d = η = C n =
0, and the functions φ u and φ v are defined via (2.21). Weadditionally set the initial cell density n i = for all i . The results from the IB model correspond to the average overfive realisations of the underlying branching random walk, with the results from each realisation plotted in pale blueto demonstrate the robustness of the results obtained. A complete description of the set-up of numerical simulations isgiven in Appendix B. 8igure 2: Results of numerical simulations on a two-dimensional static domain in the presence of chemically-controlled cell proliferation . Comparison between the discrete cell density n k i obtained by averaging the resultsof computational simulations of the IB model (top row) and the continuum cell density n ( t , x ) obtained by solvingnumerically the PDE (2.17) for d = η = C n =
0, and the functions φ u and φ v are defined via (2.21). We additionally set the initialcell density n i = × for all i . The results from the IB model correspond to the average over five realisations ofthe underlying branching random walk. The plots of the corresponding morphogen concentrations are displayed in theSupplementary Figure D1. A complete description of the set-up of numerical simulations is given in Appendix B.9igure 3: Results of numerical simulations on a one-dimensional static domain in the presence of chemotaxis . Top row.
Plots of the concentrations of morphogens at four consecutive time instants. The green lines highlight theconcentration of activator u ( t , x ) and the red lines highlight the concentration of inhibitor v ( t , x ) obtained by solvingnumerically the system of PDEs (2.18) for d =
1, complemented with (2.19) and subject to zero-flux boundaryconditions.
Bottom row.
Comparison between the discrete cell density n ki obtained by averaging the results ofcomputational simulations of the IB model (solid dark blue lines) and the continuum cell density n ( t , x ) obtainedby solving numerically the PDE (2.17) for d = η > C n >
0, and the functions φ u and φ v are defined via (2.22). We additionallyset the initial cell density n i = for all i . The results from the IB model correspond to the average over fiverealisations of the underlying branching random walk, with the results from each realisation plotted in pale blue todemonstrate the robustness of the results obtained. A complete description of the set-up of numerical simulations isgiven in Appendix B. 10igure 4: Results of numerical simulations on a two-dimensional static domain in the presence of chemically-controlled cell proliferation . Comparison between the discrete cell density n k i obtained by averaging the resultsof computational simulations of the IB model (top row) and the continuum cell density n ( t , x ) obtained by solvingnumerically the PDE (2.17) for d = η > C n >
0, and the functions φ u and φ v are defined via (2.22). We additionally set the initialcell density n i = × for all i . The results from the IB model correspond to the average over five realisations ofthe underlying branching random walk. The plots of the corresponding morphogen concentrations are displayed in theSupplementary Figure D1. A complete description of the set-up of numerical simulations is given in Appendix B.11igure 5: Emergence of possible di ff erences between cell patterns produced by the IB model and the continuummodel for low cell numbers on a one-dimensional static domain . Top row.
Plots of the concentrations ofmorphogens at four consecutive time instants. The green lines highlight the concentration of activator u ( t , x ) andthe red lines highlight the concentration of inhibitor v ( t , x ) obtained by solving numerically the system of PDEs (2.18)for d =
1, complemented with (2.19) and subject to zero-flux boundary conditions.
Bottom row.
Comparison betweenthe discrete cell density n ki obtained by averaging the results of computational simulations of the IB model (solid darkblue lines) and the continuum cell density n ( t , x ) obtained by solving numerically the PDE (2.17) for d = η = C n =
0, and thefunctions φ u and φ v are defined via (2.21). We additionally set the initial cell density n i = × for all i . The resultsfrom the IB model correspond to the average over five realisations of the underlying branching random walk, withthe results from each realisation plotted in pale blue. The parameter setting is the same as that of Figure 1 but with asmaller initial cell density and a smaller local carrying capacity of the cell population n max . A complete description ofthe set-up of numerical simulations is given in Appendix B.12igure 6: Emergence of possible di ff erences between cell patterns produced by the IB model and the continuummodel for low cell numbers on a two-dimensional static domain . Comparison between the discrete cell density n k i obtained by averaging the results of computational simulations of the IB model (top row) and the continuum celldensity n ( t , x ) obtained by solving numerically the PDE (2.17) for d = η = C n =
0, and the functions φ u and φ v are definedvia (2.21). We additionally set the initial cell density n i = × for all i . The results from the IB model correspond tothe average over five realisations of the underlying branching random walk. The plots of the corresponding morphogenconcentrations are displayed in the Supplementary Figure D1. The parameter setting is the same as that of Figure 2but with a smaller initial cell density and a smaller local carrying capacity of the cell population n max . A completedescription of the set-up of numerical simulations is given in Appendix B.13 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 ofcellular patterns occurs on growing spatial domains, for example, in morphogenesis and embryogenesis [38, 40, 41,50]. Therefore, in this section, we extend the hybrid model developed in the previous section to the case of growingspatial 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 previoussection, the deterministic continuum limit of the model is provided (see Section 3.2) and the results of numericalsimulations demonstrating a good match between the cellular patterns produced by the stochastic IB model and itsdeterministic continuum counterpart are presented (see Section 3.3).
Building upon the modelling framework presented in the previous section, we let cells and morphogens be distributedacross a d -dimensional growing domain represented by the interval [0 , L ( t )] when d = , L ( t )] when d =
2. The real, positive function L ( t ), with L ( · ) ≥
1, determines the growth of the right-hand and upper boundaryof the spatial domain ( i.e. , we consider the case where the domain grows equally in both the x and y directions). Inanalogy with the previous section, we use the notation x ∈ [0 , L ( t )] and x = ( x , y ) ∈ [0 , L ( t )] . Moreover, we make thechange of variables [15, 44] x (cid:55)→ ˆ x and x (cid:55)→ ˆ x = ( ˆ x , ˆ y ) with ˆ x : = x L ( t ) , ˆ y : = y L ( t ) (3.1)which allows one to describe the spatial position of the cells and the molecules of morphogens by means of the variableˆ x ∈ [0 ,
1] when d = x = ( ˆ x , ˆ y ) ∈ [0 , when d =
2. We then discretise the time variable t and spacevariables ˆ x and ˆ y in a similar way to the static domain case, as described in Section 2.1. Throughout this sectionwe 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. We also use thenotation L k = L ( t k ). 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 via (2.1). As in the caseof static domains, we present the model for d = d = The dynamics of u k i and v k i are governed by the following coupled system of di ff erence equations u k + i = u k i + τ D u L k χ (cid:16) δ i u k i + δ j u k i (cid:17) + τ P ( u k i , v k i ) − g i ( u k i , L k ) , v k + i = v k i + τ D v L k χ (cid:16) δ i v k i + δ j v k i (cid:17) + τ Q ( u k i , v k i ) − g i ( v k i , L k ) , ( k , i ) ∈ N × (0 , I ) , (3.2)subject to zero-flux boundary conditions. Here, δ i is the second-order central di ff erence operator on the lattice { ˆ x i } i and δ j is the second-order central di ff erence operator on the lattice { ˆ y j } j , which are defined via (2.3). Moreover, D u ∈ R ∗ + and D v ∈ R ∗ + represent the di ff usion coe ffi cients of the morphogens, which are rescaled by L k for consistency withthe 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 localreactions, which satisfy assumptions (2.4) and (2.5), as in the case of static domains. Finally, the last terms on theright-hand sides of (3.2) represent the rates of change of the concentrations of morphogens due to variation in the sizeof the domain. In the case of uniform domain growth, the following definition holds [15, 16] g i ( w k i , L k ) ≡ g ( w k i , L k ) : = d w k i L k + − L k L k . (3.3)Definiton (3.3) captures the e ff ects of dilution of the concentrations of the morphogens due to local volume changesof the spatial domain [15, 40]. On the other hand, when apical growth of the domain occurs one has [16, 44] g i ( w k i , L k ) : = (cid:104) i (cid:16) w ki + , j − w ki , j (cid:17) + j (cid:16) w ki , j + − w ki , j (cid:17)(cid:105) L k + − L k L k . (3.4)14 .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 thosedescribed in Section 2.1.2 for the case of static domains. In summary, definitions (2.6) and (2.7) are modified as T k L( i , j ) : = θ L k , T k R( i , j ) : = θ L k for ( i , j ) ∈ [1 , I − × [0 , I ] , (3.5) T k L(0 , j ) : = , T k R( I , j ) : = j ∈ [0 , I ] , T k D( i , j ) : = θ L k , T k U( i , j ) : = θ L k for ( i , j ) ∈ [0 , I ] × [1 , I − , (3.6) T k D( i , : = , T k U( i , I ) : = i ∈ [0 , I ] . Moreover, definitions (2.8) and (2.9) are modified as J k L( i , j ) : = η (cid:16) u k ( i − , j ) − u k ( i , j ) (cid:17) + u max L k , J k R( i , j ) : = η (cid:16) u k ( i + , j ) − u k ( i , j ) (cid:17) + u max L k for ( i , j ) ∈ [1 , I − × [0 , I ] , (3.7) J k L(0 , j ) : = , J k R( I , j ) : = j ∈ [0 , I ] , J k D( i , j ) : = η (cid:16) u k ( i , j − − u k ( i , j ) (cid:17) + u max L k , J k U( i , j ) : = η (cid:16) u k ( i , j + − u k ( i , j ) (cid:17) + u max L k for ( i , j ) ∈ [0 , I ] × [1 , I − , (3.8) J k D( i , : = , J k U( i , I ) : = i ∈ [0 , I ] . Finally, definitions (2.11)-(2.13) are modified as P b (cid:16) n k i , u k i , L k (cid:17) : = τ α n (cid:16) ψ ( n k i ) (cid:17) + φ u ( u k i ) + (cid:16) g i ( n k i , L k ) (cid:17) − , (3.9) P d (cid:16) n k i , u k i , v k i , L k (cid:17) : = τ (cid:16) α n (cid:16) ψ ( n k i ) (cid:17) − φ u ( u k i ) + β n φ v ( v k i ) (cid:17) + (cid:16) g i ( n k i , L k ) (cid:17) + , (3.10)and P q (cid:16) n k i , u k i , v k i , L k (cid:17) : = − τ (cid:16) α n (cid:12)(cid:12)(cid:12) ψ ( n k i ) (cid:12)(cid:12)(cid:12) φ u ( u k i ) + β n φ v ( v k i ) (cid:17) − (cid:12)(cid:12)(cid:12) g i ( n k i , L k ) (cid:12)(cid:12)(cid:12) . (3.11)Here, the function g i ( n k i , L k ) is defined via (3.3) in the case of uniform domain growth and via (3.4) in the case ofapical growth. The functions ψ , φ u and φ v satisfy assumptions (2.14) and (2.15), and we assume τ and L k to be suchthat that 0 < P h < h ∈ { b , d , q } . Similarly to the case of static domains, letting the time-step τ → χ → n ( t , ˆ x ) ∂ t n − ∇ ˆ x · (cid:18) D n L ∇ ˆ x n − C n L n ∇ ˆ x u (cid:19) = (cid:16) α n ψ ( n ) φ u ( u ) − β n φ v ( v ) (cid:17) n + G (ˆ x , n , L ) , ( t , ˆ x ) ∈ R ∗ + × (0 , d (3.12)subject to zero-flux boundary conditions, with either G (ˆ x , w , L ) ≡ G ( w , L ) : = − d w L d L d t , (3.13)15n the case of uniform domain growth, or G (ˆ x , w , L ) : = ˆ x · ∇ ˆ x w L d L d t , (3.14)in the case of apical growth. Here, D n ∈ R ∗ + defined via (2.16) is the rescaled di ff usion coe ffi cient ( i.e. , the rescaledmotility) of the cells, while C n ∈ R + defined via (2.16) represents the chemotactic sensitivity of cells to the activator.In (3.12), the concentration of the activator u ( t , ˆ x ) and the concentration of the inhibitor v ( t , ˆ x ) are governed by thecontinuum counterpart of the di ff erence equations (3.2) complemented with zero-flux boundary conditions, that is, thefollowing system of PDEs subject to zero-flux boundary conditions ∂ t u − D u L ∆ ˆ x u = P ( u , v ) + G (ˆ x , u , L ) ,∂ t v − D v L ∆ ˆ x v = Q ( u , v ) + G (ˆ x , v , L ) , ( t , ˆ x ) ∈ R ∗ + × (0 , d (3.15)which can be formally obtained by letting τ → χ → In this section, we carry out a systematic quantitative comparison between the results of numerical simulations ofthe hybrid model presented in Section 3.1 and numerical solutions of the corresponding continuum model given inSection 3.2, both in one and in two spatial dimensions. All simulations are performed in M atlab and the final time ofsimulations is chosen such that the essential features of the pattern formation process are evident.
We define the functions P , Q , ψ , φ u and φ v as in the case of static domains. In more detail, P and Q are definedvia (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 statedistributions obtained in the case of static domains with (cid:96) : =
1, and we assume the domain to slowly grow linearlyover time, that is, L ( t ) : = + . t . (3.16)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 inAppendix C. The plots in the top rows of Figures 7 and 9 and in the Supplementary Figure D2summarise the dynamics of the continuum concentrations of morphogens u ( t , ˆ x ) and v ( t , ˆ x ) obtained by solvingnumerically the system of PDEs (3.15) subject to zero-flux boundary conditions and with G (ˆ x , u , L ) and G (ˆ x , v , L )defined via (3.13), while the plots in the top rows of Figures 11 and 13 and in the Supplementary Figure D3 referto the case where G (ˆ x , u , L ) and G (ˆ x , v , L ) are defined via (3.14). Identical results hold for the discrete morphogenconcentrations u k i and v k i obtained by solving the system of di ff erence equations (3.2) (results not shown). These plotsdemonstrate 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 thecase 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 8 summarise the dynamics ofthe cell density in the case where there is no chemotaxis, chemically-controlled cell proliferation occurs – i.e. , when η = 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 ) are16igure 7: Results of numerical simulations on a one-dimensional uniformly growing domain in the presence ofchemically-controlled cell proliferation . Top row.
Plots of the concentrations of morphogens at four consecutivetime instants. The green lines highlight the concentration of activator u ( t , ˆ x ) and the red lines highlight theconcentration of inhibitor v ( t , ˆ x ) obtained by solving numerically the system of PDEs (3.15) for d = Bottom row.
Comparison betweenthe discrete cell density n ki obtained by averaging the results of computational simulations of the IB model (solid darkblue lines) and the continuum cell density n ( t , ˆ x ) obtained by solving numerically the PDE (3.12) for d = η = C n =
0, and the functions φ u and φ v are defined via (2.21). We additionally set theinitial cell density n i = for all i . The results from the IB model correspond to the average over five realisationsof the underlying branching random walk, with the results from each realisation plotted in pale blue to demonstratethe robustness of the results obtained. A complete description of the set-up of numerical simulations is given inAppendix C.defined via (3.3) and (3.13), respectively. On the other hand, the plots in the bottom row of Figure 11 and the plotsin 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.4) and (3.14). Theseplots indicate that, when the spatial domain grows over time, spots of high cell density stretch either throughout thedomain (uniform growth) or at the growing edge (apical growth) before splitting. This process causes cell patterns torescale and repeat across the domain at a smaller scale. These plots demonstrate also that there is a good quantitativematch between the discrete cell density n k i given by (2.1), with N k i obtained through computational simulations of theIB model, and the continuum cell density n ( t , ˆ x ) obtained by solving numerically the PDE (3.12) subject to zero-fluxboundary conditions and complemented with either (3.13) or (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 andchemotactic movement of the cells up the concentration gradient of the activator occurs – i.e. , when the functions φ u and φ v are defined via (2.22), η > C n > Conclusions
We have developed a hybrid discrete-continuum modelling framework that can be used to describethe formation of cellular patterns, specifically focusing on the Turing mechanism as the driving force behind thepatterns. We used reaction-di ff usion systems to describe the evolution of morphogens, which dictate the action ofcells, while cell dynamics were described by stochastic IB models. We formally derived the deterministic continuum17igure 8: Results of numerical simulations on a two-dimensional uniformly growing domain in the presenceof chemically-controlled cell proliferation . Comparison between the discrete cell density n k i obtained by averagingthe results of computational simulations of the IB model (top row) and the continuum cell density n ( t , ˆ x ) obtained bysolving numerically the PDE (3.12) for d = η = C n =
0, and the functions φ u and φ v are defined via (2.21). We additionally set the initial cell density n i = × for all i . The results from the IBmodel correspond to the average over five realisations of the underlying branching random walk. The plots of thecorresponding morphogen concentrations are displayed in the Supplementary Figure D2. A complete description ofthe set-up of numerical simulations is given in Appendix C.18igure 9: Results of numerical simulations on a one-dimensional uniformly growing domain in the presenceof chemotaxis . Top row.
Plots of the concentrations of morphogens at four consecutive time instants. The greenlines highlight the concentration of activator u ( t , ˆ x ) and the red lines highlight the concentration of inhibitor v ( t , ˆ x )obtained by solving numerically the system of PDEs (3.15) for d = Bottom row.
Comparison between the discrete cell density n ki obtainedby averaging the results of computational simulations of the IB model (solid dark blue lines) and the continuum celldensity n ( t , ˆ x ) obtained by solving numerically the PDE (3.12) for d = η > C n >
0, andthe functions φ u and φ v are defined via (2.22). We additionally set the initial cell density n i = for all i . The resultsfrom the IB model correspond to the average over five realisations of the underlying branching random walk, with theresults from each realisation plotted in pale blue to demonstrate the robustness of the results obtained. A completedescription of the set-up of numerical simulations is given in Appendix C.19igure 10: Results of numerical simulations on a two-dimensional uniformly growing domain in the presenceof chemically-controlled cell proliferation . Comparison between the discrete cell density n k i obtained by averagingthe results of computational simulations of the IB model (top row) and the continuum cell density n ( t , ˆ x ) obtained bysolving numerically the PDE (3.12) for d = η > C n >
0, and the functions φ u and φ v are defined via (2.22). We additionally set the initial cell density n i = × for all i . The results from the IBmodel correspond to the average over five realisations of the underlying branching random walk. The plots of thecorresponding morphogen concentrations are displayed in the Supplementary Figure D2. A complete description ofthe set-up of numerical simulations is given in Appendix C.20igure 11: Results of numerical simulations on a one-dimensional apically growing domain in the presence ofchemically-controlled cell proliferation . Top row.
Plots of the concentrations of morphogens at four consecutivetime instants. The green lines highlight the concentration of activator u ( t , ˆ x ) and the red lines highlight theconcentration of inhibitor v ( t , ˆ x ) obtained by solving numerically the system of PDEs (3.15) for d = Bottom row.
Comparison between thediscrete cell density n ki obtained by averaging the results of computational simulations of the IB model (solid darkblue lines) and the continuum cell density n ( t , ˆ x ) obtained by solving numerically the PDE (3.12) for d = η = C n =
0, and the functions φ u and φ v are defined via (2.21). We additionally set theinitial cell density n i = for all i . The results from the IB model correspond to the average over five realisationsof the underlying branching random walk, with the results from each realisation plotted in pale blue to demonstratethe robustness of the results obtained. A complete description of the set-up of numerical simulations is given inAppendix C. 21igure 12: Results of numerical simulations on a two-dimensional apically growing domain in the presence ofchemically-controlled cell proliferation . Comparison between the discrete cell density n k i obtained by averagingthe results of computational simulations of the IB model (top row) and the continuum cell density n ( t , ˆ x ) obtained bysolving numerically the PDE (3.12) for d = η = C n =
0, and the functions φ u and φ v are defined via (2.21). We additionally set the initial cell density n i = × for all i . The results from the IBmodel correspond to the average over five realisations of the underlying branching random walk. The plots of thecorresponding morphogen concentrations are displayed in the Supplementary Figure D3. A complete description ofthe set-up of numerical simulations is given in Appendix C.22igure 13: Results of numerical simulations on a one-dimensional apically growing domain in the presence ofchemotaxis . Top row.
Plots of the concentrations of morphogens at four consecutive time instants. The green lineshighlight the concentration of activator u ( t , ˆ x ) and the red lines highlight the concentration of inhibitor v ( t , ˆ x ) obtainedby solving numerically the system of PDEs (3.15) for d = Bottom row.
Comparison between the discrete cell density n ki obtained by averagingthe results of computational simulations of the IB model (solid dark blue lines) and the continuum cell density n ( t , ˆ x )obtained by solving numerically the PDE (3.12) for d = η > C n >
0, and the functions φ u and φ v are defined via (2.22). We additionally set the initial cell density n i = for all i . The results from the IBmodel correspond to the average over five realisations of the underlying branching random walk, with the results fromeach realisation plotted in pale blue to demonstrate the robustness of the results obtained. A complete description ofthe set-up of numerical simulations is given in Appendix C.23igure 14: Results of numerical simulations on a two-dimensional apically growing domain in the presence ofchemically-controlled cell proliferation . Comparison between the discrete cell density n k i obtained by averagingthe results of computational simulations of the IB model (top row) and the continuum cell density n ( t , ˆ x ) obtained bysolving numerically the PDE (3.12) for d = η > C n >
0, and the functions φ u and φ v are defined via (2.22). We additionally set the initial cell density n i = × for all i . The results from the IBmodel correspond to the average over five realisations of the underlying branching random walk. The plots of thecorresponding morphogen concentrations are displayed in the Supplementary Figure D3. A complete description ofthe set-up of numerical simulations is given in Appendix C.24ounterparts of the IB models, which were formulated as PDEs for the cell density, and compared the two modellingapproaches through numerical simulations both in the case of stationary spatial domains and in the case of two typesof growing domains, corresponding to uniform and apical growth. Numerical simulations demonstrated that in thecase of su ffi ciently large cell numbers there was an excellent quantitative match between the spatial patterns producedby 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 di ff erences between the spatial patternsproduced by the two modelling approaches could emerge in the regime of su ffi ciently low cell numbers. In fact, lowercell numbers correlated with both lower regularity of the cell density and demographic stochasticity, which may causea reduction in the quality of the approximations employed in the formal derivation of the deterministic continuummodel from the stochastic IB model. Hence, having both types of models available allows one to use IB models in theregime of low cells numbers – i.e. , when stochastic e ff ects associated with small cell population levels, which cannotbe captured by PDE models, are particularly relevant – and then turn to their less computationally expensive PDEcounterparts when large cell numbers need to be considered – i.e. , when stochastic e ff ects associated with small cellpopulation levels are negligible. Research perspectives
There are a number of additional elements that would be relevant to incorporate into themodelling framework presented here in order to further broaden its spectrum of applications.For instance, as was recognised by Turing himself, exogenous di ff using chemicals are not the only vehicle ofcoordination between cells. In particular, it is known that long range cell-cell interactions can be mediated by signalproteins produced by the cells themselves and also by mechanical forces between cells and components of the cellularmicroenvironment. For example, vascular endothelial growth factor signalling has been shown to control neuralcrest cell migration [53, 54, 55], and mechanical interactions between cells and the extra cellular matrix can controlcell aggregation [62]. Moreover, cellular patterning leading to the emergence of spatial structures often requires theinterplay between non-di ff usible species, transcription factors and cell signalling – viz. the process underlying digitformation in tetrapods [80]. In this regard, it would be interesting to extend the modelling framework by allowingthe cells to consume and / or produce chemicals required for successful coordination of their actions [88], and byincorporating more complex cellular processes such as anoikis [25, 24] and cell deformation [18, 66]. In the situationwhere local production and / or consumption of the chemicals by the cells occurs, for particular cases such as thoseconsidered in [5], we would still expect it to be possible to derive an e ff ective deterministic continuum limit of the IBmodel for the dynamics of the cells through formal procedures analogous to the one used here. However, there couldalso be cases in which PDE models derived using similar formal procedures might not be able to faithfully reproducethe dynamics of the branching random walk underlying the IB model, due to the interplay between stochastic e ff ectsand nonlinear dynamical interactions between the cells and the chemicals.To date, only few biological systems have been demonstrated to satisfy the necessary conditions required forthe formation of Turing pre-patterns via reaction-di ff usion systems. Since mathematical models formulated as scalarintegro-di ff erential 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 [37, 39], it would also beinteresting to explore possible ways of integrating such alternative modelling strategies into our framework. Acknowledgments (All sources of funding of the study must be disclosed)
MAJC gratefully acknowledges support of EPSRC Grant No. EP / N014642 / Conflict of interest
All authors declare no conflicts of interest in this paper. 25 eferences [1] W Alt. Biased random walk models for chemotaxis and related di ff usion approximations. J. Math. Biol. ,9(2):147–177, 1980.[2] P Arcuri and J D Murray. Pattern sensitivity to boundary and initial conditions in reaction-di ff usion models. J.Math. Biol. , 24(2):141–165, 1986.[3] Ruth E Baker, Andrew Parker, and Matthew J Simpson. A free boundary model of epithelial dynamics.
J. Theor.Biol. , 481:61–74, 2019.[4] Benjamin J Binder and Kerry A Landman. Exclusion processes on a growing domain.
J. Theor. Biol. , 259(3):541–551, 2009.[5] Federica Bubba, Tommaso Lorenzi, and Fiona R Macfarlane. From a discrete model of chemotaxis with volume-filling to a generalized patlak–keller–segel model.
Proc. R. Soc. A , 476(2237):20190871, 2020.[6] Martin Burger, Peter Markowich, and Jan-Frederik Pietschmann. Continuous limit of a crowd motion and herdingmodel: analysis and numerical simulations.
Kinet. Relat. Models , 4(4):1025–1047, 2011.[7] Andreas Buttenschoen, Thomas Hillen, Alf Gerisch, and Kevin J Painter. A space-jump derivation for non-localmodels of cell–cell adhesion and non-local chemotaxis.
J. Math. Biol. , 76(1-2):429–456, 2018.[8] Helen M Byrne and Dirk Drasdo. Individual-based and continuum models of growing cell populations: acomparison.
J. Math. Biol. , 58(4-5):657, 2009.[9] Nicolas Champagnat and Sylvie M´el´eard. Invasion and adaptive evolution for individual-based spatiallystructured populations.
J. Math. Biol. , 55(2):147, 2007.[10] M A J Chaplain, T Lorenzi, and F R Macfarlane. Bridging the gap between individual-based and continuummodels of growing cell populations.
J. Math. Biol. , pages 1–29, 2019.[11] Mark A J Chaplain, Mahadevan Ganesh, and Ivan G Graham. Spatio-temporal pattern formation on sphericalsurfaces: numerical simulation and application to solid tumour growth.
J. Math. Biol. , 42(5):387–423, 2001.[12] R Chaturvedi, C Huang, B Kazmierczak, T Schneider, J A Izaguirre, T Glimm, H G E Hentschel, J A Glazier,S A Newman, and M S Alber. On multiscale approaches to three-dimensional modelling of morphogenesis.
J.R. Soc. Interface , 2(3):237–253, 2005.[13] S Christley, M S Alber, and S A Newman. Patterns of mesenchymal condensation in a multiscale, discretestochastic model.
PLoS Comput. Biol. , 3(4):e76, 2007.[14] T M Cickovski, C Huang, R Chaturvedi, T Glimm, H G E Hentschel, M S Alber, J A Glazier, S A Newman,and J A Izaguirre. A framework for three-dimensional simulation of morphogenesis.
Trans. Comput. Biol.Bioinform. , 2(4):273–288, 2005.[15] E J Crampin, E A Ga ff ney, and P K Maini. Reaction and di ff usion on growing domains: Scenarios for robustpattern formation. Bull. Math. Biol. , 61(6):1093–1120, 1999.[16] E J Crampin, W W Hackborn, and P K Maini. Pattern formation in reaction-di ff usion models with nonuniformdomain growth. Bull. Math. Biol. , 64(4):747–769, 2002.[17] Christophe Deroulers, Marine Aubert, Mathilde Badoual, and Basil Grammaticos. Modeling tumor cellmigration: from microscopic to macroscopic models.
Phys. Rev. E , 79(3):031917, 2009.[18] D Drasdo, S Hoehme, and M Block. On the role of physics in the growth and pattern formation of multi-cellularsystems: What can we learn from individual-cell based models?
J. Stat. Phys. , 128(1-2):287, 2007.[19] Dirk Drasdo. Coarse graining in simulated cell populations.
Adv. Complex Syst. , 8(02n03):319–363, 2005.[20] B Duggan and J Metzcar. Generating Turing-like patterns in an o ff -lattice agent based model: Handout. 2019.2621] Louise Dyson, Philip K Maini, and Ruth E Baker. Macroscopic limits of individual-based models for motile cellpopulations with volume exclusion. Phys. Rev. E , 86(3):031903, 2012.[22] B Ermentrout. Stripes or spots? Nonlinear e ff ects in bifurcation of reaction-di ff usion equations on the square. Proc. R. Soc. Lond. Ser. A: Math. Phys. Sci. , 434(1891):413–417, 1991.[23] Anthony E Fernando, Kerry A Landman, and Matthew J Simpson. Nonlinear di ff usion and exclusion processeswith contact interactions. Phys. Rev. E , 81(1):011903, 2010.[24] J Galle, G Aust, G Schaller, T Beyer, and D Drasdo. Individual cell-based models of the spatial-temporalorganization of multicellular systems—Achievements and limitations.
Cytom. A , 69(7):704–710, 2006.[25] J Galle, M Loe ffl er, and D Drasdo. Modeling the e ff ect of deregulated proliferation and apoptosis on the growthdynamics of epithelial cell populations in vitro. Biophys. J. , 88(1):62–75, 2005.[26] A Gierer and H Meinhardt. A theory of biological pattern formation.
Kybernetik , 12:30–39, 1972.[27] C M Glen, M L Kemp, and E O Voit. Agent-based modeling of morphogenetic systems: Advantages andchallenges.
PLoS Comp. Biol. , 15(3), 2019.[28] Thomas Hillen and Hans G Othmer. The di ff usion limit of transport equations derived from velocity-jumpprocesses. SIAM J. Appl. Math. , 61(3):751–775, 2000.[29] Thomas Hillen and Kevin J Painter. A user’s guide to PDE models for chemotaxis.
J. Math. Biol. , 58(1-2):183,2009.[30] Barry D Hughes.
Random walks and random environments: Random walks , volume 1. Oxford University Press,1995.[31] Masaaki Inoue. Derivation of a porous medium equation from many markovian particles and the propagation ofchaos.
Hiroshima Math. J. , 21(1):85–110, 1991.[32] J A Izaguirre, R Chaturvedi, C Huang, T Cickovski, J Co ffl and, G Thomas, G Forgacs, M Alber, G Hentschel,S A Newman, et al. CompuCell, a multi-model framework for simulation of morphogenesis. Bioinformatics ,20(7):1129–1137, 2004.[33] Stuart T Johnston, Ruth E Baker, DL Sean McElwain, and Matthew J Simpson. Co-operation, competition andcrowding: a discrete framework linking allee kinetics, nonlinear di ff usion, shocks and sharp-fronted travellingwaves. Sci. Rep. , 7:42134, 2017.[34] Stuart T Johnston, Matthew J Simpson, and Ruth E Baker. Mean-field descriptions of collective migration withstrong adhesion.
Phys. Rev. E , 85(5):051922, 2012.[35] D Karig, K M Martini, T Lu, N A DeLateur, N Goldenfeld, and R Weiss. Stochastic Turing patterns in a syntheticbacterial population.
Proc. Nat. Acad. Sci. , 115(26):6572–6577, 2018.[36] M A Kiskowski, M S Alber, G L Thomas, J A Glazier, N B Bronstein, J Pu, and S A Newman. Interplay betweenactivator–inhibitor coupling and cell-matrix adhesion in a cellular automaton model for chondrogenic patterning.
Dev. Biol. , 271(2):372–387, 2004.[37] S Kondo. Turing pattern formation without di ff usion. In Conference on Computability in Europe , pages 416–421.Springer, 2012.[38] S Kondo and R Asai. A reaction-di ff usion wave on the skin of the marine angelfish Pomacanthus. Nature ,376(6543):765–768, 1995.[39] Shigeru Kondo. An updated kernel-based turing model for studying the mechanisms of biological patternformation.
J. Theor. Biol. , 414:120–127, 2017.[40] A L Krause, M A Ellis, and R A Van Gorder. Influence of curvature, growth, and anisotropy on the evolution ofTuring patterns on growing manifolds.
Bull. Math. Biol. , 81(3):759–799, 2019.2741] A L Krause, V Klika, T E Woolley, and E A Ga ff ney. From one pattern into another: Analysis of Turing patternsin heterogeneous domains via WKBJ. J. R. Soc. Interface , 17(162):20190621, 2020.[42] Kerry A Landman and Anthony E Fernando. Myopic random walkers and exclusion processes: Single andmultispecies.
Phys. A Stat. Mech. Appl. , 390(21-22):3742–3753, 2011.[43] S S Lee, E A Ga ff ney, and R E Baker. The dynamics of Turing patterns for morphogen-regulated growingdomains with cellular response delays. Bull. Math. Biol. , 73(11):2527–2551, 2011.[44] G Lolas. Spatio-temporal pattern formation and reaction di ff usion systems. Master’s thesis, University of DundeeScotland, 1999.[45] Tommaso Lorenzi, Philip J Murray, and Mariya Ptashnyk. From individual-based mechanical models ofmulticellular systems to free-boundary problems. Interface Free Bound. , In press, 2019.[46] Pavel M Lushnikov, Nan Chen, and Mark Alber. Macroscopic dynamics of biological cells interacting viachemotaxis and direct contact.
Phys. Rev. E , 78(6):061904, 2008.[47] A Madzvamuse, A J Wathen, and P K Maini. A moving grid finite element method applied to a model biologicalpattern generator.
J Comp. Phys. , 190(2):478–500, 2003.[48] Anotida Madzvamuse, Andy H W Chung, and Chandrasekhar Venkataraman. Stability analysis and simulationsof coupled bulk-surface reaction–di ff usion systems. Proc. R. Soc. A , 471(2175):20140546, 2015.[49] P K Maini and T E Woolley. The Turing model for biological pattern formation. In
The Dynamics of BiologicalSystems , pages 189–204. Springer, 2019.[50] P K Maini, T E Woolley, R E Baker, E A Ga ff ney, and S S Lee. Turing’s model for biological pattern formationand the robustness problem. Interface Focus , 2(4):487–496, 2012.[51] Philip´aK Maini, Kevin´aJ Painter, and Helene´aNguyen Phong Chau. Spatial pattern formation in chemical andbiological systems.
Journal of the Chemical Society, Faraday Transactions , 93(20):3601–3610, 1997.[52] Luciano Marcon and James Sharpe. Turing patterns in development: what about the horse part?
Curr. Opin.Genet. , 22(6):578–584, 2012.[53] R McLennan, L Dyson, K W Prather, J A Morrison, R E Baker, P K Maini, and P M Kulesa. Multiscalemechanisms of cell migration during development: Theory and experiment.
Development , 139(16):2935–2944,2012.[54] R McLennan, L J Schumacher, J A Morrison, J M Teddy, D A Ridenour, A C Box, C L Semerad, H Li,W McDowell, D Kay, et al. Neural crest migration is driven by a few trailblazer cells with a unique molecularsignature narrowly confined to the invasive front.
Development , 142(11):2014–2025, 2015.[55] R McLennan, L J Schumacher, J A Morrison, J M Teddy, D A Ridenour, A C Box, C L Semerad, H Li,W McDowell, D Kay, et al. VEGF signals induce trailblazer cell identity that drives neural crest migration.
Dev. Biol. , 407(1):12–25, 2015.[56] H Meinhardt.
Models of Biological Pattern Formation . Academic Press, London, 1982.[57] H Meinhardt. Tailoring and coupling of reaction-di ff usion systems to obtain reproducible complex patternformation during development of the higher organisms. Appl. Math. Comput. , 32(2-3):103–135, 1989.[58] J Moreira and A Deutsch. Pigment pattern formation in zebrafish during late larval stages: A model based onlocal interactions.
Dev. Dyn. , 232(1):33–42, 2005.[59] Sebastien Motsch and Diane Peurichard. From short-range repulsion to hele-shaw problem in a model of tumorgrowth.
J. Math. Biol. , 76(1-2):205–234, 2018.[60] J D Murray. A pre-pattern formation mechanism for animal coat markings.
J. Theor. Biol. , 88:161–199, 1981.2861] J D Murray.
Mathematical Biology: I. An Introduction , volume 17. Springer Science & Business Media, 2007.[62] J D Murray, G F Oster, and A K Harris. A mechanical model for mesenchymal morphogenesis.
J. Math. Biol. ,17(1):125–129, 1983.[63] Philip J Murray, Carina M Edwards, Marcus J Tindall, and Philip K Maini. From a discrete to a continuum modelof cell dynamics in one dimension.
Phys. Rev. E , 80(3):031912, 2009.[64] Philip J Murray, Carina M Edwards, Marcus J Tindall, and Philip K Maini. Classifying general nonlinear forcelaws in cell-based models via the continuum limit.
Phys. Rev. E , 85(2):021921, 2012.[65] M R Myerscough, P K Maini, and K J Painter. Pattern formation in a generalized chemotactic model.
Bull. Math.Biol. , 60(1):1–26, 1998.[66] M P Neilson, D M Veltman, P J M van Haastert, S D Webb, J A Mackenzie, and R H Insall. Chemotaxis: Afeedback-based computational model robustly predicts multiple aspects of real cell behaviour.
PLoS Biol. , 9(5),2011.[67] Karl Oelschl¨ager. On the derivation of reaction-di ff usion equations as limit dynamics of systems of moderatelyinteracting stochastic processes. Probab. Theory Related Fields , 82(4):565–586, 1989.[68] Karl Oelschl¨ager. Large systems of interacting particles and the porous medium equation.
J. Di ff . Eq. , 88(2):294–346, 1990.[69] S Okuda, T Miura, Y Inoue, T Adachi, and M Eiraku. Combining Turing and 3D vertex models reproducesautonomous multicellular morphogenesis with undulation, tubulation, and branching. Sci. Rep. , 8(1):1–15, 2018.[70] H G Othmer and T Hillen. The di ff usion limit of transport equations II: Chemotaxis equations. SIAM J. Appl.Math. , 62(4):1222–1250, 2002.[71] Hans G Othmer, Steven R Dunbar, and Wolfgang Alt. Models of dispersal in biological systems.
J. Math. Biol. ,26(3):263–298, 1988.[72] Hans G Othmer, Philip K Maini, and James D Murray.
Experimental and theoretical advances in biologicalpattern formation , volume 259. Springer Science & Business Media, 2012.[73] Hans G Othmer, Kevin Painter, David Umulis, and Chuan Xue. The intersection of theory and application inelucidating pattern formation in developmental biology.
Math. Model. Nat. Phenom. , 4(4):3–82, 2009.[74] K J Painter, P K Maini, and H G Othmer. Stripe formation in juvenile Pomacanthus explained by a generalizedTuring mechanism with chemotaxis.
Proc. Nat. Acad. Sci. , 96(10):5549–5554, 1999.[75] Kevin J Painter and T Hillen. Volume-filling and quorum-sensing in models for chemosensitive movement.
Can.Appl. Math. Quart. , 10(4), 2002.[76] Kevin J Painter and Jonathan A Sherratt. Modelling the movement of interacting cell populations.
J. Theor. Biol. ,225(3):327–339, 2003.[77] Catherine J Penington, Barry D Hughes, and Kerry A Landman. Building macroscale models from microscaleprobabilistic models: a general probabilistic approach for nonlinear di ff usion and multispecies phenomena. Phys.Rev. E , 84(4):041120, 2011.[78] Catherine J Penington, Barry D Hughes, and Kerry A Landman. Interacting motile agents: Taking a mean-fieldapproach beyond monomers and nearest-neighbor steps.
Phys. Rev. E , 89(3):032714, 2014.[79] J Schnakenberg. Simple chemical reaction systems with limit cycle behaviour.
J. Theor. Biol. , 81(3):389–400,1979.[80] F Schweisguth and F Corson. Self-organization in pattern formation.
Develop. Cell , 49(5):659–677, 2019.[81] J A Sherratt and J D Murray. Mathematical analysis of a basic model for epidermal wound healing.
J. Math.Biol. , 29(5):389–404, 1991. 2982] H Shoji, Y Iwasa, and S Kondo. Stripes, spots, or reversed spots in two-dimensional Turing systems.
J. Theor.Biol. , 224(3):339–350, 2003.[83] Matthew J Simpson, Kerry A Landman, and Barry D Hughes. Cell invasion with proliferation mechanismsmotivated by time-lapse data.
Phys. A Stat. Mech. Appl. , 389(18):3779–3790, 2010.[84] Matthew J Simpson, Alistair Merrifield, Kerry A Landman, and Barry D Hughes. Simulating invasion withcellular automata: connecting cell-scale and population-scale properties.
Phys. Rev. E , 76(2):021918, 2007.[85] A Stevens. The derivation of chemotaxis equations as limit dynamics of moderately interacting stochastic many-particle systems.
SIAM J. Appl. Math , 61(1):183–212, 2000.[86] Angela Stevens and Hans G Othmer. Aggregation, blowup, and collapse: the abc’s of taxis in reinforced randomwalks.
SIAM J. Appl. Math. , 57(4):1044–1081, 1997.[87] A M Turing. The chemical basis of morphogenesis.
Philos. Trans. R. Soc. Lond. B , 237:37–72, 1952.[88] L Tweedy, D A Knecht, G M Mackay, and R H Insall. Self-generated chemoattractant gradients: attractantdepletion extends the range and robustness of chemotaxis.
PLoS Biol. , 14(3), 2016.[89] A Volkening and B Sandstede. Modelling stripe formation in zebrafish: An agent-based approach.
J. R. Soc.Interface , 12(112):20150812, 2015.
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. Similarmethods can be used in the case where d = i , j ) ∈ [1 , I − × [1 , I − n k + i , j ) = n k ( i , j ) + θ L k (cid:104) n k ( i + , j ) + n k ( i − , j ) + n k ( i , j + + n k ( i , j − − n k ( i , j ) (cid:105) + η u max L k (cid:104)(cid:16) u k ( i , j ) − u k ( i − , j ) (cid:17) + n k ( i − , j ) + (cid:16) u k ( i , j ) − u k ( i + , j ) (cid:17) + n k ( i + , j ) (cid:105) + η u max L k (cid:104)(cid:16) u k ( i , j ) − u k ( i , j − (cid:17) + n k ( i , j − + (cid:16) u k ( i , j ) − u k ( i , j + (cid:17) + n k ( i , j + (cid:105) − η u max L k (cid:104)(cid:16) u k ( i − , j ) − u k ( i , j ) (cid:17) + + (cid:16) u k ( i + , j ) − u k ( i , j ) (cid:17) + (cid:105) n k ( i , j ) − η u max L k (cid:104)(cid:16) u k ( i , j − − u k ( i , j ) (cid:17) + + (cid:16) u k ( i , j + − u k ( i , j ) (cid:17) + (cid:105) n k ( i , j ) + τ (cid:16) α n ψ ( n k ( i , j ) ) φ u ( u k ( i , j ) ) − β n φ v ( v k ( i , j ) ) (cid:17) n k ( i , j ) − g ( i , j ) ( n k ( i , j ) , L k ) . (A.1)Using the fact that the following relations hold for τ and χ su ffi ciently small t k ≈ t , t k + ≈ t + τ, ˆ x i ≈ ˆ x , ˆ x i ± ≈ ˆ x ± χ, ˆ y j ≈ ˆ y , ˆ y j ± ≈ ˆ y ± χ n k ( i , j ) ≈ n ( t , ˆ x , ˆ y ) , n k + i , j ) ≈ n ( t + τ, ˆ x , ˆ y ) , n k ( i ± , j ) ≈ n ( t , ˆ x ± χ, ˆ y ) , n k ( i , j ± ≈ n ( t , ˆ x , ˆ y ± χ ) , u k ( i , j ) ≈ u ( t , ˆ x , ˆ y ) , u k + i , j ) ≈ u ( t + τ, ˆ x , ˆ y ) , u k ( i ± , j ) ≈ u ( t , ˆ x ± χ, ˆ y ) , u k ( i , j ± ≈ u ( t , ˆ x , ˆ y ± χ ) , v k ( i , j ) ≈ v ( t , ˆ x , ˆ y ) , v k + i , j ) ≈ v ( t + τ, ˆ x , ˆ y ) , v k ( i ± , j ) ≈ v ( t , ˆ x ± χ, ˆ y ) , v k ( i , j ± ≈ v ( t , ˆ x , ˆ y ± χ ) , L k ≈ L ( t ) , L k + ≈ L ( t + τ ) , n ( t + τ, ˆ x , ˆ y ) = n + θ L (cid:104) n ( t , ˆ x + χ, ˆ y ) + n ( t , ˆ x − χ, ˆ y ) + n ( t , ˆ x , ˆ y + χ ) + n ( t , ˆ x , ˆ y − χ ) − n (cid:105) + η u max L (cid:104)(cid:16) u − u ( t , ˆ x − χ, ˆ y ) (cid:17) + n ( t , ˆ x − χ, ˆ y ) + (cid:16) u − u ( t , ˆ x + χ, ˆ y ) (cid:17) + n ( t , ˆ x + χ, ˆ y ) (cid:105) + η u max L (cid:104)(cid:16) u − u ( t , ˆ x , ˆ y − χ ) (cid:17) + n ( t , ˆ x , ˆ y − χ ) + (cid:16) u − u ( t , ˆ x , ˆ y + χ ) (cid:17) + n ( t , ˆ x , ˆ y + χ ) (cid:105) − η u max L (cid:104)(cid:16) u ( t , ˆ x − χ, ˆ y ) − u (cid:17) + + (cid:16) u ( t , ˆ x + χ, ˆ y ) − u (cid:17) + (cid:105) n − η u max L (cid:104)(cid:16) u ( t , ˆ x , ˆ y − χ ) − u (cid:17) + + (cid:16) u ( t , ˆ x , ˆ y + χ ) − u (cid:17) + (cid:105) n + τ (cid:16) α n ψ ( n ) φ u ( u ) − β n φ v ( v ) (cid:17) n − Γ ( ˆ x , ˆ y , n , L ) , (A.2)with Γ ( ˆ x , ˆ y , n , L ) : = n L ( t + τ ) − L ( t ) L ( t ) , if g ( i , j ) ( n k ( i , j ) ) is defined via (3.3), (cid:34) ˆ x χ (cid:16) n ( t , ˆ x + χ, ˆ y ) − n (cid:17) + ˆ y χ (cid:16) n ( t , ˆ x , ˆ y + χ ) − n (cid:17)(cid:35) L ( t + τ ) − L ( t ) L ( t ) , if g ( i , j ) ( n k ( i , j ) ) is defined via (3.4), where n ≡ n ( t , ˆ x , ˆ y ), u ≡ u ( t , ˆ x , ˆ y ), v ≡ v ( t , ˆ x , ˆ y ) and L ≡ L ( t ). Dividing both sides of (A.2) by τ gives n ( t + τ, ˆ x , ˆ y ) − n τ = θ L τ (cid:104) n ( t , ˆ x + χ, ˆ y ) + n ( t , ˆ x − χ, ˆ y ) + n ( t , ˆ x , ˆ y + χ ) + n ( t , ˆ x , ˆ y − χ ) − n (cid:105) + η u max L τ (cid:104)(cid:16) u − u ( t , ˆ x − χ, ˆ y ) (cid:17) + n ( t , ˆ x − χ, ˆ y ) + (cid:16) u − u ( t , ˆ x + χ, ˆ y ) (cid:17) + n ( t , ˆ x + χ, ˆ y ) (cid:105) + η u max L τ (cid:104)(cid:16) u − u ( t , ˆ x , ˆ y − χ ) (cid:17) + n ( t , ˆ x , ˆ y − χ ) + (cid:16) u − u ( t , ˆ x , ˆ y + χ ) (cid:17) + n ( t , ˆ x , ˆ y + χ ) (cid:105) − η u max L τ (cid:104)(cid:16) u ( t , ˆ x − χ, ˆ y ) − u (cid:17) + + (cid:16) u ( t , ˆ x + χ, ˆ y ) − u (cid:17) + (cid:105) n − η u max L τ (cid:104)(cid:16) u ( t , ˆ x , ˆ y − χ ) − u (cid:17) + + (cid:16) u ( t , ˆ x , ˆ y + χ ) − u (cid:17) + (cid:105) n + (cid:16) α n ψ ( n ) φ u ( u ) − β n φ v ( v ) (cid:17) n − τ Γ ( ˆ x , ˆ y , n , L ) . (A.3)If n ( t , ˆ x , ˆ y ) is a twice continuously di ff erentiable function of ˆ x and ˆ y and a continuously di ff erentiable function of t , u ( t , ˆ x , ˆ y ) is a twice continuously di ff erentiable function of ˆ x and ˆ y , and the function L ( t ) is continuously di ff erentiable,for χ and τ su ffi ciently small we can use the Taylor expansions n ( t , ˆ x ± χ, ˆ y ) = n ± χ ∂ n ∂ ˆ x + χ ∂ n ∂ ˆ x + O ( χ ) , n ( t , ˆ x , ˆ y ± χ ) = n ± χ ∂ n ∂ ˆ y + χ ∂ n ∂ ˆ y + O ( χ ) , n ( t + τ, ˆ x , ˆ y ) = n + τ ∂ n ∂ t + O ( τ ) , u ( t , ˆ x ± χ, ˆ y ) = u ± χ ∂ u ∂ ˆ x + χ ∂ u ∂ ˆ x + O ( χ ) , u ( t , ˆ x , ˆ y ± χ ) = u ± χ ∂ u ∂ ˆ y + χ ∂ u ∂ ˆ y + O ( χ ) , L ( t + τ ) = L + τ d L d t + O ( τ ) . Substituting into (A.3), using the elementary property ( a ) + − ( − a ) + = a for a ∈ R and letting τ → χ → ∂ n ∂ t = D n L (cid:32) ∂ n ∂ ˆ x + ∂ n ∂ ˆ y (cid:33) + C n L (cid:34)(cid:32) ∂ u ∂ ˆ x + ∂ u ∂ ˆ y (cid:33) n − (cid:32) ∂ u ∂ ˆ x ∂ n ∂ ˆ x + ∂ u ∂ ˆ y ∂ n ∂ ˆ y (cid:33)(cid:35) + (cid:16) α n ψ ( n ) φ u ( u ) − β n φ v ( v ) (cid:17) n − G ( ˆ x , ˆ y , n , L ) , ( t , ˆ x , ˆ y ) ∈ R ∗ + × (0 , × (0 , , (A.4)31here G ( ˆ x , ˆ y , n , L ) is given by (3.13) in the case where g ( i , j ) ( n k ( i , j ) ) is defined via (3.3) and by (3.14) in the case where g ( i , j ) ( n k ( i , j ) ) is defined via (3.4). The PDE (A.4) can be easily rewritten in the form (3.12). Moreover, zero-flux boundaryconditions easily follow from the fact that [ cf. definitions (3.5)-(3.8)] T k L(0 , j ) : = , T k R( I , j ) : = , J k L(0 , j ) : = , J k R( I , j ) : = j ∈ [0 , I ]and T k D( i , : = , T k U( i , I ) : = , J k D( i , : = , J k U( i , I ) : = i ∈ [0 , I ] . Remark.
The derivation of the continuum limit for the static domain case can be carried out in a similar way byassuming L k to be constant, which implies that g i ≡ and results in G ≡ . B Set-up of numerical simulations on static domains
We let x ∈ [0 , y ∈ [0 ,
1] and χ : = .
005 ( i.e. , I = τ : = × − . Dynamics of the morphogens
For the dynamics of the morphogens, we consider the parameter setting used in [50],that is, D u : = × − , D v : = × − , α u : = . , β : = , γ : = , α v : = . . (B.1)Moreover, we assume the initial distributions to be small perturbations of the homogeneous steady state ( u ∗ , v ∗ ) ≡ (1 , . u i = u ∗ − ρ + ρ R and v i = v ∗ − ρ + ρ R where ρ : = .
001 and R is either a vector for d = d = , atlab function rand . Thesechoices of the initial distributions of morphogens are such that u ∗ − ρ ≤ u i ≤ u ∗ + ρ and v ∗ − ρ ≤ v i ≤ v ∗ + ρ for all i , that is, the parameter ρ determines the level of perturbation from the homogeneous steady state. Since the di ff erenceequations (2.2) governing the dynamics of the morphogens are independent from the dynamics of the cells, suchequations are solved first for all time-steps and the solutions obtained are then used to evaluate both the probabilitiesof 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 k , i u k i . Computational implementation of the rules underlying the dynamics of the cells
At each time-step, each cellundergoes a three-phase process: Phase 1) undirected, random movement according to the probabilities definedvia (2.6) and (2.7); Phase 2) chemotaxis according to the probabilities defined via (2.8) and (2.9); Phase 3) divisionand death according to the probabilities defined via (2.11)-(2.13). For each cell, during each phase, a random numberis drawn from the standard uniform distribution on the interval (0 ,
1) using the built-in M atlab function rand . It isthen 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 n i ≡ when d = n i ≡ × when d = . 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 = θ : = . , η : = , α n : = , β n : = , n max : = × . and the following ones when d = θ : = . , η : = , α n : = , β n : = . , n max : = × . d = n i ≡ × and n max : = . × and when d = n i ≡ × and n max : = × . In the case where cells undergo chemotaxis and cell proliferation is not chemically-controlled, unless statedotherwise, we use the following parameter values when d = θ : = . , η : = , α n : = . , β n : = . , n max : = × . and the following ones when d = θ : = . , η = , α n : = . , β n : = . , n max : = × . Numerical solutions of the corresponding continuum models
Numerical solutions of the PDE (2.17) and thesystem of PDEs (2.18) subject to zero-flux boundary conditions are computed through standard finite-di ff erenceschemes using initial conditions and parameter values that are compatible with those used for the IB model andthe system of di ff erence equations (2.2). In particular, the values of the parameters D n and C n in the PDE (2.17) aredefined via (2.23). C Set-up of numerical simulations on growing domains
We let x ∈ [0 , y ∈ [0 ,
1] and χ : = .
005 ( i.e. , I = τ : = × − and we define L according to (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 (B.1).Moreover, we define the initial distributions as the numerical equilibrium distributions obtained in the case of staticdomains. Similarly to the case of static domains, since the di ff erence equations (3.2) governing the dynamics of themorphogens are independent from the dynamics of the cells, such equations are solved first for all time-steps and thesolutions obtained are then used to evaluate both the probabilities of cell movement (3.5)-(3.8) and the probabilities ofcell division and death (3.9)-(3.11). The parameter u max in (3.7) and (3.8) 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 staticdomains, at each time-step, each cell undergoes a three-phase process: Phase 1) undirected, random movementaccording to the probabilities defined via (3.5) and (3.6); Phase 2) chemotaxis according to the probabilities definedvia (3.7) and (3.8); Phase 3) division and death according to the probabilities defined via (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 thebuilt-in M atlab function rand . It is then evaluated whether this number is lower than the probability of the eventoccurring 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 usedin the static domain case.
Numerical solutions of the corresponding continuum models
Numerical solutions of the PDE (3.12) and thesystem of PDEs (3.15) subject to zero-flux boundary conditions are computed thorugh standard finite-di ff erenceschemes using initial conditions and parameter values that are compatible with those used for the IB model andthe system of di ff erence equations (3.2). In particular, the values of the parameters D n and C n in the PDE (3.12) aredefined via (2.23). D Supplementary figures
Dynamics of the morphogens on a two-dimensional static domain.
Plots of the concentration ofactivator u ( t , x ) (top row) and the concentration of inhibitor v ( t , x ) (bottom row) at four consecutive time instants,obtained by solving numerically the system of PDEs (2.18) for d = Dynamics of the morphogens on a two-dimensional uniformly growing domain.
Plots of theconcentration of activator u ( t , ˆ x ) (top row) and the concentration of inhibitor v ( t , ˆ x ) (bottom row) at four consecutivetime instants, obtained by solving numerically the system of PDEs (3.15) for d =
2, subject to zero-flux boundaryconditions, complemented with (2.19), (3.13) and (3.16). A complete description of the set-up of numericalsimulations is given in Appendix B. 34igure D3:
Dynamics of the morphogens on a two-dimensional apically growing domain.
Plots of theconcentration of activator u ( t , x ) (top row) and the concentration of inhibitor v ( t , x ) (bottom row) at four consecutivetime instants, obtained by solving numerically the system of PDEs (3.15) for d ==