Unbiased on-lattice domain growth

TL;DR

This paper introduces an unbiased method for modeling on-lattice domain growth in biological systems, effectively preventing boundary particle buildup and improving accuracy over previous methods, especially when diffusion is low.

Contribution

A new unbiased on-lattice domain growth method that corrects boundary particle buildup and extends the applicability of models to low diffusion scenarios.

Findings

The new method prevents unphysical boundary particle accumulation.

It remains accurate even when diffusion is low.

The original method is still feasible when diffusion dominates domain growth.

Abstract

Domain growth is a key process in many areas of biology, including embryonic development, the growth of tissue, and limb regeneration. As a result, mechanisms for incorporating it into traditional models for cell movement, interaction, and proliferation are of great importance. A previously well-used method in order to incorporate domain growth into on-lattice reaction-diffusion models causes a build up of particles on the boundaries of the domain, which is particularly evident when diffusion is low in comparison to the rate of domain growth. Here, we present a new method which addresses this unphysical build up of particles at the boundaries, and demonstrate that it is accurate even for scenarios in which the previous method fails. Further, we discuss for which parameter regimes it is feasible to continue using the original method due to diffusion dominating the domain growth mechanism.

Equations

$${\alpha }_0=\sum _{i=1}^k\left({\alpha }_i^L+{\alpha }_i^R+{\alpha }_i^G\right).$$

