###### Abstract

Reaction–diffusion mechanism are a robust paradigm that can be used to
represent many biological and physical phenomena over multiple spatial scales.
Applications include intracellular dynamics, the migration of cells and the
patterns formed by vegetation in semi-arid landscapes. Moreover, domain growth
is an important process for embryonic growth and wound healing. There are many
numerical modelling frameworks capable of simulating such systems on growing
domains, however each of these may be well suited to different spatial scales
and particle numbers. Recently, spatially extended hybrid methods on static
domains have been produced in order to bridge the gap between these different
modelling paradigms in order to represent multiscale phenomena. However, such
methods have not been developed with domain growth in mind. In this paper, we
develop three hybrid methods on growing domains, extending three of the
prominent static domain hybrid methods. We also provide detailed algorithms to
allow others to employ them. We demonstrate that the methods are able to
accurately model three representative reaction-diffusion systems accurately
and without bias.

* AlgorithmexpList of Algorithms

Incorporating domain growth into hybrid methods for reaction-diffusion systems

Cameron A. Smith1,∗, Christian A. Yates1

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

Key words: Reaction–diffusion, domain growth, hybrid methods.

##  1 Introduction

The reaction–diffusion paradigm can be employed to model a range of biological
and physical scenarios over multiple length scales, from representing
vegetation patterns in semi-arid landscapes (Sherratt, 2005) and the study of
epidemics (Volpert and Petrovskii, 2009), to intracellular dynamics (Khan et
al., 2011; Andasari et al., 2012; ZhuGe et al., 2000). These systems couple
the random movement of particles (which when considered at a continuum level
manifests as the movement of particle density down the concentration gradient)
and the interaction of particles with each other and potentially with the
domain boundaries.

Domain growth is a process which underpins many biological processes, and it
is therefore important that we have accurate and efficient modelling methods
to represent it. Examples span many biological applications, including the
growth and shrinkage of tissue (Yates, 2014; Wolpert et al., 2015) and neural
crest cell migration (Mort et al., 2016; Kulesa et al., 2010; McLennan et al.,
2012). Domain growth has also been shown to play an important role in
theoretical studies of pattern formation in reaction-diffusion systems
(Crampin et al., 1999; Taylor et al., 2015; Woolley et al., 2011).

Reaction–diffusion systems can be modelled in several different ways, each of
which has different suitabilities depending on the scale of the system being
modelled. The coarsest of the three methods that we focus on is the
macroscale, which uses partial differential equations (PDEs) in order to
represent the system. PDEs model how the continuum density of particles evolve
in time, and are only suitable if the number of particles is high enough to
consider a continuum limit. The reaction–diffusion PDE on a growing domain
consists of four components (see equation (1)): a first order differential of
density with respect to time which describes the change of concentration in
time; a second-order differential of density with respect to space which
represents diffusion; a first-order differential of density with respect to
space representing domain growth, and finally a term which represents
reactions, if any are present. The macroscale is generally quick to implement
using a number of established techniques (for example, see Brenner and
Carstensen (2004); Eymard et al. (2000); Morton and Mayers (2005); Smith
(1985)), and there are often analytical approaches that can be employed to
investigate such systems. However, if particle numbers are too low, the
assumption that the continuum limit holds may break down, and stochastic
fluctuations may be found to play a more pivotal role. Moreover, the
deterministic mean-field PDE may not fully agree with its stochastic
counterparts described below. Attempts to derive a deterministic equivalent to
a non-linear stochastic model may result in an infinite hierarchy of
interrelated equations. This results in the need for moment closure, which
necessarily leads to the loss of some of the information encapsulated in the
higher-order moments.

The second modelling paradigm that we consider is the mesoscale, where we
split the spatial domain into a series of compartments, within which particles
reside. These particles are able to jump between neighbouring compartments,
and are able to interact/react with others within their own compartment. These
events are given exponential waiting times whose rates dictate the evolution
of the system. There are a large number of algorithms to simulate such
systems, both exact (Gillespie, 1977; Gibson and Bruck, 2000; Anderson, 2007)
and approximate (Gillespie, 2001; Auger et al., 2006). This middle scale is
generally slower than the macroscale, but provides more fine-grained accuracy
when particle numbers are lower. In order to incorporate domain growth at
least two methods that have been proposed. The first is a local method
proposed by Baker et al. (2010), which chooses a compartment uniformly at
random to instantaneously double in size and then divide. The particles in the
parent compartment are distributed into the two daughter compartments
according to some symmetric probability distribution. This method causes a
build-up of particles at the boundaries when growth occurs on a faster time-
scale than diffusion (Smith et al., 2019). The second method (and the method
we will employ in this work), is a global method introduced by Smith et al.
(2019). Compartments grow uniformly until a growth event is due to occur. At
this point, the boundaries between compartments are redrawn to include one
extra compartment, and the particles are redistributed appropriately. This
method works well in both high and low diffusion regimes.

Finally, at the microscale, we investigate Brownian-based dynamics. At this,
the smallest scale that we employ, individual particles are tracked and
updated in continuous space. There are several techniques that can be employed
in order to simulate a system at this scale, including the time-driven
mechanisms of Brownian motion for purely diffusive processes and Smoluchowski
dynamics (Smoluchowski, 1917) when reactions are involved, or the event-driven
Green’s function reaction dynamics (GFRD) (van Zon and ten Wolde, 2005). Under
time-driven algorithms, particles diffuse and are diluted according to an
appropriate stochastic differential equation (SDE). As well as each particle’s
location, if the system requires second- or higher-order reactions, we also
need to calculate the pairwise distances between all particles at every time-
step, which means that while this modelling paradigm is the most accurate, it
is also the slowest. Added to this, if the system in question is diffusion
limited, a very small time-step is required to accurately resolve the
dynamics. More efficient time-stepping is employed by the event-driven GFRD.
This sets a maximum time-step and the solution to the Smoluchowski equation in
order to combine diffusion and interactions, whilst accounting for the
additional error which is introduced in doing so.

There is another, even finer, scale of spatially resolved model known as
molecular dynamics, which we do not consider in this paper. We direct the
interested reader to (Holley, 1971; Dürr et al., 1981) for more information.

Many biological problems of interest are genuinely multiscale (Markevich et
al., 2004; Black and McKane, 2012; Gillespie et al., 2013; Robinson et al.,
2014). Consequently we require methods that are able to resolve the dynamics
at the appropriate scale. Spatially extended hybrid methods (Smith and Yates,
2018b) are one such class of techniques that are able to do this. These
methods employ two or more reaction-diffusion modelling paradigms (described
above) to represent the dynamics in different areas of the domain in the most
appropriate way. In regions with low particle numbers, for example, one of the
finer-grained stochastic methods might be employed at the cost of reduced
simulation efficiency. However, in regions of high particle numbers — enough
to consider a continuum limit — the more computationally efficient PDE may be
used. There are many examples of spatially extended hybrid methods (Yates and
Flegg, 2015; Moro, 2004; Spill et al., 2015; Flegg et al., 2015, 2012;
Robinson et al., 2014; Smith and Yates, 2018a; Franz et al., 2013; Alexander
et al., 2002; Erban, 2014, 2016), all of which focus on a static (non-growing)
domain. This paper extends three of these methods onto uniformly growing
domains. The pseudo-compartment method (PCM) (Yates and Flegg, 2015) is a
macroscopic-to-mesoscopic method which uses an interface to divide the domain
into two subdomains. Particles are able to jump between the two subdomains via
a ‘pseudo-compartment’ adjacent to the interface within the PDE subdomain. The
ghost cell method (GCM) proposed by Flegg et al. (2015) is a mesoscopic-to-
microscopic method that makes similar use of an extra compartment, coined the
ghost cell, adjacent to the interface in the microscopic subdomain in order to
couple compartment-based and Brownian dynamics. Finally, the auxiliary region
method (ARM) (Smith and Yates, 2018a) is a macroscopic-to-microscopic method
which, similar to the PCM and GCM, which employs a mesoscale auxiliary region
at the interface in order to allow particles to jump between the two
subdomians.

The rest of this paper will be set out as follows. In Section 2 we briefly
explain how the macroscopic, mesoscopic and microscopic can be considered
equivalent to each other. In Sections 3–5 we introduce the PCM, GCM and ARM,
explain the key differences between the algorithms for a static and growing
domain, and present the algorithms in full. We present representative results
for multiple test problems in Section 6, and discuss our findings in Section
7.

##  2 Equivalence Framework

In this section we present each of the three modelling paradigms that are the
constituent parts of our three hybrid methods. In Section 2.1 we introduce the
PDE approach. We present the mesoscopic approach in Section 2.2 and briefly
demonstrate how we can consider it to be equivalent to the PDE in the
appropriate limit. Finally, in Section 2.3 we do the same for the individual,
particle-based dynamics. All numerical algorithms can be found in the
Supplementary Material.

###  2.1 Macroscopic modelling

Firstly we consider the macroscale PDE representation. Consider a population
with density $u(x,t)$ undergoing diffusion at a position $x$ on an
exponentially growing one-dimensional domain $(0,L(t))\subseteq\mathbb{R}$ at
time $t>0$. We consider just a single dimension and exponential growth (with
rate $\rho$) for ease of description here, but the methods outlined can be
extended to higher dimensions, and will work for any form of uniform growth.
Under these assumptions, the length of the domain is $L(t)=L_{0}\exp\\{\rho
t\\}$, where $L_{0}$ is the initial length of the domain, and the
concentration evolves according to the following PDE:

$\frac{\partial u}{\partial t}(x,t)=D\frac{\partial^{2}u}{\partial
x^{2}}(x,t)-\rho\frac{\partial(xu(x,t))}{\partial x},$ (1)

which holds for $x\in(0,L(t))$ and for $t>0$. This description of domain
growth is known as the Eulerian representation in Eulerian coordinates, $x$
and $t$. The first term on the right-hand side represents the spread of
particles due to diffusion (with Fickian diffusion coefficient $D$) and the
second term is the dilution and spread of concentration caused by the
stretching of the spatial domain. We will employ this PDE in order to
demonstrate the equivalence of the two finer-scale methods to this first one.

In order to simulate this PDE, we need to switch from the Eulerian coordinates
above, in which the domain grows in time, to Lagrangian coordinates, in which
the domain remains static in time. We do this through the following change of
coordinates:

$\displaystyle x$ $\displaystyle=X\exp\\{\rho t\\},$ (2) $\displaystyle t$
$\displaystyle=\tau.$ (3)

Note now that the length of the domain in the Lagrangian spatial variable $X$
ranges between $0$ and the fixed value $L_{0}$. Further, the Lagrangian and
Eulerian temporal variables, $\tau$ and $t$ respectively, coincide here. There
are times when it is useful to rescale time, such as when finding analytical
solutions to the diffusion equation on uniformly growing domains (Simpson,
2015). For more details on the Lagrangian PDE and how we solve it numerically,
please see Section LABEL:SUPP-sect:Macro_Num of the supplementary material
(SM).

###  2.2 Mesoscopic modelling

In this section, we describe the mesoscale representation, which is the
“middle” scale that we will consider. In order to model at this level, we will
divide the spatial domain $(0,L(t))$ into a number of compartments labelled
$C_{i}(t)$ for $i\in\\{1,...,K(t)\\}$, where $K(t)$ is the (time-dependent)
number of compartments at time $t$. Particles lie within these compartments
and are able to jump between neighbouring compartments mimicking diffusion,
and can interact with one another within the same compartment through $R$
reaction channels. We will define the state of the system at time $t$ to be
$\boldsymbol{N}(t)$, where $N_{i}(t)$ is the number of particles at time $t$
in compartment $i$. Throughout this paper, we will implement the modified next
reaction method (Anderson, 2007) in order to advance the system forwards in
time. This method is used for explicit time-varying propensity functions (the
propensity function is a proxy for the rate of that event occurring), which
are important due to the domain growth. For an explanation of the method and
the algorithm, see the Supplementary Material, Section LABEL:SUPP-
sect:Meso_MNRM.

To extend the domain, we utilise the stretching method of Smith et al. (2019).
As the domain grows through a series of discrete fixed size extensions, the
compartments grow in length uniformly. Once a stretch event has been
determined to occur, the compartment boundaries are redrawn, making each
compartment smaller and making way for a new compartment. The particles are
then appropriately redistributed assuming a uniform distribution of particles
across each compartment. This algorithm can be included in the main mesoscale
algorithm in two ways. Firstly, propensity functions for the splitting events
(events where a new compartment is added) can be added to the list of
propensity functions, making this stochastic in time and space. Alternatively
(and the method we utilise here), we calculate the deterministic time at which
new compartments should be added to effect a particular domain growth pattern,
and add them at this time. In this way, the growth events are still stochastic
in space, but are now time-deterministic.

It can be shown, in the limit of small compartment size and fast inter-
compartment jumping rate, that the average behaviour of this mesoscopic model
can be described by the Eulerian PDE (1). We briefly present the main steps of
the calculation in the Supplementary Material (Section LABEL:SUPP-
sect:Meso_Equiv). For a more complete calculation please see (Smith et al.,
2019) for the mean equations and Section LABEL:SUPP-sect:Meso_Equiv in the
Supplementary Material for the small compartment limit.

###  2.3 Microscopic modelling

Within this section, we introduce the finest scale modelling paradigm that we
will employ for the hybrid methods presented in this paper. The microscale
tracks the individual locations of particles which diffuse and are
repositioned (due to domain growth) according to a stochastic differential
equation (SDE) as well as interacting depending on proximity to one another.
Particles diffuse through a Brownian motion, and the growth is implemented
using a deterministic drift term. In practice, this means that each particles
is “pulled along” as the domain stretches.

Let $X_{t}\in(0,L(t))$ be the position of a particle at time $t$. Then this
evolves according to:

$\text{d}X_{t}=\rho X_{t}\ \text{d}t+\sqrt{2D}\ \text{d}W_{t}.$ (4)

Here, the term on the left hand side denotes the change in position, the first
term on the right hand side is the drift term representing the repositioning
of the particles caused by the stretching of the domain, and the second term
on the right hand side is the diffusion. Note also that $\text{d}W_{t}$ is a
standard Weiner process.

In order to show that the density of particles evolving according to this SDE
is described by the Eulerian PDE (1), we look to the Fokker-Planck equation
(FPE) (otherwise known as the Kolmogorov forward equation (KFE)) (Risken,
1996). As with the mesoscopic case, we briefly outline the derivation in the
SM (Section LABEL:SUPP-sect:Micro_Num), and refer the interested reader to the
aforementioned reference.

##  3 The pseudo-compartment method

The first hybrid method we will adapt to incorporate uniform domain growth is
the macroscopic-to-mesoscopic pseudo-compartment method (PCM) (Yates and
Flegg, 2015). At time $t$ we decompose the domain as follows. The PDE
subdomain occupies the region $\Omega_{\text{{\tiny P}}}(t)=(0,I(t))$, where
$I(t)$ is the location of the interface at time $t$, and the mesoscopic
subdomain is $\Omega_{\text{{\tiny C}}}(t)=(I(t),L(t))$, where $L(t)$ is the
total length of the domain. The values of $I(t)$ and $L(t)$ will be calculated
deterministically from the initial position of these boundaries, and the
growth process, via the following pair of equations

$\displaystyle I(t)$ $\displaystyle=I(0)\exp\\{\rho t\\},$ $\displaystyle
L(t)$ $\displaystyle=L(0)\exp\\{\rho t\\}.$

Figure 1: A schematic of the static PCM (Yates and Flegg, 2015). The green
line denotes the density of particles in the macroscopic domain, while the
blue rectangles represent particles within each compartment. The red line is
the interface, and the green rectangles represent the number of “pseudo-
particles”, obtained by direct integration of the PDE solution over the
pseudo-compartment. The arrows crossing the interface denote the movement of
the pseudo-particles between the two subdomains.

In order to couple the macroscale and mesoscale, we allow particles to jump
between the first mesoscopic compartment (next to the interface (vertical red
line) in Figure 1) and the so-called “pseudo-compartment”, which is a region
of space the width of one compartment next to the interface on the PDE side of
the domain. Mass may only cross the interface through this jumping mechanism
and not though continuum PDE mass flowing over the boundary. As such, a zero-
flux boundary condition is imposed on the PDE at the interface.