$$m_{j}=\left\{\begin{array}[]{ll}r_{j},&\text{if }j\in\{1,...,i-1\},\\ b,&\text{if }j=i,\\ r_{i}-b,&\text{if }j=i+1,\\ r_{j-1},&\text{if }j\in\{i+2,...,k\}.\end{array}\right.$$

$$\frac{\text{d}p}{\text{d}t}(\mathbf{m},k,t)=\rho \sum _{j=1}^{k-1}\pi (m_j,m_{j+1}|m_j+m_{j+1})p(G_j\mathbf{m},k-1,t)-\rho kp(\mathbf{m},k,t).$$

$$G_j:(m_1,...,m_j,m_{j+1},...,m_k{)}^T\mapsto (m_1,...,m_j+m_{j+1},...,m_k{)}^T.$$

$$\pi (x,y|z)=(\genfrac{}{}{0ex}{}{z}{x}){\left(\frac{1}{2}\right)}^z,$$

$$M_i^k(t)=\sum _{\mathbf{m}\in {\mathbb{N}}_0^k}{m}_{i}p(\mathbf{m},k,t),$$

$$\frac{\text{d}{M}_i^k}{\text{d}t}(t)=\rho \left[\left(i-\frac{3}{2}\right)M_{i-1}^{k-1}(t)+\left(k-\frac{1}{2}-i\right)M_i^{k-1}(t)-k{M}_i^k(t)\right],$$

$$u_i^{k+1}=\underset{\begin{array}{c}\text{Compartment}\\ \text{to right of }\\ i\text{ splits}\end{array}}{\underbrace{u_i^k\text{mathbbm}{1}_{[I^k>i]}}}+\underset{\begin{array}{c}\text{Compartment }\\ i\text{ splits}\end{array}}{\underbrace{q{u}_i^k\text{mathbbm}{1}_{[I^k=i]}}}+\underset{\begin{array}{c}\text{Compartment }\\ i-1\text{ splits}\end{array}}{\underbrace{(1-q)u_{i-1}^k\text{mathbbm}{1}_{[I^k=i-1]}}}+\underset{\begin{array}{c}\text{Compartment}\\ \text{to left of }\\ i-1\text{ splits}\end{array}}{\underbrace{u_{i-1}^k\text{mathbbm}{1}_{[I^k<i-1]}}},i\in \{2,...,k\},$$

$$u_1^{k+1}=u_1^k\text{mathbbm}{1}_{[I^k>1]}+q{u}_1^k\text{mathbbm}{1}_{[I^k=1]}.$$

$$M_1^{k+1}$$

$$=M_1^k\left(\frac{k-1}{k}+\frac{1/2}{k}\right),$$

$$=M_1^k\left(1-\frac{1/2}{k}\right).$$

$$M_1^{k+1}=M_1^1\prod _{n=1}^k\left(1-\frac{1/2}{n}\right)=\prod _{n=1}^k\left(1-\frac{1/2}{n}\right),$$

$$\prod _{n=a}^b\left(1+\frac{c}{n}\right)\approx \frac{\Gamma (a)}{\Gamma (a+c)e^c}{b}^c,$$

$$M_1^{k+1}=\prod _{n=1}^k\left(1-\frac{1/2}{n}\right)\approx \frac{\Gamma (1)}{\Gamma (1/2)e^{-1/2}}{k}^{-1/2}=\frac{1}{\sqrt{\pi }{e}^{-1/2}}{k}^{-1/2}.$$

$$u_2^{k+1}=u_2^k\text{mathbbm}{1}_{[I^k>2]}+\frac{1}{2}{u}_2^k\text{mathbbm}{1}_{[I^k=2]}+\frac{1}{2}{u}_1^k\text{mathbbm}{1}_{[I^k=1]}.$$

$$M_2^{k+1}$$

$$⋮$$

$$=M_2^2\underset{\text{Term 1}}{\underbrace{\prod _{n=2}^k\left(1-\frac{3/2}{n}\right)}}+\underset{\text{Term 2}}{\underbrace{\sum _{n=2}^{k-1}\left(\prod _{j=n+1}^k\left(1-\frac{3/2}{j}\right)\right)M_1^n\left(\frac{1/2}{n}\right)}}+\underset{\text{Term 3}}{\underbrace{M_1^k\left(\frac{1/2}{k}\right)}}.$$

$$\prod _{n=2}^k\left(1-\frac{3/2}{n}\right)\approx \frac{\Gamma (2)}{\Gamma (1/2)e^{-3/2}}{k}^{-3/2}.$$

$$\sum _{n=2}^{k-1}\left(\prod _{j=n+1}^k\left(1-\frac{3/2}{j}\right)\right)M_1^n\left(\frac{1/2}{n}\right)\approx \frac{1}{2}\frac{1}{\sqrt{\pi }{e}^{-1/2}}{k}^{-1/2}-\frac{1}{\sqrt{\pi }{e}^{-1/2}}{k}^{-3/2}.$$

$$M_1^k\left(\frac{1/2}{k}\right)\approx \frac{1}{2}\frac{1}{\sqrt{\pi }{e}^{-1/2}}{k}^{-3/2}$$

$$M_2^{k+1}\approx M_2^2\frac{\Gamma (2)}{\Gamma (1/2)e^{-3/2}}{k}^{-3/2}+\frac{1}{2}\frac{1}{\sqrt{\pi }{e}^{-1/2}}{k}^{-1/2}-\frac{1}{\sqrt{\pi }{e}^{-1/2}}{k}^{-3/2}+\frac{1}{2}\frac{1}{\sqrt{\pi }{e}^{-1/2}}{k}^{-3/2},$$

$$M_2^{k+1}\approx \frac{1}{2}\frac{1}{\sqrt{\pi }{e}^{-1/2}}{k}^{-1/2}.$$

$$M_i^k\approx c_i{k}^{-1/2},$$

$$c_1$$

$$c_i$$

$${\delta }_i^{k-1}\frac{k}{k-1}h$$

$$=ih-\frac{k}{k-1}ih+\frac{k}{k-1}h$$

$$=\frac{k-i}{k-1}h.$$

$$⟹{\delta }_i^{k-1}$$

$$\mathbf{m}=\sum _{i=1}^{k-1}{b}_i{\mathbf{e}}_i+(r_i-b_i){\mathbf{e}}_{i+1},$$

$${\alpha }_0=\sum _{i=1}^k\left({\alpha }_i^L+{\alpha }_i^R\right)+{\alpha }^G.$$

$$m_{j}=\left\{\begin{array}[]{ll}b_{1},&\text{if }j=1,\\ b_{j}+(r_{j-1}-b_{j-1}),&\text{if }j\in\{2,...,k-1\},\\ r_{k-1}-b_{k-1},&\text{if }j=k.\end{array}\right.$$

Full text

\nobibliography

\newlistofAlgorithmexpList of Algorithms

**Unbiased on-lattice domain growth **

Cameron A. Smith1,∗, Cécile Mailler2, Christian A. Yates1

1**Centre for Mathematical Biology, Department of Mathematical Sciences, University of Bath, Claverton Down, Bath, BA2 7AY, United Kingdom

2Probability Laboratory, Department of Mathematical Sciences, University of Bath, Claverton Down, Bath, BA2 7AY, United Kingdom

∗E-mail: c.smith3@bath.ac.uk**

Key words: domain growth, mesoscopic modelling, reaction-diffusion master equation, morphogen gradient, uniform growth, on-lattice modelling

1 Introduction

Domain growth is an inherent feature of many biological systems, from neural crest cell migration [1, 2] to the growth and shrinkage of tissue [3, 4], and it has been investigated in the context of pattern formation in reaction-diffusion systems [5, 6]. It is therefore important that we have reliable mathematical tools in order to model such systems.

One method of modelling general reaction-diffusion systems is to compartmentalise the spatial domain into a lattice of small regions, in which particles are considered to reside. Typically, particles are permitted to jump between neighbouring compartments (although non-local jumping can also be incorporated [7]), and particles may react with others in their current compartment (similarly, in some methods, particles may be allowed to interact with particles in neighbouring compartments [8]). These events, under the assumption that updates happen in continuous time, are considered Markovian and, as a result, have associated exponentially distributed waiting times. The most commonly used method to simulate such systems is the Gillespie algorithm [9], however there are many others that are also used within the literature (see, for example, the next reaction method [10], the next subvolume method [11, 12], and the sorting direct method [13]).

The word ‘particle’ in this context can refer to a multitude of biological or physical entities; a particle may be a cell when modelling biological systems such as neural crest cell migration in embryonic development, a chemical species when considering pattern formation or different classes of individual in the case of the spread of an epidemic. The discretisation of the domain does not necessarily correspond to any physical attribute of the space being modelled, but is usually used as a mathematical tool through which we are able to model stochasticity in particle positions and domain growth. As an example, consider the formation of a morphogen gradient on a growing domain [14]. In this example, the particles are the morphogen molecules, and the domain may be discretised to arbitrary accuracy, with no physical meaning necessarily being attached to the compartments. This case study will be investigated in more detail in Section 3.1.

One possible option for incorporating domain growth into these on-lattice position-jump processes has been previously suggested by Baker et al. [15]. Growth is achieved by choosing compartments to divide uniformly at random at a given rate. Once chosen, a compartment instantaneously doubles in length and splits down the middle to produce a new compartment (pushing all compartments to the right of the one chosen by one compartment’s width). Particles that resided in the original compartment are then redistributed into the newly created compartments by placing each particle into one or the other with equal probability. This is illustrated schematically in Figure 1, and a method for its implementation is detailed in Algorithm 1. A commonly implemented simulation allows particles to diffuse by jumping between neighbouring boxes while, at the same time, attempting to grow the domain uniformly [15, 3, 16]. When simulating this scenario using the method of Baker et al. [15] and averaging over multiple repeats, it becomes apparent that there is an issue with particles building up at the ends of the domain, an effect that can be expected when diffusion is low in comparison to the domain growth rate. When diffusion is larger, there is enough time for particles to relax to equilibrium (particles have enough time to occupy newly created sites) before the next domain growth event.

Algorithm \theAlgorithm: The original method [15]

• (0a)

Initialise time $t=0$ and set the final time $T>0$. Initialise the number of compartments $k$ and the particle numbers in each compartment $m_{i},\ i\in\{1,...,k\}$. Specify the size of compartment $h$, the jump rate $d=D/h^{2}$ and the growth rate $\rho$.

• (0b)

Calculate the propensity functions $\alpha_{i}^{L}=dm_{i}$, $\alpha_{i}^{R}=dm_{i}$ and $\alpha_{i}^{G}=\rho$ for each compartment $i\in\{1,...,k\}$ for left jumps, right jumps and growth events respectively, with $\alpha_{1}^{L}=0=\alpha_{k}^{R}$. Calculate the sum of the propensity functions:

[TABLE]

• (0c)

Calculate the time until the next event by firstly drawing a uniform random variable between zero and one, $u_{1}\sim\text{Unif}(0,1)$, and setting $\tau=1/\alpha_{0}\log\left(1/u_{1}\right)$. Update the time $t\leftarrow t+\tau$.

• (0d)

Determine which event, amongst all compartments, is next to occur by choosing one at random with probability proportional to the propensity function.

If the event corresponds to a left (resp. right) jump event from compartment $i$, remove a particle from compartment $i$ and add one to compartment $i-1$ (resp. $i+1$). 2. 2.

If the event corresponds to a growth event in compartment $i$: Record the pregrowth particle numbers $\boldsymbol{r}\leftarrow\boldsymbol{m}$, create an extra compartment at the right end of the postgrowth domain (increasing $k$ by 1 by setting $k\leftarrow k+1$), draw a binomial random variable with $r_{i}$ trials and probability of success 1/2, $b\sim\text{Bin}(r_{i},1/2)$, and redistribute particles as follows (for $j\in\{1,...,k\}$):

[TABLE]

• (0e)

If $t<T$ then return to step (\theAlgorithmb), else end.

In order to rectify this problem we have designed a novel mechanism which enables the implementation of unbiased domain growth — one that prevents the accumulation of particles at the boundaries. This method enacts a more continuous approach to domain growth compared to Baker et al. [15], in which compartments are first stretched and then renormalised with concurrent redistribution of their particles in to neighbouring compartments. In general, we will consider growth dynamics in isolation, without considering any other mechanisms such as diffusion and reactions. This is in order to demonstrate the possible issues without confounding or hiding growth induced phenomena with other dynamics. All models will be presented in one dimension, but extensions to higher dimensions will be discussed in Section 4.

This paper will be set out in the following way. In Section 2 we investigate and explain the problems with the commonly employed domain growth mechanism in more detail. In Section 3, we present a novel, but distinct and importantly unbiased domain growth method that we use to address this issue. Finally, we present our conclusions in Section 4.

2 Problems with an existing domain growth mechanism

Within this section, we describe and demonstrate problems with an existing domain growth mechanism at the mesoscale [15], which we will refer to as the “original method”. We then postulate why the phenomenon has not been noticed before and confirm the existence of the problem using three different techniques.

The original method causes an accumulation of particles at the ends of the domain. This build up is caused by the inherent bias in the way in which growth is implemented: compartments are always shifted to the right when one is chosen to split. When a single compartment is chosen to split, most postgrowth compartments will either retain their pregrowth contents (if a pregrowth compartment to its right splits) or will take on the contents of the pregrowth compartment to its left (if the split is to the left of that pregrowth compartment). The notable exceptions to these rules are when postgrowth compartments retain or gain (respectively) only half of the particles (on average) from its or a neighbour’s (respectively) pregrowth compartment. These events only affect two compartments per splitting event, but are nevertheless important.

The only way a postgrowth compartment retains half of its particles (on average) is if it is the compartment chosen to split. Similarly, the only way a postgrowth compartment gains half of the particles (on average) of a neighbouring compartments is if the pregrowth neighbour immediately to the left splits. All the postgrowth compartments, therefore, gain or retain half a compartments worth of particles (on average) in two different ways. The only exceptions are the first and last postgrowth compartments. The first postgrowth compartment can only retain half of its particles (on average) if it is chosen to split. It can never gain half of the particles from another compartment as there are no pregrowth compartments to its left that can split. Similarly, the last postgrowth compartment can never retain half of its particles since it did not exist on the pregrowth domain. Instead it can only gain half of the particles from a neighbouring compartment when the pregrowth compartment in the final position splits.

If particles are initially distributed uniformly, then an unbiased domain growth method will maintain this uniform distribution (on average, and provided no particles enter or leave the domain). For the method of Baker et al. [15] (see Figure 1 and Algorithm 1), splitting events will redistribute the particles into the postgrowth compartments which correspond to the pregrowth compartment chosen to split and its neighbour to the right. For example, in Figure 1, the particle redistribution event affects postgrowth compartments 5 and 6 when pregrowth compartment 5 is chosen to split. If these splitting events affected all compartments equally then the domain growth method would be unbiased and the particle profile would remain uniform as the domain grew. However, since the two end compartments suffer (on average) only half of the particle-reducing splitting events that the non-end compartments suffer, this leads to particles accumulating at the ends of the domains.

The build up of particles at either end of the domain has previously been hidden by the smoothing effect of fast diffusion [15]. If diffusion is large in comparison to the rate of domain growth, particles are able to diffuse away from the high concentration regions at the end of the domain, leading to a near-uniform particle profile (given a uniform initial condition and zero-flux boundary conditions).

We will illustrate the problems with the compartment-based domain growth method of Baker et al. [15] in three different ways. For the first we undertake a stochastic simulation of the original method using the Gillespie algorithm [9] and average particle densities over many repeats. For the second demonstration, we calculate the numerical solution of the mean-field equations that stem from the corresponding domain growth master equation of the stochastic process in the absence of diffusion. Thirdly, we employ an analytical mathematical argument based on local redistribution of particles to assess the particle distribution in the limit of large numbers of boxes.

For the stochastic simulation, we initialise a total of 1000 particles uniformly in the domain $[0,4]$ which grows exponentially in time with rate $\rho$ (we will use exponential growth as our primary test simulation, however other types of growth exhibit similar problems). We set the compartment width to be $0.4$ so that initially we have 10 compartments. The growth rate, $\rho$, is set to be 0.01 (note that all units here are arbitrary in both space and time). We choose compartment splitting times to be deterministic, meaning that we pre-calculate the times at which a compartment should split (on average) and always enact a splitting event at those times, while the Gillespie algorithm handles the stochastic events between these pre-determined times. Alternatively, growth events could be incorporated stochastically as part of the Gillespie algorithm, but for the purposes of demonstrating the problem, and for ease of visualisation, we use deterministic growth in this exposition. The compartment to divide at each predetermined splitting time is chosen uniformly at random from amongst the current compartments. Diffusion is set to be 0, so that the only effect on compartment occupancy is domain growth. Using this domain growth configuration and averaging over 50,000 repeats, we see the clear build up of particles at both ends of the domain (see Figure 2).

To further corroborate these results and explain where this behaviour originates from, we consider the master equation for the splitting algorithm outlined above, which was first derived in [15]. Let $p(\boldsymbol{m},k,t)=\mathbb{P}(\boldsymbol{N}(t)=\boldsymbol{m},K(t)=k)$ be the probability that the state variable (the number of particles in each compartment) at time $t$ is $\boldsymbol{m}=(m_{1},...,m_{k})^{T}$ and the number of compartments at time $t$ is $k$. This quantity evolves according to:

[TABLE]

Here $\rho$ is the splitting rate, which is the rate at which each compartment is chosen to divide, $\pi(x,y|z)$ is a distribution describing the probability that, given there are $z$ particles in a compartment before splitting, there are $x$ and $y$ particles in the two post-split compartments (where $x+y=z$), and $G_{j}:\mathbb{R}^{k}\rightarrow\mathbb{R}^{k-1}$ is an operator that combines the contents of compartments $j$ and $j+1$ (the opposite process to splitting), so that

[TABLE]

Note that the growth considered for the master equation is stochastic in both position and timing.

We can calculate the mean-field equations for the evolution of the mean number of particles in each compartment, under the splitting probability corresponding to a symmetric binomial distribution

[TABLE]

by multiplying both sides of equation (1) by $m_{i}$ for each index $i\in\{1,...,k\}$, and summing over the entire postgrowth state variable. Define the quantity

[TABLE]

which is the mean number of particles in compartment $i$ at time $t$, given that there are $k$ compartments in the system overall. These corresponding mean-field equations are given by:

[TABLE]

which holds for $k\in\mathbb{N},\ j\in\{1,...,k\}$ and $t>0$ [15], noting that $M_{k}^{k-1}\equiv 0$. The full derivation is omitted here, however we direct the interested reader to Baker et al. [15], or Appendix D.1 for a similar calculation. We plot the solutions to equation (2) in order to illustrate the average behaviour of the system under the biased domain growth algorithm, in Figure 3. We show that in this case, we exhibit the same build up of particles at the boundaries that is evident in the stochastic simulation, shown in Figure 2.

Finally, we gain a more quantitative insight into the bias engendered by the domain growth method of Baker et al. [15] by using an analytical approach. In particular, we derive coefficients which describe the densities of particles in each compartment. We let $u_{i}^{k}$ denote the normalised density of particles in compartment $i$ when there are $k$ compartments in total (a random variable), $M_{i}^{k}=\mathbb{E}[u_{i}^{k}]$, and $I^{k}$ is the index of the compartment which splits when there are $k$ compartments in total (chosen uniformly at random for each event). Finally, we define $q$ to be a sample from some symmetric distribution with mean $1/2$. This will denote the random proportion of density that is placed in the left-hand postgrowth compartment when there are $k$ compartments in total, and is independent of the particle density and the compartment chosen to split. We begin with all of the density in the first (and only) compartment, so that $u_{1}^{1}=1$. We will set up recursion relationships between the $u_{i}^{k}$’s for different $i$ and $k$ in order to approximate $u_{i}^{k}$ for large $k$. Specifically, we express $u_{i}^{k+1}$ in terms of $u_{j}^{k}$ for $j\in\{1,...,k\}$:

[TABLE]

with similar expressions for $i=1$ and $k+1$. Here $\mathbbm{1}_{\text{[condition]}}$ is the indicator function, which is unity when the subscripted condition is satisfied and zero otherwise.

Considering the first compartment, relationship (3) stipulates that

[TABLE]

Taking expectations of this expression, noting that the choice of compartment, the density in each compartment and the proportion of the density placed in the left pregrowth compartment are independent (so that $u_{i}^{k}$, $I^{k}$ and $q$ are independent for every $i\in\{1,...,k\}$) we find that:

[TABLE]

Applying this recursive relation $k$ times:

[TABLE]

since $M_{1}^{1}=\mathbb{E}[u_{1}^{1}]=1$. In order to simplify this expression, we can use the following relationship, which is derived in Appendix A:

[TABLE]

where $a,b\in\mathbb{N}$ such that $a<b$ and $b$ is large, and $c\in\mathbb{R}$. Applying approximation (5) to equation (4) for large $k$ gives:

[TABLE]

Now consider the second compartment. As with the first, we write the recursion relation (3):

[TABLE]

Once again, we take expectations and simplify by applying relation (3) recursively (as in equation (4)):

[TABLE]

We consider each of the three terms in equation (7) sequentially. Using equation (5), we can approximate term 1, for large $k$, as

[TABLE]

It can be shown (see Appendix B) that for large $k$, term 2 can be approximated by

[TABLE]

To simplify term 3 we make use of result (6):

[TABLE]

Substituting the resultant expressions (8), (9) and (10) into (7) gives, for large $k$:

[TABLE]

For large $k$, the $O(k^{-1/2})$ terms dominate leaving us with the following approximation:

[TABLE]

Following the same procedure, we can find the approximate expressions for the asymptotic particle density in each of the compartments. In particular, it can be shown that

[TABLE]

where

[TABLE]

To assess the accuracy of this approximation, we undertake a stochastic simulation initialised with a single compartment containing 100 particles. At time $t=100$, under our time-deterministic splitting mechanism each simulation finishes with 20 compartments. We then compare the particle numbers in each compartment, averaged over 10,000 repeats, to equation (12). Since the simulation domain has only finitely many compartments but our mathematical analysis considers an infinite number of compartments, we average $M_{i}^{20}$ and $M_{21-i}^{20}$ when we plot compartment $i$ under the assumption that densities are initialised, and subsequently remain, symmetric. Finally, we scale each of the plots so that they have the same number of particles. Although a quantitative agreement is not expected, since our results hold strictly only on an infinite domain, the results in Figure 4 demonstrate that our mathematical analysis matches the simulation results qualitatively.

The original domain growth method proposed by Baker et al. [15] has been used in many compartment-based studies of domain growth [17, 18, 19, 3]. Despite not having previously been evident, we have been able to demonstrate in three distinct ways, that this domain growth method is biased. The consequence of this bias is that particles tends to accumulate at the extreme ends of the domain. In the next section, we introduce the stretching method, which prevents this build up of particles and gives genuinely uniform, unbiased domain growth.

3 Stretching method

We now introduce the stretching method. This differs from the original method because it is a global method as opposed to a local one. That is, instead of choosing a single compartment to instantaneously grow to twice its length and split, we stretch all compartments by a small amount and redistribute particles amongst all compartments (for a brief discussion of local growth mechanisms, please see Appendix C). We will firstly explain the method, before demonstrating its effectiveness. We do this by showing that it can correctly maintain a uniform particle profile on a uniformly growing domain. We then look at a case study, the formation of a morphogen gradient on an exponentially growing domain, in order to directly compare the original and stretching methods using an example with its roots in biology (see Section 3.1). Finally, in Section 3.2, we investigate the parameter regimes in which the spatial inhomogeneities in the original method are negligible.

We begin by describing the method. Assume the number of compartments before a growth event is $k-1$ for some $k>2$, and define each compartment to be of width $h$ (see Figure 5(a)). We will define the state variable before growth to be $\boldsymbol{r}\in\mathbb{N}_{0}^{k-1}$ and after growth to be $\boldsymbol{m}\in\mathbb{N}_{0}^{k}$ in order to be consistent with the original method. The method proceeds as follows:

When a growth event is chosen to occur, we stretch the domain to be of size $kh$ rather than $(k-1)h$ (see Figure 5(b)). We do this uniformly across the entirety of the domain, so that each compartment on the stretched domain is now of width $kh/(k-1)$. 2. 2.

In the second step we add a compartment to the right-hand end of the pre-stretched domain (see Figure 5(c)). It is on this postgrowth domain that we define the state variable $\boldsymbol{m}$. Note that we now have two domains, each with a different number of compartments, but both of the same length. 3. 3.

For the third step, we compare the two meshes. Note that for every stretched compartment, exactly two of the postgrowth compartments intersect it (see Figure 5(d)). Assuming particles are uniformly distributed across each stretched compartment, we can calculate the proportion of particles, $\delta_{i}^{k-1}$, that should be placed in the left overlapping compartment. If we denote the right-hand end of compartment $i$ in the renormalised domain as $x_{i}=ih$ and use $\tilde{x}_{i}=\left[ki/(k-1)\right]h$ the quantity for the stretched domain, then:

[TABLE] 4. 4.

Finally we calculate the new state $\boldsymbol{m}$ by drawing $k-1$ binomial random variables $b_{i}\sim\text{Bin}(r_{i},\delta_{i}^{k-1})$ and calculating $\boldsymbol{m}$ as:

[TABLE]

where the $\boldsymbol{e}_{i}$ are the standard $k$-dimensional basis vectors.

Algorithm \theAlgorithm: The stretching method

• (0a)

Initialise time $t=0$ and set the final time $T>0$. Initialise the number of compartments, $k$, and the particle numbers in each compartment $m_{i},\ i\in\{1,...,k\}$. Specify the size of compartment $h$, the jump rate $d=D/h^{2}$ and the growth rate $\rho$.

• (0b)

Calculate the propensity functions $\alpha_{i}^{L}=dm_{i}$, $\alpha_{i}^{R}=dm_{i}$ for left and right jumps from each compartment for $i=1,...,k$ with $\alpha_{1}^{L}=0=\alpha_{k}^{R}$, and set the propensity function for a growth event to be $\alpha^{G}=\rho k$. Calculate the sum of the propensity functions:

[TABLE]

• (0c)

Calculate the time until the next event by firstly drawing $u_{1}\sim\text{Unif}(0,1)$ and setting $\tau=1/\alpha_{0}\log\left(1/u_{1}\right)$. Set the time to be the next event that occurs $t\leftarrow t+\tau$.

• (0d)

Determine which event is next to occur by choosing one at random with probability proportional to the propensity function.

If the event corresponds to a left (resp. right) jump event from box $i$, remove a particle from compartment $i$ and add one to compartment $i-1$ (resp. $i+1$). 2. 2.

If the event corresponds to a growth event: firstly, define the pregrowth state to be $\boldsymbol{r}\leftarrow\boldsymbol{m}$, calculate the overlap proportions $\delta_{i}^{k}=(k+1-i)/(k+1)$, and use these in order to draw $k$ binomial random variables $b_{i}\sim\text{Bin}(m_{i},\delta_{i}^{k})$. Create an extra compartment at the right end of the postgrowth domain (increasing $k$ by 1 by setting $k\leftarrow k+1$). Then set, for $j\in\{1,...,k\}$:

[TABLE]

• (0e)

If $t<T$ then return to step (\theAlgorithmb), else end.

We assess the stretching method by initialising a uniform profile and test to see whether uniformity is maintained under the stretching domain growth method. In Figure 6, we can see that the stretching method performs very well in comparison to the original method.

3.1 Case study: Morphogen gradient

For our case study, we apply the original and stretching methods to the formation of a morphogen gradient on an exponentially growing domain [20, 21]. We begin with an initial domain of length $L_{0}$, which grows with rate $\rho$, and on which particles with density $u(x,t)$ move and interact. Particles diffuse with diffusion coefficient $D$, and they decay uniformly at a rate $\mu$. There is also an influx of particles at the left-hand boundary, at rate $D\lambda$. We will only allow particles to jump between adjacent compartments and to interact within their own compartments. In order to compare the results of the stochastic simulation, we compute the solution to the associated mean-field PDE, which is given by:

[TABLE]

The PDE comprises three terms on the right-hand side — the first is the diffusive term, the second is dilution due to domain growth and the third is degradation of particles over the spatial domain [5]. The first boundary condition is the in-flux of particles at the left of the domain, while the second boundary condition is a reflective boundary. For our on-lattice simulations, diffusion is implemented by allowing particles to jump between neighbouring compartments with rate $d=D/h^{2}$. Particle degradation is achieved by allowing each particle in a compartment to be removed with rate $\mu$, while production is included through a production reaction in the compartment adjacent to the left-hand boundary, as is described by Taylor et al. [7].

We simulate this system with exponential growth rate $\rho=0.01$ on a domain of initial length $L_{0}=4$. We set the diffusion coefficient to be $D=0.0025$, the influx rate is specified by setting $\lambda=200$ and the degradation rate $\mu=0.005$. We simulate until a final time, $t=100$ and average over 50,000 repeats. The results are displayed in Figure 7.

As with the case of maintaining a uniform gradient, particle densities for the stretching method agree well with the associated mean-field PDE [15] (see Figures 7-7), however we still observe the same collection of particles at the boundaries for the original method (see Figures 7-7). This is particularly evident at the right-hand side of the domain, indicating that we are able to correctly simulate a reaction-diffusion system which incorporates first-order reactions using the stretching method. We anticipate that the extension to second- and higher-order reactions will yield similar results since the domain growth mechanisms is decoupled from the reaction mechanism.

3.2 Comparison of methods

In this section, we will investigate the two methods, and the parameter regimes in which the errors from the original method are acceptable due to the interplay between diffusion, the growth rate and initial domain length. We explore this in two ways. The first is through a heuristic argument.

The mean squared displacement describes how the variance in position of a Brownian particle changes in time. If multiple particles are initialised at the origin and diffuse for a time, $t$, then $\langle x^{2}\rangle=2Dt$, where the angled brackets denote an ensemble average of the squared distances from the origin. If a particle is to diffuse over the entire domain before the domain grows, then the squared distance from the origin would be $x^{2}\sim L_{0}^{2}$. Likewise, the typical time frame for growth is given by setting $t\sim 1/\rho$. Substituting these into the expression for the mean squared displacement yields $D\sim L_{0}^{2}\rho$. Therefore, we say we are in a ‘high diffusion’ parameter regime when $D>L_{0}^{2}\rho$. When in the low diffusion regime ($D<L_{0}^{2}\rho$), particles are unable to spread and equilibrate before the domain grows, which is when we see the build up of particles in the original method.

We verify this heuristic result with a second argument. In order to do this, we have simulated a non-dimensionalised stochastic system and compared it to the solution of the mean-field continuum diffusion equation in order to determine a threshold value for diffusion. The results can be seen in Figure 8, where we plot the histogram distance error (HDE), which is given by the expression

[TABLE]

Here, $p_{n}^{s}(D^{*})$ is the value of the normalised solution of the stochastic simulation (original or stretching method) at the final time with non-dimensional diffusion $D^{*}$, $p_{n}^{m}(D^{*})$ is the solution for the associated non-dimensional mean-field PDE, and $n$ indexes a common mesh on which we compare the two solutions. The non-dimensional diffusion parameter is equal to $D^{*}=D/(L_{0}^{2}\rho)$ where, as before, $L_{0}$ is the initial domain length and $\rho$ is the exponential growth rate. From Figure 8, a value $D^{*}$ greater than 1 yields similar HDE values for both the original and stretching methods. This indicates that, if $D>L_{0}^{2}\rho$ then the original method should have a similar performance to the stretching method. However, if the inequality is reversed, so that $D<L_{0}^{2}\rho$, then the stretching method should be used. We also note that the same pattern is apparent when comparing the solutions of the mean equations (2) and (24), which can be seen in Figure 8. The HDE for the stretching method is exactly zero because it maintains uniformity precisely in the mean field.

4 Discussion

Domain growth mechanisms for on-lattice models are of importance for the accurate representation of many biological processes. We have demonstrated beyond doubt that the original method, suggested by Baker et al. [15], causes a build up of particles at the boundaries of the spatial domain (see Figure 2). Consequently, we have developed a method for implementing domain growth when modelling reaction-diffusion systems at the mesoscale in order to correct this build-up. This technique involves stretching all compartments by a small amount (leading to the creation of a new compartment) and the appropriate re-distribution of the particles. We have demonstrated that this method agrees with the corresponding mean-field equations derived in the continuum limit, while maintaining a uniform profile, and have shown that it correctly models morphogen gradient formation on a growing domain.

The stretching method will be particularly useful when developing spatially extended hybrid methods on growing domains. These methods split the spatial domain into subdomains, in which different modelling paradigms are used, separated by an interface or overlap region [22]. We envisage that the stretching method can be used in the compartment-based subdomain of a hybrid model for reaction-diffusion on a growing domain, without causing a build-up of particles at the interface.

There are still several open questions regarding on-lattice domain growth whose answers go beyond the scope of this paper. The first of these concerns modelling domain growth in higher dimensions. To induce on-lattice domain growth in higher dimensions we can employ the following method by Binder et al. [23]. Consider, for the purposes of this example, a two-dimensional domain (although higher dimensional growth is straightforward to implemented by analogy). In order to maintain a rectangular domain, a “growth event” must increase the total number of rows or columns by one. For example, when a growth event is chosen to occur in the vertical direction, we must increase in the number of rows. To do this, we temporarily treat each column as its own one-dimensional domain, and implement a single vertical growth event in each column, independent of the others. Doing this for every column results in the whole domain increasing in height by a single row.

We have simulated such a domain growth process using the original method of Baker et al. [15] to implement the independent row or column elongations when carrying out a horizontal or vertical (respectively) growth event. Specifically, for clarity, we carry out horizontal and vertical growth events simultaneously, which maintains the aspect ratio of the initially square domain we begin with. Diffusion of particles is turned off in order to clearly demonstrate the bias induced by this two-dimensional version of domain growth as illustrated in Figure 9. The same effect that we have observed in one dimension (namely a preponderance of particles towards the boundaries of the domain) is also present in higher dimensions. Extending the stretching method will provide a straightforward fix to this problem in higher dimensions. Overlap fractions will be calculated as ratios of (hyper-)volumes as opposed to ratios of lengths, and particles in a pregrowth compartment will be distributed between multiple overlapping postgrowth compartments using multinomial distributions (the natural generalisation of the binomial distributions used in the one-dimensional case).

Whilst domain growth is important, equally, domain shrinkage is also of significant biological interest. Domain shrinkage through directed apoptosis (programmed cell death) is an important component of many wound healing processes [24, 25] for example. Without further investigation it is not immediately clear whether domain shrinkage implemented by the removal of a randomly chosen compartment (in analogy to the domain growth method of Baker et al. [15] and as implemented by Yates [3]) would induce bias in cell densities. In contrast we are confident that implementing domain shrinkage into the stretching method by considering pre- and post-shrink overlap regions compartments will not induce bias. Nevertheless these hypotheses remain to be tested.

Finally, it is of interest to adapt the stretching method to account for non-uniform growth, which has been shown to be important in biological scenarios [26]. All of the examples we have presented have implemented uniform growth — in which all regions of space grow at the same rate. The stretching method could be adapted to incorporate non-uniform growth by splitting the domain into groups of compartments each of which have a different growth rate. The stretching method can then be used on each of the groups individually.

Many authors have used the method presented by Baker et al. [15] to incorporate growth into on-lattice simulations of reaction-diffusion processes. For example, Woolley et al. [17] investigate the role that domain growth plays on modelling stochastic reaction-diffusion systems. Thompson et al. [18] explore cell migration and adhesion during biological development, while tissue growth and shrinkage are studied by Yates [3]. However, as demonstrated in this paper, particularly in the case of low diffusivity, the inherent bias in the domain growth method suggests that the conclusions drawn from these studies may require re-evaluation. We suggest that the ‘stretching’ domain growth method we propose in this paper is an appropriate alternative with which to re-evaluate these results and which should be employed in future studies of reaction-diffusion processes on growing domains.

Appendix A Deriving equation (5)

Within this appendix, we derive the expression in equation (5):

[TABLE]

which holds for $a,b\in\mathbb{N}$ such that $a<b$ and $b$ is large, and for $c\in\mathbb{R}$. We start by representing the product on the left-hand side of equation (5) using gamma functions:

[TABLE]

We apply Stirling’s approximation, which says that for large z, $\Gamma(z)\approx\sqrt{2\pi(z-1)}(z-1)^{z-1}e^{-(z-1)}$. Applying this to equation (A1), we obtain:

[TABLE]

We assume that $b$ is large, and thus, $b+c\approx b$. Applying this:

[TABLE]

Appendix B Deriving equation (9)

Within this appendix, we derive the expression in equation (9):

[TABLE]

where the approximation is for large values of $k$.

[TABLE]

where, in the final approximation, we have used equation (5) in order to simplify the two products. Simplifying yields:

[TABLE]

We now utilise equation (6) (approximation is better for large $k$), in order to obtain:

[TABLE]

Appendix C Discussion on local methods

Throughout Section 2, we have demonstrated that at low diffusion levels, the method presented by Baker et al. [15] fails to correctly model the growth of a uniform profile. This is a local method, which means that on every occasion that the domain is due to grow, we choose a compartment at random at which the growth event occurs. Over many repeats of the same process, different compartments will be chosen and so, when averaging over these repeats, each compartment is chosen an equal number of times (when considering only a single growth event).

There are three elements that define a local method:

The probability of choosing each compartment, 2. 2.

The direction of growth, 3. 3.

The redistribution of particles.

In the case of Baker et al. [15], we (1) choose each compartment uniformly at random, (2) always grow to the right and (3) redistribute the particles using a symmetric binomial distribution. We now discuss which of these we can change in order to create a local method that correctly maintains a uniform profile on a growing domain.

The simplest way of creating a different local method is to change one of the three elements. However, it can be shown that changing only a single element does not yield the expected uniform growth. As a result, the next simplest is to change two of the elements, whilst fixing one. One such way would be to (1) set the probability of choosing each compartment to be general, (2) set the direction of growth to be left or right with a probability of a half each, and (3) redistribute the particles using a binomial distribution with a generic probability of success. These properties yield an algorithm which has $k$ degrees of freedom (where there are $k-1$ pregrowth compartments), which can be used in order to solve a series of equations to ensure that a uniform profile is maintained.

In order to calculate the values for the probability distribution, and the probability for success in the binomial distribution when redistributing particles post-split, we can write a series of equations that relate the pre- and postgrowth states, on average. Using these, we are able to use numerical optimisation techniques in order to find the optimum values for the unknown parameters. However, while this method fixes the main issue with the original method, it introduces some new ones. Firstly, the centre of growth is no longer fixed at one end of the domain. This may cause problems for certain methodological applications (e.g. employing the method for spatial hybrid methods [27, 28, 29]) or for some biological applications in which tissues genuinely grow from a fixed origin at one end of the domain grow from one end of the domain (e.g. hyphal tip growth [30]). Secondly, we have to solve an overdetermined system for every possible number of compartments that might occur. This can be computationally expensive, especially for large compartment numbers. In practice, the unknown parameters can be computed and stored a priori, although, if the timing of domain growth events is stochastic, it may not be clear in advance exactly how many compartments the domain will comprise.

We also note that we have changed the direction of growth from being always to the right in the original method, to being left or right with equal probability. This is because there is no probability distribution for the compartments together with redistribution probability that maintains a uniform profile when the direction of growth is always the same.

Appendix D Justifying the stretching method

In Section 3, we introduced the stretching method as an alternative to the original method presented by Baker et al. [15]. In order to demonstrate that the original method fails, we have simulated the corresponding mean equations, and also analytically calculated the density of particles in each compartment with a large number of compartments. In this section, we will conduct a similar analysis of the stretching method.

D.1 The mean equations

In this section, we will calculate the mean equations for the stretching method by firstly considering the master equation for the process. We will then simulate the solutions to the mean equations, demonstrating that the solution remains uniform.

We begin by defining the sets $\mathcal{M}_{k}^{N}=\left\{\boldsymbol{m}\in\mathbb{N}^{k}:\sum_{i=1}^{k}{m_{i}}=N\right\}$. This is the set of all state vectors when there are $k$ compartments and $N$ particles in the system in total. We would like an expression for the probability that there are $k$ compartments in total, and the state variable is $\boldsymbol{m}\in\mathcal{M}_{k}^{N}$ with $N$ a fixed integer. We call this probability $p(\boldsymbol{m},k,t)$ as we did in Section 2. Then, defining $\pi(\boldsymbol{m}|\boldsymbol{r})$ to be the transition probability from state $\boldsymbol{r}$ to state $\boldsymbol{m}$:

[TABLE]

Here, the summand in the first term of the right hand side represents the rate of moving from a state $\boldsymbol{r}\in\mathcal{M}_{k-1}^{N}$ to the state $\boldsymbol{m}\in\mathcal{M}_{k}^{N}$. The second term is the rate at which the process leaves the state $\boldsymbol{m}$.

We next derive the mean equations. Similarly to the original method, define:

[TABLE]

Multiplying the CME (16) by $m_{i}$, summing over all possible $\boldsymbol{m}\in\mathcal{M}_{k}^{N}$, and applying equation (17), we obtain:

[TABLE]

From now on, we will drop the range of the sums, and simply write $\boldsymbol{m}$ or $\boldsymbol{r}$ and implicitly assume we are summing over the correct sets in order to simplify the notation. Recall that for each $\boldsymbol{r}$, there is an associated vector of binomial random variables $\boldsymbol{b}$ that re-distributes the particles from the pregrowth state to the postgrowth state. Further, for any $j\in\{1,...,k-1\}$, the probability of drawing $b_{j}$ is given by the probability mass function:

[TABLE]

where $\delta_{j}^{k-1}$ is the ‘overlap’ proportion defined in equation (13). We can then use the relationships set out in Algorithm 3, namely:

[TABLE]

Assuming that $i\in\{2,...,k-1\}$ (a similar argument can be applied to the case $i\in\{1,k\}$) and substituting the relationships in (20) into the sum in equation (18), we obtain:

[TABLE]

where in the final step, we have split the sum and also switched the order of summation. We will now make use of the binomial random variables (the $b_{j}$’s). We note that in order to transition from the pregrowth state $\boldsymbol{r}\in\mathcal{M}_{k-1}^{N}$ to the postgrowth state $\boldsymbol{m}\in\mathcal{M}_{k}^{N}$, we need to find a vector $\boldsymbol{b}$. However, only a small subset of these $\boldsymbol{b}$ vectors have a non-zero probability of occurring. Therefore, the sum over $\boldsymbol{m}$ can be re-written as a sum over the possible $\boldsymbol{b}$ values that have a non-zero probability. Using this, and substituting for $\pi(\boldsymbol{m}|\boldsymbol{r})$ the specific probability $\mathbb{P}(\boldsymbol{b}|\boldsymbol{r})$ of the binomial redistribution which takes us from state $\boldsymbol{r}$ to state $\boldsymbol{m}$, equation (21) becomes:

[TABLE]

Here $\mathbb{P}(\boldsymbol{b}|\boldsymbol{r})=\prod_{j=1}^{k-1}{\mathbb{P}(b_{j}|\boldsymbol{r})}$ (where $\mathbb{P}(b_{j}|r_{j})$ are given in equation (19)), and the sum over $\boldsymbol{b}$ is a sum over the set $\left\{\boldsymbol{b}:b_{j}\in\{0,...,r_{j}\}\ \forall j\in\{1,...,k-1\}\right\}.$ We now re-arrange the two sums in expression (22):

[TABLE]

The sum in the square brackets in the first term is equal to 1 as it is the sum of a probability distribution and the sum in the square brackets in the second term is a difference of two expectations:

[TABLE]

Finally, using the standard expectation for the binomial distribution, and bringing terms together, we find that:

[TABLE]

where the final equality uses the definition (17). Substituting this expression into equation (18) yields the mean-field density evolution equations for the stretching method:

[TABLE]

D.2 Analytical density

Here we will use a similar approach to Section 3 in order to calculate the average density of particles when the domain grows according to the stretching method. We let $u_{i}^{k}$ be the random density of particles in compartment $i$ when there are $k$ compartments in total, $M_{i}^{k}=\mathbb{E}[u_{i}^{k}]$, and $q(u_{i}^{k},\delta_{i}^{k})$ a realisation from a general probability distribution with mean value $M_{i}^{k}\delta_{i}^{k}$. The $q(u_{i}^{k},\delta_{i}^{k})$ values denote the density that is redistributed into postgrowth compartment $i$ from pregrowth compartment $i$, analogous to the $\mathbb{P}(b_{i}|r_{i})$ in equation (19) for the discretised process. Then the $u_{i}^{k}$’s are related by the following recursive relation:

[TABLE]

with an initial condition $u_{1}^{1}=1$. Here, the first term on the right-hand side is the fraction of the density that is provided to compartment $i$ on the postgrowth domain by the pregrowth compartment to the right, while the second term is the density provided from the left pregrowth compartment. We now take expectations to yield the recursive equations:

[TABLE]

We show that, under this recursion relation, $M_{i}^{k}=1/k$ by induction on the number of compartments, $k$. Clearly, using the initial condition, we have the base case $M_{1}^{1}=1$. Now assume that $M_{i}^{k-1}=1/(k-1)$ for all $i\in\{1,...,k-1\}$. Then, using (26):

[TABLE]

Here, the second line uses the inductive hypothesis and the third line employs the definition of $\delta_{i}^{k-1}$. This completes the inductive step, and hence the proof that $M_{i}^{k}=1/k$.