Suppose that the number of particles in compartment $C_{1}(t)$ at time $t$ is
$n_{1}(t)$, and that it is of width $h_{c}(t)$. Then the region of space
occupied by the pseudo-compartment of length $h_{PC}(t)$ is denoted
$C_{-1}(t)$ and contains $n_{-1}(t)$ particles, where

$n_{-1}(t)=\int_{C_{-1}(t)}{u(x,t)\ dx}.$ (5)

In order to incorporate jumps into and out of the pseudo-compartment, we
introduce two events into the list of mesoscopic events, one for a jump into
the pseudo-compartment from compartment 1, and one for a jump in the other
direction. Since, on the growing domain these compartments are potentially of
different widths, the appropriate jumping rates would be $D/(h_{c}(t)^{2})$
and $D/(h_{PC}(t)h_{c}(t))$ respectively. The second of these takes into
account the differing compartment sizes. In the case that $h_{c}(t)$ and
$h_{PC}(t)$ are the same, this collapses to the usual diffusive jump rate.

The reasons that the length of a “regular” compartment $h_{c}(t)$ is typically
different from the length of the pseudo-compartment $h_{PC}(t)$ are subtle and
due to the differences between the Lagrangian (static) and Eulerian (growing)
coordinates. The PDE is being solved using a fixed mesh-width in Lagrangian
coordinates, and as such, the mesh grows when considered in Eulerian
coordinates. On the other hand, the compartment size is set in Eulerian
coordinates, and as a result, the number of Lagrangian PDE points per Eulerian
compartment decreases in time.

On the growing domain, the lengths of these compartments are calculated using
other simulation parameters. The length of the compartments $C_{i}(t)$ for
$i\in\\{1,...,K(t)\\}$ is given by the length of the mesoscopic part of the
domain, divided by the number of compartments $K(t)$. Therefore

$h_{c}(t)=\frac{L(t)-I(t)}{K(t)}.$ (6)

The pseudo-compartment size is calculated from the current compartment size as
follows. As well as having compartment properties we also need to solve the
PDE in the region occupied by the pseudo-compartment. As such, we will set its
length to be an integer number of PDE mesh points. The PDE mesh, when
considered from the Eulerian perspective, has width
$h_{p}(t)=h_{p}(0)\exp\\{\rho t\\}$, and we calculate the number of PDE mesh
points that would be in the pseudo-compartment, if it was of length
$h_{c}(t)$, to be

$p_{c}(t)=\frac{h_{c}(t)}{h_{p}(t)}.$

The value of $p_{c}(t)$ is generally not an integer, so in order to obtain an
integer we round (up or down or to the nearest integer) the value of
$p_{c}(t)$ and multiply this by $h_{p}(t)$ in order to find $h_{PC}(t)$.

If a particle is chosen to jump out of $C_{1}(t)$ and into $C_{-1}(t)$, we
firstly reduce $n_{1}(t)$ by one, and then add a particle’s worth of mass to
the pseudo-compartment according to:

$\boldsymbol{U}_{PC}=\boldsymbol{U}_{PC}+\frac{1}{h_{PC}(t)}\boldsymbol{1},$
(7)

where the vector $\boldsymbol{U}_{PC}$ is the numerical approximation to the
solution $u(x,t)$ at the PDE nodes contained within the pseudo-compartment,
and $\boldsymbol{1}$ is a vector of ones of the appropriate size. If a
particle jumps out of the pseudo-compartment and into $C_{1}(t)$, we add one
particle to $n_{1}(t)$ and remove a particle’s worth of mass uniformly across
the pseudo-compartment

$\boldsymbol{U}_{PC}=\boldsymbol{U}_{PC}-\frac{1}{h_{PC}(t)}\boldsymbol{1}.$
(8)

The algorithm for the implementation of the growing pseudo-compartment method
(gPCM) can be found in Algorithm 3, and a schematic for the static case is
given in in Figure 1. The algorithm is given for diffusion only, however
reactions may be incorporated through any appropriate method (see for example
(Erban and Chapman, 2009)).

Algorithm : The growing pseudo-compartment method (Diffusion only)

  * Initialise: Initial time — $t=0$; Final time — $t_{f}$; Initial compartment size — $h_{c}$; Initial pseudo-compartment size — $h_{pc}$; PDE solution — $\boldsymbol{U}$; Number of pseudo-particles — $n_{-1}$; Number of compartments — $K$; Compartment particle numbers — $\boldsymbol{n}$; Propensity functions — $a_{i}$ for $i\in\\{1,...,2K+1\\}$; Internal clock times — $T_{i}$ for $i\in\\{1,...,2K+1\\}$; Next firing times — $P_{i}$ for $i\in\\{1,...,2K+1\\}$; Times until next event — $\Delta t_{i}$ for $i\in\\{1,...,2K+1\\}$; PDE update step — $\Delta t$; Time until next PDE update — $t_{p}$; Time until next split event — $t_{s}$; Time until next PDE re-mesh event — $t_{r}$.

  * (0a)

At time $t>0$:

    1. 1.

Calculate $\Delta=\min\left\\{i\in\\{1,...,2K+1\\}:\Delta t_{i}\right\\}$ and
$\beta=\operatornamewithlimits{argmin}\left\\{i\in\\{1,...,2K+1\\}:\Delta
t_{i}\right\\}$. Set $t_{\Delta}=t+\Delta$.

    2. 2.

If $\min\\{t_{\Delta},t_{p},t_{s},t_{r}\\}=t_{\Delta}$:

      1. (a)

For every $i\in\\{1,...,2K+1\\}$, update $T_{i}$ according to

$T_{i}\leftarrow
T_{i}+\frac{1}{2\rho}a_{i}\left(1-\exp\\{-2\rho\Delta\\}\right).$

      2. (b)

For event $\beta$, set

$P_{\beta}\leftarrow P_{\beta}+\ln\left(\frac{1}{u_{1}}\right)\text{, where
}u_{1}\sim\text{Unif}(0,1).$

      3. (c)

Enact the event $\beta$:

         * •

If the event $\beta$ corresponds to a jump from the pseudo-compartment to the
first compartment, set $n_{1}\leftarrow n_{1}+1$ and set
$\boldsymbol{U}_{pc}\leftarrow\boldsymbol{U}_{pc}-1/h_{pc}\boldsymbol{1}$,
where $\boldsymbol{U}_{pc}$ are the pseudo-compartment nodes of the PDE
solution, and $\boldsymbol{1}$ is a vector of ones of the appropriate size.

         * •

If the event $\beta$ corresponds to a jump from the first compartment to the
pseudo-compartment, set $n_{1}\leftarrow n_{1}-1$ and set
$\boldsymbol{U}_{pc}\leftarrow\boldsymbol{U}_{pc}+1/h_{pc}\boldsymbol{1}$.

         * •

Otherwise, set
$\boldsymbol{n}\leftarrow\boldsymbol{n}+\boldsymbol{\nu}_{\beta}$, where
$\boldsymbol{\nu}_{\beta}$ is the stoichiometric vector for the event $\beta$.

      4. (d)

Set $t=t_{\Delta}$.

    3. 3.

Else if $\min\\{t_{\Delta},t_{p},t_{s},t_{r}\\}=t_{s}$:

      1. (a)

Enact a growth event according to Smith et al. (2019) (see Algorithm
LABEL:SUPP-alg:Stretch). Set $K^{\prime}=K+1$.

      2. (b)

For every $i\in\\{1,...,2K+1\\}$, update $T_{i}$ according to

$T_{i}\leftarrow
T_{i}+\frac{1}{2\rho}a_{i}\left(1-\exp\\{-2\rho(t_{s}-t)\\}\right).$

      3. (c)

For $i\in\\{2K+2,2K+3\\}$, set $T_{i}=0$ and
$P_{i}=\ln\left(1/u_{2,i}\right)$, where $u_{2,i}\sim\text{Unif}(0,1)$.

      4. (d)

Set $t=t_{s}$. Update $t_{s}$.

      5. (e)

Set $K\leftarrow K^{\prime}$.

    4. 4.

Else if $\min\\{t_{\Delta},t_{p},t_{s},t_{r}\\}=t_{p}$:

      1. (a)

For every $i\in\\{1,...,2K+1\\}$, update $T_{i}$ according to

$T_{i}\leftarrow
T_{i}+\frac{1}{2\rho}a_{i}\left(1-\exp\\{-2\rho(t_{p}-t)\\}\right).$

      2. (b)

Enact a PDE update step using Algorithm LABEL:SUPP-alg:PDE.

      3. (c)

Set $t\leftarrow t_{p}$. Set $t_{p}\leftarrow t_{p}+\Delta t$.

    5. 5.

Else:

      1. (a)

For every $i\in\\{1,...,2K+1\\}$, update $T_{i}$ according to

$T_{i}\leftarrow
T_{i}+\frac{1}{2\rho}a_{i}\left(1-\exp\\{-2\rho(t_{r}-t)\\}\right).$

      2. (b)

Re-mesh the PDE solution according to Algorithm LABEL:SUPP-alg:Remesh

      3. (c)

Set $h_{p}\leftarrow h_{p}/2$.

      4. (d)

Set $t\leftarrow t_{r}$. Update $t_{r}$ according to Algorithm LABEL:SUPP-
alg:Remesh.

    6. 6.

Update $h_{c}$ and $h_{pc}$.

    7. 7.

Update all propensity functions $a_{i}$, for $i\in\\{1,...,2K+1\\}$.

    8. 8.

Update $\Delta t_{i}$ according to

$\Delta
t_{i}=-\frac{1}{2\rho}\ln\left(1-\frac{2\rho(P_{i}-T_{i})}{a_{i}}\right),$

for $i\in\\{1,...,2K+1\\}$.

  * (0b)

If $t<t_{f}$, return to (a). Otherwise, end.

##  4 The ghost-cell method

Within this section, we describe the growing ghost cell method (gGCM), the
static counterpart of which was proposed by Flegg et al. (2015). This hybrid
method couples the mesoscopic and microscopic descriptions of reaction-
diffusion systems. We define the two spatial domains to be
$\Omega_{\text{{\tiny B}}}(t)=(0,I(t))$ for the microscopic, Brownian-based
dynamics, and $\Omega_{\text{{\tiny C}}}(t)=(I(t),L(t))$ for the compartment-
based subdomain. The definitions of $I(t)$ and $L(t)$ are the same as for
gPCM.

The coupling is implemented in a similar way to the PCM. A ghost cell (which
is analogous to the pseudo-compartment for the PCM) is created within the
microscopic subdomain, adjacent to the interface. Transport of mass over the
interface is implemented using the mesoscopic approach. As such, the
microscopic subdomain has a reflective boundary at the interface to ensure
that no particles are able to move across through that medium.

Figure 2: A schematic for the static GCM (Flegg et al., 2015). The blue
rectangles and red line are the same as in Figure 1. The yellow dots, denote
the positions of the individual particles, with arrows denoting the next jump
size and direction. Note that we have given each particle a different height
to aid clarity, but all particles lie on the axis in reality. The yellow
rectangles are the number of ghost cell particles, while the arrows over the
interface denote the direction of travel for the ghost cell particles. We
further note that this diagram has the two subdomains in the opposite order
compared to the description in the text.

Unlike in the PCM, since we don’t have two different discretisation lengths
being used in the numerical realisation of the algorithm (PDE and compartment
meshes), we are free to choose the ghost-cell to be the same size as the other
compartments in the simulation (as they change in size), and we set it to
$h_{c}(t)$ as defined in equation (6). The ghost cell therefore occupies the
region $C_{\text{{\tiny GC}}}(t)=(I(t)-h_{c}(t),I(t))$. The propensity
function for particles to jump from the ghost cell into the first compartment
in the compartment-based regime is $D/(h_{PC}(t)^{2})$, multiplied by the
number of particles $n_{\text{{\tiny GC}}}(t)$ in the ghost cell.
$n_{\text{{\tiny GC}}}(t)$ is calculated simply by counting the number of
particles within the ghost cell, so

$n_{\text{{\tiny GC}}}(t)=\sum_{i=1}^{n_{\text{{\tiny
B}}}(t)}{\mathbbm{1}_{[y_{i}(t)\in C_{\text{{\tiny GC}}}(t)]}},$

where $\mathbbm{1}_{[A]}$ is the indicator function that is one if $A$ is true
and 0 otherwise, $y_{i}(t)$ is the location of Brownian particle $i$ at time
$t$, and $n_{\text{{\tiny B}}}(t)$ is the total number of particles in the
Brownian-based subdomain at time $t$.

When a particle jumps from the ghost cell to the first compartment, we
increase $n_{1}(t)$ by one and remove one of the ghost cell particles
uniformly at random whilst simultaneously reducing $n_{\text{{\tiny GC}}}(t)$
by 1. When a jump occurs from the first compartment into the ghost cell, we
reduce $n_{1}(t)$ by one and add a new particle to the ghost cell by sampling
a uniform position within the ghost cell (that is, add a new particle at a
position $y^{*}\sim\text{Unif}(I(t)-h_{c}(t),I(t))$) and subsequently increase
$n_{\text{{\tiny GC}}}(t)$ by 1. Domain growth is implemented using the
stretching method (Smith et al., 2019) for the mesoscale, and via the SDE for
the microscale. The algorithm for the growing ghost-cell method can be found
in Algorithm 4, and a schematic for the static case is given in Figure 2.

Algorithm : The growing ghost cell method (Diffusion only)

  * Initialise: Initial time — $t=0$; Final time — $t_{f}$; Compartment size — $h_{c}$; Positions of particles — $\boldsymbol{y}$; Number of compartments — $K$; Compartment particle numbers — $\boldsymbol{n}$; Propensity functions — $a_{i}$ for $i\in\\{1,...,2K+1\\}$; Internal clock times — $T_{i}$ for $i\in\\{1,...,2K+1\\}$; Next firing times — $P_{i}$ for $i\in\\{1,...,2K+1\\}$; Times until next event — $\Delta t_{i}$ for $i\in\\{1,...,2K+1\\}$; Brownian update step — $\Delta t$; Time until next Brownian update — $t_{b}$; Time until next splitting event — $t_{s}$.

  * (0a)

At time $t>0$:

    1. 1.

Calculate $\Delta=\min\left\\{i\in\\{1,...,2K+1\\}:\Delta t_{i}\right\\}$ and
$\beta=\operatornamewithlimits{argmin}\left\\{i\in\\{1,...,2K+1\\}:\Delta
t_{i}\right\\}$. Set $t_{\Delta}=t+\Delta$.

    2. 2.

If $\min\\{t_{\Delta},t_{b},t_{s}\\}=t_{\Delta}$:

      1. (a)

For every $i\in\\{1,...,2K+1\\}$, update $T_{i}$ according to

$T_{i}\leftarrow
T_{i}+\frac{1}{2\rho}a_{i}\left(1-\exp\\{-2\rho\Delta\\}\right).$

      2. (b)

For event $\beta$, set

$P_{\beta}\leftarrow P_{\beta}+\ln\left(\frac{1}{u_{1}}\right)\text{, where
}u_{1}\sim\text{Unif}(0,1).$

      3. (c)

Enact the event $\beta$:

         * •

If the event $\beta$ corresponds to a jump from the ghost cell to the first
compartment, set $n_{1}\leftarrow n_{1}+1$ and remove a particle from the
ghost cell uniformly at random.

         * •

If the event $\beta$ corresponds to a jump from the first compartment to the
ghost cell, set $n_{1}\leftarrow n_{1}-1$ and add a new particle to the ghost
cell by drawing $u_{2}\sim\text{Unif}(0,1)$ and set the new particle’s
position $y^{*}$ to be $y^{*}=I(t)-u_{2}h_{c}(t)$.

      4. (d)

Set $t=t_{\Delta}$.

    3. 3.

Else if $\min\\{t_{\Delta},t_{b},t_{s}\\}=t_{s}$:

      1. (a)

Enact a growth event according to Smith et al. (2019) (see Algorithm
LABEL:SUPP-alg:Stretch in the SM). Set $K^{\prime}=K+1$.

      2. (b)

For every $i\in\\{1,...,2K+1\\}$, update $T_{i}$ according to

$T_{i}\leftarrow
T_{i}+\frac{1}{2\rho}a_{i}\left(1-\exp\\{-2\rho(t_{s}-t)\\}\right).$

      3. (c)

For $i\in\\{2K+2,2K+3\\}$, set $T_{i}=0$ and
$P_{i}=\ln\left(1/u_{3,i}\right)$, where $u_{3,i}\sim\text{Unif}(0,1)$.

      4. (d)

Set $t=t_{s}$. Update $t_{s}$.

      5. (e)

Set $K\leftarrow K^{\prime}$.

    4. 4.

Else:

      1. (a)

For every $i\in\\{1,...,2K+1\\}$, update $T_{i}$ according to

$T_{i}\leftarrow
T_{i}+\frac{1}{2\rho}a_{i}\left(1-\exp\\{-2\rho(t_{b}-t)\\}\right).$

      2. (b)

Enact a Brownian update step using algorithm LABEL:SUPP-alg:Ind-Update.

      3. (c)

Set $t\leftarrow t_{b}$. Set $t_{b}\leftarrow t_{b}+\Delta t$.

    5. 5.

Update $h_{c}(t)$.

    6. 6.

Update all propensity functions $a_{i}$, for $i\in\\{1,...,2K+1\\}$.

    7. 7.

Update $\Delta t_{i}$ according to

$\Delta
t_{i}=-\frac{1}{2\rho}\ln\left(1-\frac{2\rho(P_{i}-T_{i})}{a_{i}}\right),$

for $i\in\\{1,...,2K+1\\}$.

  * (0b)

If $t<t_{f}$, return to (a). Otherwise, end.

##  5 The auxiliary region method

Within this section, we describe the third and final of our growing hybrid
methods, the growing auxiliary region method (gARM) (Smith and Yates, 2018a).
This method couples the macro and microscales, using a similar methodology to
both the PCM and the GCM. Auxiliary regions are set up on either side of the
interface, which act as compartments for the purpose of allowing particles to
move between the two subdomains. A schematic for the static version of the
method is in Figure 3.

Figure 3: A schematic for the static ARM (Smith and Yates, 2018a). The green
and red lines are the same as in Figure 1, while the yellow dots and arrows
are the same as in Figure 2. The green and yellow rectangles denote the number
of PDE and Brownian auxiliary particles respectively. The arrows over the
interface denote the movement of these auxiliary particles.

The two auxiliary regions are equally sized compartments (although it is
possible to have them at different sizes), where particle numbers in the PDE
auxiliary region are calculated as in the pseudo-compartment of the gPCM for
the macroscopic subdomain, and particles in the Brownian auxiliary region are
calculated as in the ghost cell of the gGCM for the microscopic domain. As is
the case in each of those other hybrid methods, compartment-based jumping via
the auxiliary regions is the only mechanism by which particles may pass
between the subdomains. As such, the PDE requires a zero-flux boundary
condition at the interface, and the Brownian-based dynamics need an equivalent
reflective boundary.

In order to calculate the size of the auxiliary regions, we use a similar idea
to that employed in the gPCM. We set the size of a compartment to be $h$
initially, and use this to find the actual auxiliary region size,
$h_{\text{{\tiny AR}}}(t)$ by firstly calculating how many PDE mesh points lie
within this initial size, via:

$p_{c}(t)=\frac{h}{h_{p}(t)},$

where $h_{p}(t)$ is as in Section 3, the mesh-width in Eulerian coordinates.
We convert this to an integer number of PDE mesh points by rounding it to the
nearest integer, and convert it back to a length by multiplying by $h_{p}(t)$.
Using this, we can write the regions occupied by the auxiliary regions to be
$\Omega_{\text{{\tiny PA}}}(t)=(I(t)-h_{\text{{\tiny AR}}}(t),I(t))$ for the
PDE auxiliary region and $\Omega_{\text{{\tiny
BA}}}(t)=(I(t),I(t)+h_{\text{{\tiny AR}}}(t))$ for the Brownian-based
auxiliary region.

Particle numbers for the PDE auxiliary region, $n_{\text{{\tiny PA}}}(t)$, and
the Brownian auxiliary region, $n_{\text{{\tiny BA}}}(t)$, are calculated as
in the PCM and GCM respectively, and are given by

$\displaystyle n_{\text{{\tiny PA}}}(t)=\int_{\Omega_{\text{{\tiny
PA}}}(t)}{u(x,t)\ dx},$ (9) $\displaystyle n_{\text{{\tiny
BA}}}(t)=\sum_{i=1}^{N_{\text{{\tiny
B}}}(t)}{\mathbbm{1}_{[y_{i}(t)\in\Omega_{\text{{\tiny BA}}}(t)]}}.$ (10)

The implementation of the gARM is described in Algorithm 5 below.

Algorithm : The growing auxiliary region method (Diffusion only)

  * Initialise: Initial time — $t=0$; Final time — $t_{f}$; PDE solution — $\boldsymbol{U}$; Positions of particles — $\boldsymbol{y}$; Propensity functions — $a_{i}$ for $i=1,2$; Internal clock times — $T_{i}$ for $i=1,2$; Next firing times — $P_{i}$ for $i=1,2$; Times until next event — $\Delta t_{i}$ for $i=1,2$; PDE/Brownian update step — $\Delta t$; Time until next PDE/Brownian update — $t_{b}$; Time until next re-mesh event — $t_{r}$; Auxiliary region size $h_{\text{{\tiny AR}}}$.

  * (0a)

At time $t>0$:

    1. 1.

Calculate $\Delta=\min\left\\{i\in\\{1,2\\}:\Delta t_{i}\right\\}$ and
$\beta=\operatornamewithlimits{argmin}\left\\{i\in\\{1,2\\}:\Delta
t_{i}\right\\}$. Set $t_{\Delta}=t+\Delta$.

    2. 2.

If $\min\\{t_{\Delta},t_{b},t_{r}\\}=t_{\Delta}$:

      1. (a)

For every $i\in\\{1,2\\}$, update $T_{i}$ according to

$T_{i}\leftarrow
T_{i}+\frac{1}{2\rho}a_{i}\left(1-\exp\\{-2\rho\Delta\\}\right).$

      2. (b)

For event $\beta$, set

$P_{\beta}\leftarrow P_{\beta}+\ln\left(\frac{1}{u_{1}}\right)\text{, where
}u_{1}\sim\text{Unif}(0,1).$

      3. (c)

Enact the event $\beta$:

         * •

If the event $\beta$ corresponds to a jump from the PDE auxiliary region to
the microscopic auxiliary region, draw
$u_{2}\sim\text{Unif}(I(t)-h_{\text{{\tiny AR}}},I(t))$ and place a new
particle at that position, and set
$\boldsymbol{U}_{pc}\leftarrow\boldsymbol{U}_{pc}-1/h_{pc}\boldsymbol{1}$,
where $\boldsymbol{U}_{pc}$ are the pseudo-compartment nodes of the PDE
solution, and $\boldsymbol{1}$ is a vector of ones of the appropriate size.

         * •

If the event $\beta$ corresponds to a jump from the microscopic auxiliary
region to the PDE auxiliary region, choose one of the particles contained in
the microscopic auxiliary region uniformly at random and remove it, and set
$\boldsymbol{U}_{pc}\leftarrow\boldsymbol{U}_{pc}+1/h_{pc}\boldsymbol{1}$.

      4. (d)

Set $t=t_{\Delta}$.

    3. 3.

Else if $\min\\{t_{\Delta},t_{b},t_{r}\\}=t_{b}$:

      1. (a)

For $i\in\\{1,2\\}$, update $T_{i}$ according to

$T_{i}\leftarrow
T_{i}+\frac{1}{2\rho}a_{i}\left(1-\exp\\{-2\rho(t_{b}-t)\\}\right).$

      2. (b)

Enact a Brownian update step using algorithm LABEL:SUPP-alg:Ind-Update.

      3. (c)

Enact a PDE update step using algorithm LABEL:SUPP-alg:PDE.

      4. (d)

Set $t\leftarrow t_{b}$. Set $t_{b}\leftarrow t_{b}+\Delta t$.

    4. 4.

Else:

      1. (a)

For $i\in\\{1,2\\}$, update $T_{i}$ according to

$T_{i}\leftarrow
T_{i}+\frac{1}{2\rho}a_{i}\left(1-\exp\\{-2\rho(t_{r}-t)\\}\right).$

      2. (b)

Re-mesh the PDE solution according to Algorithm LABEL:SUPP-alg:Remesh

      3. (c)

Set $h_{p}\leftarrow h_{p}/2$.

      4. (d)

Set $t\leftarrow t_{r}$. Update $t_{r}$ according to Algorithm LABEL:SUPP-
alg:Remesh.

    5. 5.

Update $h_{\text{{\tiny AR}}}$.

    6. 6.

Update all propensity functions $a_{i}$, for $i\in\\{1,2\\}$.

    7. 7.

Update $\Delta t_{i}$ according to

$\Delta
t_{i}=-\frac{1}{2\rho}\ln\left(1-\frac{2\rho(P_{i}-T_{i})}{a_{i}}\right),$

for $i\in\\{1,2\\}$.

  * (0b)

If $t<t_{f}$, return to (a). Otherwise, end.

##  6 Results

Within this section, we present results from three test problems, for all
three of the methods described in Sections 3-5. The three test problems are
designed to evaluate the performance of the algorithms in comparison to the
corresponding PDE solutions. As such, the examples are relatively simple so
that the PDE is in exact correspondence with the expected behaviour of the
individual-based methods (i.e. no reactions of order higher than one). These
choices mean that discrepancies between the mean behaviour of the hybrid
methods and the PDE solutions can be be attributed directly to the
hybridisation. All examples will be on a one-dimensional, exponentially
growing domain, but can be straightforwardly extended to higher dimensions on
Cartesian domains, and to other forms of uniform domain growth.

The next three subsections will be devoted to the three test problems which
will each assess a different aspect of the performance of the three hybrid
algorithms.

###  6.1 Test problem 1: Maintaining uniformity

The first test problem verifies that the algorithms are able to maintain a
uniform particle distribution under pure diffusion. For this, we will use the
growing domain diffusion equation:

$\displaystyle\frac{\partial u}{\partial t}(x,t)$
$\displaystyle=D\frac{\partial^{2}u}{\partial
x^{2}}(x,t)-\rho\frac{\partial(xu(x,t))}{\partial x}$ $\displaystyle
x\in(0,2\exp\\{\rho t\\}),\ t>0,$ (11) $\displaystyle\frac{\partial
u}{\partial x}(0,t)$ $\displaystyle=0$ $\displaystyle t>0,$ (12)
$\displaystyle\frac{\partial u}{\partial x}(2\exp\\{\rho t\\},t)$
$\displaystyle=0$ $\displaystyle t>0$ (13) $\displaystyle u(x,0)$
$\displaystyle=\frac{M}{2}$ $\displaystyle x\in[0,2].$ (14)

Here, $M$ is the number of particles in the system, and all other parameters
are as in Sections 3—5. This PDE system has an analytical solution of

$u(x,t)=\frac{M}{2}\exp\\{-\rho t\\}$ (15)

We run this example with a diffusion coefficient of $D=0.0025$, an exponential
growth rate of $\rho=0.001$ and $M=500$ particles. Each hybrid simulation is
averaged over 1000 independent repeats for comparison and error plotting
purposes.

Figure 4: Results for Test Problem 1. We present results for the gPCM (Section
3) in column 1, the gGCM (Section 4) in column 2 and the gARM (Section 5) in
column 3. Figures 4-4 are snapshots of each method at the initial (row 1) and
final ($t=500$, row 2) times of the simulation for the three methods. Green
lines denote the PDE solution of the hybrid methods, blue bars are the
particle densities for the mesoscale within the hybrid methods, and the
microscale densities are denoted by the yellow bars, where we have binned
particles onto the same mesh as the compartments. The red vertical lines on
each plot denote the position of the interface at that time, and the black
dashed line is the solution of the PDE across the whole domain, which we
consider as our ground truth. Figures 4-4 display the relative errors in the
left subdomain (row 3) and right subdomain (row 4) corresponding to the
different modelling paradigms in the hybrid methods. The red curves denoting
the relative error are given by the formula (16) for the left subdomain, with
an analogous formula for the right side. The black dashed line corresponds to
an error of 0. There is no bias in the error in the positive or negative
direction for any of the methods.

In Figure 4, we present the results for each of the hybrid methods (gPCM in
column 1, gGCM in column 2 and gARM in column 3), with snapshots of the
solution at the initial time (row 1) and final time $t_{f}=500$ (row 2), and
then the relative errors in the left (row 3) and right (row 4) halves of the
domain. The relative error is calculated as

$E_{\text{{\tiny L}}}(t)=\frac{n_{\text{{\tiny L}}}^{H}(t)-n_{\text{{\tiny
L}}}^{P}(t)}{n_{\text{{\tiny L}}}^{P}(t)},$ (16)

for the left side of the domain. Here, $E_{\text{{\tiny L}}}(t)$ is the
relative error in the left side of the domain at time $t$, $N_{\text{{\tiny
L}}}^{H}(t)$ is the number of particles in the left subdomain of the hybrid
method at time $t$ and $n_{\text{{\tiny L}}}^{P}(t)$ is the number of
particles calculated in the left-hand side of the PDE solution at time $t$. As
can be seen from the plots, each hybrid method is able to correctly maintain
uniformity with no bias in particle numbers to either subdomain.

###  6.2 Test problem 2: Testing flux

The second test problem is designed to assess the ability of each of the
hybrid methods to cope with a non-zero flux acrosss the interface. For this,
we use the same PDE from test problem 1, but change the boundary conditions to
ensure a flux over the entire domain. The PDE with its boundary and initial
conditions are below:

$\displaystyle\frac{\partial u}{\partial t}(x,t)$
$\displaystyle=D\frac{\partial^{2}u}{\partial
x^{2}}(x,t)-\rho\frac{\partial(xu(x,t))}{\partial x}$ $\displaystyle
x\in(0,2\exp\\{\rho t\\}),\ t>0,$ (17) $\displaystyle-D\frac{\partial
u}{\partial x}(0,t)$ $\displaystyle=Ru(2\exp\\{\rho t\\},t)$ $\displaystyle
t>0,$ (18) $\displaystyle-D\frac{\partial u}{\partial x}(2\exp\\{\rho t\\},t)$
$\displaystyle=Ru(2\exp\\{\rho t\\},t)$ $\displaystyle t>0,$ (19)
$\displaystyle u(x,0)$ $\displaystyle=\frac{M}{2}$ $\displaystyle x\in[0,2].$
(20)

Equation (19) specifies a Robin boundary condition and equation (18) is the
corresponding condition which ensures periodic boundaries: particles that exit
the domain at the right-hand end will re-enter the domain at the left-hand
end. The resultant steady state (when the domain is not growing) is a linear
gradient in density.

Figure 5: Results for Test Problem 2. Figure descriptions are as in Figure 4.

We present our results for this test problem in Figure 5. As in the case with
the first test problem, we can see that all three of the hybrid methods
perform accurately with no discernible bias in particle numbers on either side
of the domain. This demonstrates that the hybrid methods are able to resolve
gradients over the interface.

###  6.3 Test problem 3: Morphogen gradient formation

For the final test problem, we investigate the formation of a morphogen
gradient on a growing domain. This example is designed to determine whether
the hybrid methods are able to accurately perform when zeroth- and first-order
reactions are incorporated111In this proof of principle paper we do not
include examples of second- or higher-order interactions. We know that the
mean-field PDE model does not correspond exactly to the mean of the stochastic
models in these cases due to the necessity for moment closure when deriving
the continuum equations from the individual-based models. To accurately
determine whether our hybrd coupling introduces bias we consider only examples
in which the expected behaviour of the individual-based method matches the
behaviour of the equivalent continuum model.. The PDE, boundary and initial
conditions are below:

$\displaystyle\frac{\partial u}{\partial t}(x,t)$
$\displaystyle=D\frac{\partial^{2}u}{\partial
x^{2}}(x,t)-\rho\frac{\partial(xu(x,t))}{\partial x}-\mu u$ $\displaystyle
x\in(0,2\exp\\{\rho t\\}),\ t>0,$ (21) $\displaystyle\frac{\partial
u}{\partial x}(0,t)$ $\displaystyle=-\lambda$ $\displaystyle t>0,$ (22)
$\displaystyle\frac{\partial u}{\partial x}(2\exp\\{\rho t\\},t)$
$\displaystyle=0$ $\displaystyle t>0,$ (23) $\displaystyle u(x,0)$
$\displaystyle=\frac{M}{2}$ $\displaystyle x\in[0,2].$ (24)

The regular diffusion and dilution in the PDE is augmented with a sink term,
$-\mu u$. This degradation of mass manifests as a first-order reaction of the
form:

$A\xrightarrow{\mu}\emptyset$

in both of the mesoscopic and microscopic representations, where $A$
represents a particle whose density is given by $u$. The boundary condition at
the left-hand end of the domain represents an influx of particles with rate
$D\lambda$. This can be thought of as a zeroth-order reaction at the left-hand
boundary, of the form:

$\emptyset\xrightarrow{\kappa}A,$

where $\kappa$ is the rate of introduction of new particles which is related
to $\lambda$ via $\kappa=\lambda D$.

Figure 6: Results for Test Problem 3. Figure descriptions are as in Figure 4.

We present the results for test problem 3 in Figure 6. As with each of the
previous test problems, we see that there is a very good qualitative agreement
with all three methods and the ground truth numerical PDE solution (6-6) and
this is further confirmed by the relative error plots, which show no
discernible bias in either direction for any of the hybrid methods.

##  7 Discussion

In this paper, we have developed three spatially extended hybrid methods which
are capable of simulating multi-scale reaction diffusion processes on
uniformly growing domains. Each method has been extended from previous methods
originally specified on static domains (Yates and Flegg, 2015; Flegg et al.,
2015; Smith and Yates, 2018a). We have provided descriptions, schematics and
detailed algorithms for the implementation of these methods. Furthermore, we
have demonstrated that each of these methods are accurate and unbiased under
three test problems: pure diffusion with zero-flux boundaries, diffusion with
partially absorbing periodic boundaries, and the formation of a morphogen
gradient with particle degradation and an influx boundary.

We have focussed on exponential growth, however it should be noted that the
methods set out here will work for any type of uniform growth, provided that
there is an equivalence framework between the three individual modelling
paradigms. For example, if the domain were to grow linearly with rate $\rho$,
the domain length at time $t$ (given it is of length $L_{0}$ to begin with)
would be $L(t)=L_{0}(1+\rho t)$. The advection term in equation (1) would
become $\rho(ux)_{x}/(1+\rho t)$, where the subscript here denotes
differentiation with respect to the $x$ variable. In order for the mesoscale
to be equivalent, domain growth events must occur with rate $\rho/(h_{c}K)$
(Baker et al., 2010). We must also be careful to adapt the calculation of the
internal times of the mesoscopic events in each of the algorithms for this
altered domain growth scenario (Anderson, 2007).

We have also presented each of our examples here in one dimension for
simplicity and clarity of explanation, however, it would be straightforward to
extend the methods to higher dimensions which have (hyper-)planar interfaces.
Care must be taken if domain growth is to be implemented in higher dimensions,
that the interface remains coherent. For example, when implementing domain
growth in the compartment-based method, it makes sense to employ deterministic
growth as in (Smith et al., 2019) (rather than stochastic growth (Baker et
al., 2010)) so that the rate of domain growth along the interface matches in
each of the coupled methodologies.

There are several extensions (many adapted from the static hybrid literature)
that might be included in order to render these growing hybrid methods more
versatile. The inclusion of an adaptive interface, for example, would allow
the methods to more robustly and efficiently deal with examples in which the
concentration profile changes significantly. Static interfaces (interfaces
which remain in the same place when considered in Lagrangian coordinates) are
useful when we either know something about what the solution looks like a
priori, and hence we can place the interface in an appropriate position, or if
there is a region of the spatial domain that requires a more detailed
representation (such as around ion channels when considering transmembrane
transport - see, for example (Dobramysl et al., 2015)). However, if neither of
these are true, an adaptive interface may be required. Implementing an
adaptive interface requires the use of local densities around the interface in
order to determine its position; if the density in the finer scale subdomain
is too high, the interface can move into that subdomain - extending the
coarser subdomain, and vice-versa if the local density in the coarser
subdomain is too low (Robinson et al., 2014; Harrison and Yates, 2016; Smith
et al., 2019).

While domain growth is an important phenomena in many biological scenarios
(Baker et al., 2010), domain shrinkage also has equally important
applications. As an example, directed apoptosis is an important component of
wound healing, and requires models that incorporate domain shrinkage (Yates,
2014). The hybrid methods presented in this paper would extend equally well to
a domain that uniformly shrinks in size as they do to those that grow.

The hybrid methods presented here provide accurate ways of simulating
reaction-diffusion systems on uniformly growing domains, and will be of
interest to those who model such systems that either have a scale difference
in particle numbers across their domain, or have regions of space which need
to be modelled in more detail than others. The methods developed here will
allow members of the modelling community to probe the important effects of
stochasticity in their multiscale systems without the suffering potential
simulation penalties that modelling the entire system using a fine-scale
methodology might bring.

## References

  * Alexander et al. (2002) F.J. Alexander, A.L. Garcia, and D.M. Tartakovsky.  Algorithm refinement for stochastic partial differential equations: I. linear diffusion.  _J. Comput. Phys._ , 182(1):47–66, 2002.
  * Andasari et al. (2012) V. Andasari, R.T. Roper, M.H. Swat, and M.A.J. Chaplain.  Integrating intracellular dynamics using CompuCell3D and Bionetsolver: applications to multiscale modelling of cancer cell growth and invasion.  _PLoS One_ , 7(3):e33726, 2012.
  * Anderson (2007) D.F. Anderson.  A modified next reaction method for simulating chemical systems with time dependent propensities and delays.  _J. Chem. Phys._ , 127(21):214107, 2007.
  * Auger et al. (2006) A. Auger, P. Chatelain, and P. Koumoutsakos.  R-leaping: Accelerating the stochastic simulation algorithm by reaction leaps.  _J. Chem. Phys._ , 125:084103, 2006.
  * Baker et al. (2010) R.E. Baker, C.A. Yates, and R. Erban.  From microscopic to macroscopic descriptions of cell migration on growing domains.  _Bull. Math. Biol._ , 72(3):719–762, 2010.
  * Black and McKane (2012) A.J. Black and A.J. McKane.  Stochastic formulation of ecological models and their applications.  _Trends Ecol. Evol._ , 27(6):337–345, 2012.
  * Brenner and Carstensen (2004) S.C. Brenner and C. Carstensen.  _Finite element methods_ , 2004.
  * Crampin et al. (1999) E.J. Crampin, E.A. Gaffney, and P.K. Maini.  Reaction and diffusion on growing domains: scenarios for robust pattern formation.  _Bull. Math. Biol._ , 61(6):1093–1120, 1999.
  * Dobramysl et al. (2015) U. Dobramysl, S. Rüdiger, and R. Erban.  Particle-based multiscale modeling of intracellular calcium dynamics.  _Multiscale Modeling & Simulation_, 14(3):997–1016, 2015.
  * Dürr et al. (1981) D. Dürr, S. Goldstein, and J.L. Lebowitz.  A mechanical model of brownian motion.  _Commun. Math. Phys._ , 78(4):507–530, 1981.
  * Erban (2014) R. Erban.  From molecular dynamics to Brownian dynamics.  _Proc. R. Soc. A_ , 470(2167):20140036, 2014.
  * Erban (2016) R. Erban.  Coupling all-atom molecular dynamics simulations of ions in water with brownian dynamics.  In _Proc. R. Soc. A_ , volume 472, page 20150556. The Royal Society, 2016.
  * Erban and Chapman (2009) R. Erban and S.J. Chapman.  Stochastic modelling of reaction–diffusion processes: algorithms for bimolecular reactions.  _Phys. Biol._ , 6(4):1–18, 2009.
  * Eymard et al. (2000) R. Eymard, T. Gallouët, and R. Herbin.  Finite volume methods.  _Handbook of numerical analysis_ , 7:713–1018, 2000.
  * Flegg et al. (2012) M.B. Flegg, S.J. Chapman, and R. Erban.  The two-regime method for optimizing stochastic reaction–diffusion simulations.  _J. Roy. Soc. Interface._ , 9(70):859–868, 2012\.
  * Flegg et al. (2015) M.B. Flegg, S. Hellander, and R. Erban.  Convergence of methods for coupling of microscopic and mesoscopic reaction-diffusion simulations.  _J. Comput. Phys._ , 289(C):1–17, 2015.
  * Franz et al. (2013) B. Franz, M.B. Flegg, S.J. Chapman, and R. Erban.  Multiscale reaction-diffusion algorithms: PDE-assisted Brownian dynamics.  _SIAM J. Appl. Math._ , 73(3):1224–1247, 2013\.
  * Gibson and Bruck (2000) M.A. Gibson and J. Bruck.  Efficient exact stochastic simulation of chemical systems with many species and many channels.  _J. Phys. Chem. A._ , 104(9):1876–1889, 2000\.
  * Gillespie (1977) D.T. Gillespie.  Exact stochastic simulation of coupled chemical reactions.  _J. Phys. Chem._ , 81(25):2340–2361, 1977.
  * Gillespie (2001) D.T. Gillespie.  Approximate accelerated stochastic simulation of chemically reacting systems.  _J. Chem. Phys._ , 115(4):1716–1733, 2001.
  * Gillespie et al. (2013) D.T. Gillespie, A. Hellander, and L.R. Petzold.  Perspective: Stochastic algorithms for chemical kinetics.  _J. Chem. Phys._ , 138(17):170901, 2013.
  * Harrison and Yates (2016) J.U. Harrison and C.A. Yates.  A hybrid algorithm for coupling PDE and compartment-based dynamics.  _J. Roy. Soc. Interface_ , 13:20160335, 2016.
  * Holley (1971) R. Holley.  The motion of a heavy particle in an infinite one dimensional gas of hard spheres.  _Probab. Theory. Rel._ , 17(3):181–219, 1971\.
  * Khan et al. (2011) S. Khan, Y. Zou, A. Amjad, A. Gardezi, C.L. Smith, C. Winters, and T.S. Reese.  Sequestration of camkii in dendritic spines in silico.  _J. Comput. Neurosci._ , 31(3):581–594, 2011\.
  * Kulesa et al. (2010) P.M. Kulesa, C.M. Bailey, J.C. Kasemeier-Kulesa, and R. McLennan.  Cranial neural crest migration: new rules for an old road.  _Dev. Biol._ , 344(2):543–554, 2010.
  * Markevich et al. (2004) N.I. Markevich, J.B. Hoek, and B.N. Kholodenko.  Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascades.  _J. Cell Biol._ , 164(3):353–359, 2004.
  * McLennan et al. (2012) R. McLennan, L. Dyson, K.W. Prather, J.A. Morrison, R.E. Baker, P.K. Maini, and Kulesa P.M.  Multiscale mechanisms of cell migration during development: theory and experiment.  _Development_ , 139(16):2935–2944, 2012.
  * Moro (2004) E. Moro.  Hybrid method for simulating front propagation in reaction-diffusion systems.  _Phys. Rev. E_ , 69(6):060101, 2004.
  * Mort et al. (2016) R.L. Mort, R.J.H. Ross, K.J. Hainey, O.J. Harrison, M.A. Keighren, G. Landini, R.E. Baker, K.J. Painter, I.J. Jackson, and C.A. Yates.  Reconciling diverse mammalian pigmentation patterns with a fundamental mathematical model.  _Nat. Commun._ , 7:10288, 2016.
  * Morton and Mayers (2005) K.W. Morton and D.F. Mayers.  _Numerical Solution of Partial Differential Equations_.  Cambridge University Press, 2005.
  * Risken (1996) H. Risken.  Fokker-planck equation.  In _The Fokker-Planck Equation_ , pages 63–95. Springer, 1996.
  * Robinson et al. (2014) M. Robinson, M.B. Flegg, and R. Erban.  Adaptive two-regime method: application to front propagation.  _J. Chem. Phys._ , 140(12):124109, 2014.
  * Sherratt (2005) J.A. Sherratt.  An analysis of vegetation stripe formation in semi-arid landscapes.  _J. Math. Biol._ , 51(2):183–197, 2005.
  * Simpson (2015) M.J. Simpson.  Exact solutions of linear reaction-diffusion processes on a uniformly growing domain: Criteria for successful colonization.  _PLoS One_ , 10(2):e0117949, 2015.
  * Smith and Yates (2018a) C.A. Smith and C.A. Yates.  The auxiliary region method: A hybrid method for coupling a PDE to Brownian-based dynamics for reaction-diffusion systems.  _R. Soc. Open Sci._ , 5(8), 2018a.
  * Smith and Yates (2018b) C.A. Smith and C.A. Yates.  Spatially extended hybrid methods: a review.  _J. Roy. Soc. Interface_ , 15(139), 2018b.
  * Smith et al. (2019) C.A. Smith, C Mailler, and C.A Yates.  Unbiased on-lattice domain growth.  _Phys. Rev. E_ , 100(6), dec 2019.
  * Smith (1985) G.D. Smith.  _Numerical solution of partial differential equations: finite difference methods_.  Oxford University Press, 1985.
  * Smoluchowski (1917) M. Smoluchowski.  Versuch einer mathematischen theorie der koagulationskinetik kolloider lösungen.  _Z. Phys. Chem._ , 92(129-168):9, 1917.
  * Spill et al. (2015) F. Spill, P. Guerrero, T. Alarcon, P.K. Maini, and H. Byrne.  Hybrid approaches for multiple-species stochastic reaction–diffusion models.  _J. Comput. Phys._ , 299:429–445, 2015.
  * Taylor et al. (2015) P.R. Taylor, R.E. Baker, and C.A. Yates.  Deriving appropriate boundary conditions, and accelerating position-jump simulations, of diffusion using non-local jumping.  _Phys. Biol._ , 12(1):016006, 2015.
  * van Zon and ten Wolde (2005) J.S. van Zon and P.R. ten Wolde.  Green’s-function reaction dynamics: a particle-based approach for simulating biochemical networks in time and space.  _J. Chem. Phys._ , 123(23):234910, 2005.
  * Volpert and Petrovskii (2009) V. Volpert and S. Petrovskii.  Reaction–diffusion waves in biology.  _Phys. Life Rev._ , 6(4):267–310, 2009.
  * Wolpert et al. (2015) L. Wolpert, C. Tickle, and A.M. Arias.  _Principles of development_.  Oxford University Press, 5 edition, 2015.
  * Woolley et al. (2011) T.E. Woolley, R.E. Baker, E.A. Gaffney, and P.K. Maini.  Stochastic reaction and diffusion on growing domains: understanding the breakdown of robust pattern formation.  _Phys. Rev. E_ , 84(4):046216, 2011.
  * Yates (2014) C.A. Yates.  Discrete and continuous models for tissue growth and shrinkage.  _J. Theor. Biol._ , 350:37–48, 2014.
  * Yates and Flegg (2015) C.A. Yates and M.B. Flegg.  The pseudo-compartment method for coupling partial differential equation and compartment-based models of diffusion.  _J. Roy. Soc. Interface_ , 12(106):20150141, 2015\.
  * ZhuGe et al. (2000) R. ZhuGe, K.E. Fogarty, R.A. Tuft, L.M. Lifshitz, K. Sayar, and J.V. Walsh.  Dynamics of signaling between Ca2+ sparks and Ca2+-activated K+ channels studied with a novel image-based method for direct intracellular measurement of ryanodine receptor Ca2+ current.  _J. Gen. Physiol._ , 116(6):845–864, 2000.