# Exact solutions for diffusive transport on heterogeneous growing domains

Stuart T. Johnston  School of Mathematics and Statistics, The University of
Melbourne, Australia<EMAIL_ADDRESS>Matthew J. Simpson
School of Mathematical Sciences, Queensland University of Technology,
Brisbane, Australia.

###### Abstract

From the smallest biological systems to the largest cosmological structures,
spatial domains undergo expansion and contraction. Within these growing
domains, diffusive transport is a common phenomenon. Mathematical models have
been widely employed to investigate diffusive processes on growing domains.
However, a standard assumption is that the domain growth is spatially uniform.
There are many relevant examples where this is not the case, such as the
colonisation of growing gut tissue by neural crest cells. As such, it is not
straightforward to disentangle the individual roles of heterogeneous growth
and diffusive transport. Here we present exact solutions to models of
diffusive transport on domains undergoing spatially non-uniform growth. The
exact solutions are obtained via a combination of transformation, convolution
and superposition techniques. We verify the accuracy of these solutions via
comparison with simulations of a corresponding lattice-based random walk. We
explore various domain growth functions, including linear growth, exponential
growth and contraction, and oscillatory growth. Provided the domain size
remains positive, we find that the derived solutions are valid. The exact
solutions reveal the relationship between model parameters, such as the
diffusivity and the type and rate of domain growth, and key statistics, such
as the survival and splitting probabilities.

###### keywords:

Random walk , Diffusion , Survival probability , Splitting probability ,
Growing domain , Exact solutions.

††journal: arXiv

##  1 Introduction

The expansion and contraction of spatial domains occurs throughout nature [1,
2, 3, 4]. At the cosmological scale, the universe is undergoing an
accelerating expansion [3]; at the organism scale, regular development
requires tissue growth [4]; at the cellular scale, cardiomyocytes expand and
contract with every heartbeat [1]. Diffusive processes regularly occur within
or along such domains. Tissue growth involves the migration of cell
populations, which can be considered as a diffusive process [2]. Throughout
the intracellular cardiomyocyte environment, chemical species diffuse while
driving cellular function [1]. Regulation of diffusive processes within
expanding and contracting domains can be critical to healthy function and
form. The failure of neural crest cells to colonise growing gut tissue during
development leads to Hirschsprung’s disease, which is a life-threatening birth
defect characterised by a partial or full absence of the enteric nervous
system [2].

While a common assumption, domain expansion and contraction is not necessarily
a spatially uniform process [5, 6, 7]. Binder _et al._ [5] demonstrate that
the growth rate of quail gut tissue during development is different between
the two midgut regions and the hindgut region. Accordingly, the migration of
quail neural crest cells occurs on a domain that undergoes spatially non-
uniform growth. Multilayered drug-loaded nanocapsules exhibit temperature-
dependent swelling, which alters the domain through which drug molecules
diffuse [8].

Mathematical models are widely used to investigate diffusive processes [9].
However, in the context of diffusive processes on multiple growing domains,
previous investigations have typically focused on either a single uniformly
growing domain [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20] or multiple non-
growing domains [21, 22, 23]. For example, Crampin _et al._ [10, 11]
investigate the conditions sufficient for pattern formation to occur on a
uniformly growing domain, which is relevant for biological morphogenesis.
Simpson _et al._ [24, 19], Ryabov _et al._ [17] and Yuste _et al._ [20] each
present exact results for specific models of diffusive processes on a
uniformly growing domain. The exact results allow key statistics such as the
survival probability (the probability that an individual has yet to cross the
boundary of an expanding domain) to be readily calculated [24, 19, 20].
Approaches have been presented to investigate diffusive processes in spatially
heterogeneous non-growing domains [23], including homogenisation techniques
and exact solutions for multilayered piecewise homogeneous domains [22]. One
approach for representing spatially non-uniform domain growth is via multiple
domains that exhibit piecewise uniform growth [5, 6]. Numerical results have
been presented in the context of pattern formation during non-uniform domain
growth, where the non-uniform growth is a piecewise uniform process [6]. While
both single growing domain problems and multiple non-growing domain problems
are both well-studied, there is a lack of exact results detailing the dynamics
of diffusive processes on multiple growing domains. Accordingly, it is unclear
how spatially non-uniform domain growth may inhibit or enhance diffusive
transport. Understanding this interplay between diffusive transport and non-
uniform domain growth may provide insight into why normal development
processes succeed or fail, and allow the calculation of the effective release
rate of therapeutics from nanocapsules.

Here we introduce new exact solutions for density profiles, survival
probabilities and splitting probabilities for diffusive processes on multiple
growing domains. The solutions are obtained via a combination of
transformation, convolution and superposition techniques. We demonstrate the
solutions are valid through comparison with density profiles obtained via
repeated realisations of a corresponding lattice-based random walk. We show
that the solution profiles can exhibit jump discontinuities at inter-domain
boundaries and investigate how model parameters, such as the domain growth
rate, diffusivity and initial location, influence the survival and splitting
probabilities. Critically, each of these parameters can vary between domains,
which gives rise to rich behaviour that is not possible on a single uniformly
growing domain. We investigate a suite of domain growth functions, including
linear growth [19], exponential growth and decay [25, 26, 19, 27], and
oscillatory growth [28] and show that the exact solutions are valid, provided
that the domain size does not reduce to zero.

##  2 Model

Figure 1: (a) Configuration of the domain in the discrete model. Agents exist
on a lattice of width $\Delta$. The pink line denotes the boundary between two
domains. (b) Lattice configuration before and after possible movement events
for agents in either domain. The agent that moves is highlighted in orange.
(c) Lattice configuration before and after possible positive domain growth
events. The location where the new lattice site is added is highlighted in
orange. (d) Lattice configuration before and after possible negative domain
growth events. The location where the lattice site is removed is highlighted
in orange. (e) Space-time diagram of a growing domain, highlighting the
expansion of space. The location of the boundaries are shown in grey.
Characteristics for domain one, two and three are shown in cyan, orange and
pink, respectively. For visual clarity we present characteristics for a
linearly growing domain; in practice, other types of domain growth are
possible.

###  2.1 Agent-based discrete model

We implement a one-dimensional position-jump random walk model on a growing
domain [19]. Agents in the random walk are located on discrete lattice sites
at $x=i\Delta,\ i\in\\{-\lfloor L(j\delta t)/\Delta\rfloor,\ \ldots,\ \lfloor
L(j\delta t)/\Delta\rfloor\\}$ where $L(j\delta t)$ is the half-length of the
growing domain after $j\in\mathcal{N}_{0}$ timesteps of duration $\delta t$
and $\Delta$ is the lattice size. There are no restrictions on the number of
agents that can occupy a single lattice site; we note other investigations
impose such a restriction [16, 29]. Agents located on the site at $x$ randomly
move to one of the two nearest-neighbour sites at $x\pm\Delta$ with
probability $P/2$ in a timestep. If selected, each agent will therefore move
with probability $P$ in a timestep. We select $N(j\delta t)$ agents at random,
with replacement, each timestep, where $N(j\delta t)$ is the number of agents
on the lattice. If an agent crosses the boundary at $x=\pm\Delta\lfloor
L(j\delta t)/\Delta\rfloor$ (i.e. the agent moves to a “ghost” lattice site at
$x=\pm\Delta\lfloor L(j\delta t)/\Delta+1\rfloor$) then the agent is
determined to have left the domain, and the agent is removed from the random
walk at the end of the timestep.

We first consider a single domain undergoing spatially-uniform growth. That
is, each location in space experiences uniform growth and there is a single
growth rate across the domain [19]. We implement positive domain growth in the
random walk as follows. When $\lfloor L(j\delta t)/\Delta\rfloor>\lfloor
L((j-1)\delta t)/\Delta\rfloor$, we randomly select a lattice site $k\in\\{1,\
\ldots,\ \lfloor L((j-1)\delta t)/\Delta\rfloor\\}$. We insert a new lattice
site at $x=k\Delta$, and hence all agents located at sites $i\geq k$ are
displaced a distance $\Delta$ in the positive direction [19, 27]. We repeat
this process for $k\in\\{-\lfloor L((j-1)\delta t)/\Delta\rfloor,\ \ldots,\
-1\\}$, noting that this results in a displacement of $\Delta$ in the negative
direction for agents located at sites $i\leq k$. New lattice sites initially
contain zero agents. For negative domain growth (i.e. contraction or
shrinkage), when $\lfloor L(j\delta t)/\Delta\rfloor<\lfloor L((j-1)\delta
t)/\Delta\rfloor$, we randomly select a lattice site $k\in\\{1,\ \ldots,\
\lfloor L((j-1)\delta t)/\Delta\rfloor\\}$. We remove the lattice site at
$x=k\Delta$, and hence all agents located at sites $i>k$ are displaced a
distance $\Delta$ in the negative direction [27]. Agents located at the
removed site $k$ are not displaced. This process is repeated for
$k\in\\{-\lfloor L((j-1)\delta t)/\Delta\rfloor,\ \ldots,\ -1\\}$, where there
is a displacement of $\Delta$ in the positive direction for agents located at
sites $i<k$. A schematic of movement and domain growth events is presented in
Figure 1. We choose domain growth functions such that $L(j\delta t)=0$ does
not occur.

We next consider multiple growing domains, where the growth within each
individual domain is spatially uniform. There is a series of adjacent domains
(Figure 1), each of which grows at a (potentially distinct) uniform rate.
Accordingly, the growth rate of the entire domain can be non-uniform. Here
agents randomly move in the positive and negative directions with probability
$P_{g}/2$ where $g\in\\{1,\ \ldots,\ G\\}$ and there are $2G-1$ adjacent, non-
overlapping domains in total. Here we always examine domains that are
symmetric around $x=0$ (Figure 1(e)). We now have $L(j\delta
t)=\sum_{g=1}^{G}L_{g}(j\delta t)$. New lattice sites are inserted into
individual domains at random when $\lfloor L_{g}(j\delta
t)/\Delta\rfloor>\lfloor L_{g}((j-1)\delta t)/\Delta\rfloor$ and are removed
from the individual domains at random when $\lfloor L_{g}(j\delta
t)/\Delta\rfloor<\lfloor L_{g}((j-1)\delta t)/\Delta\rfloor$ according to the
processes described above.

We calculate a number of key statistics from individual realisations of the
random walk, and perform numerous realisations of the random walk to generate
representative average behaviour. Specifically, we calculate the average agent
density

$\overline{C}(i\Delta,j\delta t)=\frac{1}{\Delta MN}\sum_{m=1}^{M}N_{m}(i,j),\
i\in\bigg{\\{}-\sum_{g=1}^{G}\bigg{\lfloor}\frac{L_{g}(j\delta
t)}{\Delta}\bigg{\rfloor},\ \ldots,\
\sum_{g=1}^{G}\bigg{\lfloor}\frac{L_{g}(j\delta
t)}{\Delta}\bigg{\rfloor}\bigg{\\}},\ j\in\mathcal{N}_{0},$

where $M$ is the number of identically-prepared realisations of the random
walk performed, $N$ is the initial number of agents in the random walk, and
$N_{m}(i,j)$ is the number of agents located at site $i$ after $j$ timesteps.
We again note that agents are removed from the domain at the end of a timestep
if they are located at the “ghost” sites at
$x=\pm\Delta(1+\sum_{g=1}^{G}\lfloor L_{g}(j\delta t)/\Delta\rfloor)$. We
calculate the average survival probability, which represents the proportion of
the initial agents that remain on the domain after $j$ timesteps,

$\overline{S}(j\delta t)=1-\frac{1}{MN}\sum_{m=1}^{M}\sum_{i}N_{m}(i,j),\
i\in\bigg{\\{}-\sum_{g=1}^{G}\bigg{\lfloor}\frac{L_{g}(j\delta
t)}{\Delta}\bigg{\rfloor},\ \ldots,\
\sum_{g=1}^{G}\bigg{\lfloor}\frac{L_{g}(j\delta
t)}{\Delta}\bigg{\rfloor}\bigg{\\}},\ j\in\mathcal{N}_{0}.$

From these time-dependent statistics we can calculate relevant long time
statistics such as the average splitting probabilities, which are the
proportion of the agents that ever cross the negative and positive boundaries

$\overline{S}^{-}(j\delta
t)=\frac{1}{MN}\sum_{k=1}^{j}\sum_{m=1}^{M}N_{m}\bigg{(}-1-\sum_{g=1}^{G}\bigg{\lfloor}\frac{L_{g}(k\delta
t)}{\Delta}\bigg{\rfloor},k\bigg{)},$

and

$\overline{S}^{+}(j\delta
t)=\frac{1}{MN}\sum_{k=1}^{j}\sum_{m=1}^{M}N_{m}\bigg{(}1+\sum_{g=1}^{G}\bigg{\lfloor}\frac{L_{g}(k\delta
t)}{\Delta}\bigg{\rfloor},k\bigg{)}.$

Note that here, due to the growing domain, there is no guarantee that all
agents will leave the domain as $j\to\infty$. This is in contrast to the non-
growing case where all agents eventually leave the domain. The average
fraction of agents remaining on the domain at steady state is

$\overline{\psi}=\lim_{j\to\infty}\overline{S}(j\delta
t)=1-\overline{\theta}^{-}-\overline{\theta}^{+}=\lim_{j\to\infty}\Big{[}1-\overline{S}^{-}(j\delta
t)-\overline{S}^{+}(j\delta t)\Big{]}$

In practice, we must eventually terminate the simulation of the random walk,
and hence $j\in\\{0,\ \ldots,\ j_{\text{max}}\\}$ where $j_{\text{max}}$ is
chosen to be sufficiently large that the long-time observations approximate
the $t\to\infty$ steady state limit.

###  2.2 Continuum model

Previous studies have derived the relationship between a velocity field
$v(x,t)$ that translates each point $0\leq x\leq L(t)$ in a domain, and the
overall growth rate of the domain due to the expansion of each infinitesimal
width of space [12],

$\frac{\text{d}L(t)}{\text{d}t}=\int_{0}^{L(t)}\frac{\partial v(x,t)}{\partial
x}\ \text{d}x.$

We note that in all cases here we consider a symmetric growing domain centred
around $x=0$, that is, $v(x,t)=-v(-x,t)$. To account for the multiple growing
domains we choose a velocity field that is piecewise defined for the separate
domains. Accordingly,

$\frac{\text{d}L(t)}{\text{d}t}=\sum_{g=1}^{G}\frac{\text{d}L_{g}(t)}{\text{d}t}=\int_{0}^{L_{1}(t)}\frac{\partial
v}{\partial x}\
\text{d}x+\sum_{g=1}^{G-1}\int_{B_{g}(t)}^{B_{g+1}(t)}\frac{\partial
v}{\partial x}\ \text{d}x,$

where, for ease of notation, we define the position of the $g$th moving
boundary

$B_{g}(t)=\sum_{i=1}^{g}L_{i}(t),\ g\in\\{1,\ \ldots,\ G\\},$ (1)

with $B_{0}(t)=0$. Here we select spatially-uniform growth within each
individual domain. That is, $v(x,t)$ is piecewise and chosen such that
$\partial v/\partial x$ is independent of $x$ within an individual domain,
though we note there may be a temporal dependence,

$\frac{\partial v(x,t)}{\partial
x}=\sigma_{g}(t)=\frac{1}{L_{g}(t)}\frac{\text{d}L_{g}(t)}{\text{d}t},\
B_{g-1}(t)<x<B_{g}(t),\ g\in\\{1,\ \ldots,\ G\\}.$

The velocity field is therefore

$v(x,t)=\frac{1}{L_{g}(t)}\Big{(}x-B_{g-1}(t)\Big{)}\frac{\text{d}L_{g}(t)}{\text{d}t}+\sum_{i=1}^{g-1}\frac{\text{d}L_{i}(t)}{\text{d}t},\
\ B_{g-1}(t)<x<B_{g}(t),\ g\in\\{1,\ \ldots,\ G\\},$ (2)

and recalling that $v(x,t)=-v(-x,t)$ with $v(0,t)=0$. We observe that the
velocity field at the boundary between domains is the sum of the domain growth
rates for all domains to that point (Figure 1(e)).

We now consider the evolution of a mass that undergoes linear diffusion on
this growing domain. The governing equation for the mass density $C(x,t)$ is
[19]

$\frac{\partial C(x,t)}{\partial t}=\frac{\partial}{\partial
x}\bigg{(}D(x)\frac{\partial C(x,t)}{\partial
x}\bigg{)}-\frac{\partial}{\partial x}\bigg{(}v(x,t)C(x,t)\bigg{)},\
-L(t)<x<L(t),$ (3)

where the diffusivity $D(x)$ is piecewise defined

$D(x)=D_{g}=\lim_{\Delta,\delta t\to 0}\frac{\Delta^{2}P_{g}}{2\delta t},\ \
B_{g-1}(t)<x<B_{g}(t)\cup-B_{g}(t)<x<-B_{g-1}(t),\ g\in\\{1,\ \ldots,\ G\\},$
(4)

on each individual domain in line with the agent-based model, and $v(x,t)$ is
defined as above. To be consistent with the agent-based model we impose

$C(-L(t),t)=C(L(t),t)=0,$ (5)

that is, any mass that reaches the external boundary is removed from the
system. Other boundary conditions could be considered, such as a no-flux
boundary condition at one domain boundary as in [19], but such extensions are
left for future work. For the boundaries between internal domains, the flux
across the boundaries must be conserved and hence

$-D_{g}\frac{\partial C}{\partial
x}\bigg{|}_{x=B_{g}(t)^{-}}=-D_{(g+1)}\frac{\partial C}{\partial
x}\bigg{|}_{x=B_{g}(t)^{+}}.$ (6)

Equations (1)-(6) therefore define the general form of the model for a yet-to-
be-specified choice of $L(t)$. Here all functional forms for $L(t)$ are chosen
to be smooth functions where $L_{g}(t)=0$ does not occur, thereby avoiding
finite-time blow-up in the solutions.

####  2.2.1 Fixed domain solutions

To aid in constructing exact solutions to the model for multiple growing
domains, we first focus on the simplest case where there is a single domain of
fixed size. Here the governing equation reduces to the classical diffusion
equation on a fixed domain

$\frac{\partial C(x,t)}{\partial t}=D\frac{\partial^{2}C(x,t)}{\partial
x^{2}},\ -L<x<L,\ C(-L,t)=C(L,t)=0.$

For a Dirac delta initial condition located at $x=x_{0}$, the solution
obtained via method of images [9] is

$\displaystyle C(x,t)=$ $\displaystyle\frac{1}{\sqrt{4\pi
Dt}}\exp\bigg{(}-\frac{(x-x_{0})^{2}}{4Dt}\bigg{)}$
$\displaystyle+\frac{1}{\sqrt{4\pi
Dt}}\Bigg{[}\sum_{j=1}^{\infty}(-1)^{j}\exp\bigg{(}-\frac{(x-(-2jL+(-1)^{j}x_{0}))^{2}}{4Dt}\bigg{)}$
$\displaystyle+(-1)^{j}\exp\bigg{(}-\frac{(x-(2jL+(-1)^{j}x_{0}))^{2}}{4Dt}\bigg{)}\Bigg{]}.$
(7)

The fluxes at the two boundaries, which represent the instantaneous rate of
mass leaving the domain, are

$\displaystyle F^{-}(t)=-D\frac{\partial C(x,t)}{\partial x}\bigg{|}_{x=-L}$
$\displaystyle=\sum_{j=0}^{\infty}(-1)^{j}\frac{(2j+1)L+(-1)^{j}x_{0}}{\sqrt{4\pi
Dt^{3}}}\exp\bigg{(}-\frac{((2j+1)L+(-1)^{j}x_{0})^{2}}{4Dt}\bigg{)},$ (8)
$\displaystyle F^{+}(t)=-D\frac{\partial C(x,t)}{\partial x}\bigg{|}_{x=L}$
$\displaystyle=\sum_{j=0}^{\infty}(-1)^{j}\frac{(2j+1)L-(-1)^{j}x_{0}}{\sqrt{4\pi
Dt^{3}}}\exp\bigg{(}-\frac{((2j+1)L-(-1)^{j}x_{0})^{2}}{4Dt}\bigg{)}.$ (9)

Note that while the solution is an infinite summation, in practice, we can
approximate this solution via a truncated summation. In all cases we choose a
summation with 10 terms and verify that the solution is not sensitive to the
inclusion of additional terms.

####  2.2.2 Single growing domain solutions

We next focus on exact solutions for the case with a single growing domain.
Here the governing equation is [19]

$\frac{\partial C(x,t)}{\partial t}=D\frac{\partial^{2}C(x,t)}{\partial
x^{2}}-\frac{\partial}{\partial x}\bigg{(}v(x,t)C(x,t)\bigg{)},\
-L(t)<x<L(t),\ C(L(t),t)=C(-L(t),t)=0,$

where $L(t)$ is yet to be specified. If we can transform this equation to a
linear diffusion equation on a fixed domain, we can exploit the various
properties of solutions to the linear diffusion equation. We consider
transformations of $C(x,t),\ x$, and $t$ to $C(\xi,T),\ \xi$, and $T$,
respectively. A natural spatial scale is the ratio of the domain size to the
spatial variable and hence we introduce

$\xi=\frac{x}{L(t)}L(0),\ -L(0)<\xi<L(0).$

Recalling that for spatially-uniform growth

$v(x,t)=\frac{x}{L(t)}\frac{\text{d}L}{\text{d}t},$

we obtain, after cancellation of the first order terms in $\xi$,

$\frac{\partial C(\xi,t)}{\partial
t}=\frac{DL(0)^{2}}{L(t)^{2}}\frac{\partial^{2}C(\xi,t)}{\partial\xi^{2}}-\frac{C(\xi,t)}{L(t)}\frac{\text{d}L(t)}{\text{d}t}.$

The new source term acts to dilute (concentrate) the mass density in the
transformed coordinate system, reflecting the stretched (contracted) nature of
the domain for $\text{d}L(t)/\text{d}t>0$ ($\text{d}L(t)/\text{d}t<0$). We
next introduce a temporal scaling such that the coefficient of the second
$\xi$ derivative term is a constant

$T(t)=\int_{0}^{t}\bigg{(}\frac{L(0)}{L(s)}\bigg{)}^{2}\ \text{d}s,$

and hence

$\frac{\partial C(\xi,T)}{\partial
T}=D\frac{\partial^{2}C(\xi,T)}{\partial\xi^{2}}-f(T)C(\xi,T),$

where

$f(T)=\frac{L(t)}{L(0)^{2}}\frac{\text{d}L(t)}{\text{d}t}.$

We now define

$C(\xi,T)=U(\xi,T)\exp\bigg{(}-\int_{0}^{T}f(s)\ \text{d}s\bigg{)}.$

Via a change of variable in the integral we obtain

$C(\xi,T)=U(\xi,T)\frac{L(0)}{L(t)},$

and hence $U(\xi,T)$ satisfies

$\frac{\partial U(\xi,T)}{\partial
T}=D\frac{\partial^{2}U(\xi,T)}{\partial\xi^{2}},\ -L(0)<\xi<L(0),\ T>0,$

which is the linear diffusion equation on a finite domain. This result is
equivalent to the result derived in Simpson _et al._ [19]. As such, there are
various well-known methods of solution for $U(\xi,T)$ and it remains to
transform back from the fixed domain solution $U(\xi,T)$ to the growing domain
solution $C(x,t)$. Consider the fundamental solution to the diffusion equation
on an infinite domain with a Dirac delta initial condition at $x=x_{0}$, which
corresponds to $\xi=x_{0}L(0)/L(t)$,

$U(\xi,T)=\frac{1}{\sqrt{4\pi DT}}\exp\bigg{(}-\frac{\Big{(}\xi-
x_{0}\frac{L(0)}{L(t)}\Big{)}^{2}}{4DT}\bigg{)},\ -\infty<\xi<\infty,$

and hence

$C(\xi,T)=\frac{1}{\sqrt{4\pi DT}}\exp\bigg{(}-\frac{\Big{(}\xi-
x_{0}\frac{L(0)}{L(t)}\Big{)}^{2}}{4DT}\bigg{)}\frac{L(0)}{L(t)}.$

Transforming back into the original space and time coordinates, we have

$C(x,t)=\frac{1}{\sqrt{4\pi
DT(t)}}\exp\bigg{(}-\frac{\Big{(}x\frac{L(0)}{L(t)}-x_{0}\Big{)}^{2}}{4DT(t)}\bigg{)}\frac{L(0)}{L(t)}.$

To satisfy the homogeneous Dirichlet boundary conditions at $x=\pm L(t)$,
additional terms are required in the solution to account for the flux through
the boundaries and hence

$\displaystyle C(x,t)=$ $\displaystyle\frac{1}{\sqrt{4\pi
DT(t)}}\frac{L(0)}{L(t)}\Bigg{[}\exp\bigg{(}-\frac{\Big{(}x\frac{L(0}{L(t)}-x_{0}\Big{)}^{2}}{4DT(t)}\bigg{)}$
$\displaystyle+\Bigg{[}\sum_{j=1}^{\infty}(-1)^{j}\exp\bigg{(}-\frac{\Big{(}x\frac{L(0)}{L(t)}-(-2jL(0)+(-1)^{j}x_{0})\Big{)}^{2}}{4DT(t)}\bigg{)}$
$\displaystyle+(-1)^{j}\exp\bigg{(}-\frac{\Big{(}x\frac{L(0)}{L(t)}-(2jL(0)+(-1)^{j}x_{0})\Big{)}^{2}}{4DT(t)}\bigg{)}\Bigg{]}.$
(10)

####  2.2.3 Multiple fixed domain solutions

Figure 2: (a) Illustration of a solution to the diffusion equation on multiple
fixed domains and (b) its constituent components. The multiple solution
profiles centred at the boundaries are an illustrative representation of the
convolution terms in Equations (11), (13) and (14).

Before considering the solutions for the model of multiple growing domains, we
now consider exact solutions for the simplified case of multiple fixed domains
(Figure 2). This allows us to demonstrate the methods for accounting for the
flux across internal boundaries, which will be employed to solve the equations
for the multiple growing domains case. Here the governing equations are

$\displaystyle\frac{\partial C_{1}(x,t)}{\partial t}$
$\displaystyle=D_{1}\frac{\partial}{\partial x}\bigg{(}\frac{\partial
C_{1}(x,t)}{\partial x}\bigg{)},\ -L_{1}<x<L_{1},$
$\displaystyle\frac{\partial C_{2}(x,t)}{\partial t}$
$\displaystyle=D_{2}\frac{\partial}{\partial x}\bigg{(}\frac{\partial
C_{2}(x,t)}{\partial x}\bigg{)},\ -L_{2}-L_{1}<x<-L_{1},$
$\displaystyle\frac{\partial C_{3}(x,t)}{\partial t}$
$\displaystyle=D_{2}\frac{\partial}{\partial x}\bigg{(}\frac{\partial
C_{3}(x,t)}{\partial x}\bigg{)},\ L_{1}<x<L_{1}+L_{2},$

where the rate of diffusion in domain two is always the same as in domain
three. The flux across the internal boundaries at $x=\pm L_{1}$ is conserved
and hence

$-D_{1}\frac{\partial C_{1}}{\partial
x}\bigg{|}_{x=-L_{1}}=-D_{2}\frac{\partial C_{2}}{\partial
x}\bigg{|}_{x=-L_{1}},\qquad-D_{1}\frac{\partial C_{1}}{\partial
x}\bigg{|}_{x=L_{1}}=-D_{2}\frac{\partial C_{3}}{\partial
x}\bigg{|}_{x=L_{1}}.$

Individuals are absorbed at the far boundaries, such that
$C_{2}(-L_{2}-L_{1},t)=0$ and $C_{3}(L_{1}+L_{2},t)=0$. We construct solutions
by first solving for $C_{1}(x,t)$ subject to $C_{1}(\pm L_{1},t)=0$. We will
denote this solution $G(x,t)$, which is the solution to a related problem
(with different boundary conditions) that provides the insight necessary to
solve the problem for the original boundary conditions. We calculate the flux
that leaves the inner domain according to $G(x,t)$, and then re-introduce a
fraction of that flux at $x=\pm L_{1}$. This fraction represents the agents
that cross the boundary but ultimately return to the original domain. Due to
the linearity property of the diffusion equation, we are able to superimpose
solutions that satisfy the governing equations. We have defined the boundary
fluxes $F^{-}(t)$ and $F^{+}(t)$ in Equations (8)-(9). The relevant fraction
to be added at the boundary can be determined from conservation of flux.
Observing that the solution to $C_{1}(x,t)$ is composed of the related
solution $G(x,t)$ (which, by definition, is zero at $x=\pm L_{1}$) and the re-
introduced sources, the solution at the boundary at $x=L_{1}$ is

$C_{1}(L_{1},t)=\lim_{s\to t}\int_{0}^{s}\frac{w_{1}^{+}}{\sqrt{\pi
D_{1}(t-\tau)}}F^{+}(\tau)\ \text{d}\tau,$

to first order; see below for a discussion on the timescales over which this
solution is accurate. Here $w_{1}^{+}$ is both the fraction of the flux that
remains in domain one, and a weight that relates the flux at $x=L_{1}$ in the
positive $x$ direction for $G(x,t)$ to the flux for $C_{1}(x,t)$, that is,

$-D_{1}(1-w_{1}^{+})\frac{\partial G(x,t)}{\partial
x}\bigg{|}_{x=L_{1}}=-D_{1}\frac{\partial C_{1}(x,t)}{\partial
x}\bigg{|}_{x=L_{1}}.$

The solution in domain three at the same boundary, which is only composed to
the introduced sources, is (to first order)

$C_{3}(L_{1},t)=\lim_{s\to t}\int_{0}^{s}\frac{w_{3}^{-}}{\sqrt{\pi
D_{2}(s-\tau)}}F^{+}(\tau)\ \text{d}\tau,$

where $w_{3}^{-}$ is the fraction of flux that enters domain three, subject to
the restrction that $w_{1}^{+}+w_{3}^{-}=1$. For conservation of flux to hold,
we require

$w_{1}^{+}=\frac{\sqrt{D_{2}}}{\sqrt{D_{1}}+\sqrt{D_{2}}},\qquad
w_{3}^{-}=\frac{\sqrt{D_{1}}}{\sqrt{D_{1}}+\sqrt{D_{2}}}$

at $x=L_{1}$. Similarly, for the boundary at $x=-L_{1}$, we require

$w_{1}^{-}=\frac{\sqrt{D_{2}}}{\sqrt{D_{1}}+\sqrt{D_{2}}},\qquad
w_{2}^{+}=\frac{\sqrt{D_{1}}}{\sqrt{D_{1}}+\sqrt{D_{2}}}.$

If $D_{1}=D_{2}$ then it is equally likely for an agent that crosses the
boundary to be located in either domain. If $D_{1}>D_{2}$ then agents reside
in the outer domains for longer than in the inner domain, and the agent
density will be higher in the outer domains; equally, if $D_{1}<D_{2}$ then
agents reside in the inner domain for longer, and the agent density will be
higher in the inner domain. For a Dirac delta initial condition located at
$x=x_{0}$ the solution on $-L_{1}<x<L_{1}$ is therefore (Figure 2)

$\displaystyle
C_{1}(x,t)=G(x,t)+\int_{0}^{t}w_{1}^{-}F^{-}(\tau)P_{1}^{-}(x,-L_{1},t-\tau)+w_{1}^{+}F^{+}(\tau)P_{1}^{+}(x,L_{1},t-\tau)\
\text{d}\tau,$ (11)

where

$\displaystyle G(x,t)=$ $\displaystyle\frac{1}{\sqrt{4\pi
D_{1}t}}\exp\bigg{(}-\frac{(x-x_{0})^{2}}{4D_{1}t}\bigg{)}$
$\displaystyle+\frac{1}{\sqrt{4\pi
D_{1}t}}\Bigg{[}\sum_{j=1}^{\infty}(-1)^{j}\exp\bigg{(}-\frac{(x-(-2jL_{1}+(-1)^{j}x_{0}))^{2}}{4D_{1}t}\bigg{)}$
$\displaystyle+(-1)^{j}\exp\bigg{(}-\frac{(x-(2jL_{1}+(-1)^{j}x_{0}))^{2}}{4D_{1}t}\bigg{)}\Bigg{]},$
(12)

is the solution to the diffusion equation with homogeneous Dirichlet boundary
conditions, $P_{i}^{k}(x,x_{0},t),\ k\in\\{-,+\\}$ is the fundamental solution
to the diffusion equation on domain $i$ with a Dirac delta initial condition
at $x=x_{0}$, a homogeneous Neumann boundary condition at the leftmost
boundary, and a homogeneous Dirichlet boundary condition at the rightmost
boundary (for $k=-$; the boundary conditions are reversed for $k=+$), and

$F^{-}(t)=-D_{1}\frac{\partial G}{\partial x}\bigg{|}_{x=-L_{1}},\qquad
F^{+}(t)=-D_{1}\frac{\partial G}{\partial x}\bigg{|}_{x=L_{1}},$

as before. For example,

$\displaystyle P_{1}^{-}(x,x_{0},t)=$ $\displaystyle\frac{1}{\sqrt{4\pi
D_{1}t}}\Bigg{[}\exp\bigg{(}-\frac{(x-x_{0})^{2}}{4D_{1}t}\bigg{)}$
$\displaystyle+\sum_{j=1}^{\infty}(-1)^{\lfloor(j+1)/2\rfloor}\exp\bigg{(}-\frac{(x-(-2jL_{1}+(-1)^{j}x_{0})^{2}}{4D_{1}t}\bigg{)}$
$\displaystyle+\sum_{j=1}^{\infty}(-1)^{\lfloor
j/2\rfloor}\exp\bigg{(}-\frac{(x-(2jL_{1}+(-1)^{j}x_{0})^{2}}{4D_{1}t}\bigg{)}\Bigg{]}.$

Note that the mass contained in the solution components corresponding to the
reintroduced flux will eventually reach the opposing boundary (i.e. mass
reintroduced at $x=L_{1}$ will reach $x=-L_{1}$). Due to the homogeneous
Dirichlet boundary imposed at $x=-L_{1}$ for the solution component
$P_{1}^{+}(x,L_{1},t-\tau)$, this implies that there will be an additional
flux through the boundary at $x=-L_{1}$. This could be accounted for in a
similar way to the original boundary flux, where the appropriate fraction of
the flux is reintroduced on either side of the boundary via a convolution
term. The solution would therefore involve a double convolution component.
However, here we consider problems where the average timescale for an
individual to cross both internal boundaries is beyond the timescale of
interest, so we neglect this additional term. The solutions for
$-L_{1}-L_{2}<x<-L_{1}$ and $L_{1}<x<L_{1}+L_{2}$ are

$\displaystyle C_{2}(x,t)$
$\displaystyle=\int_{0}^{t}w_{2}^{+}F^{-}(\tau)P_{2}^{+}(x,-L_{1},t-\tau)\
\text{d}\tau,$ (13) $\displaystyle C_{3}(x,t)$
$\displaystyle=\int_{0}^{t}w_{3}^{-}F^{+}(\tau)P_{3}^{-}(x,L_{1},t-\tau)\
\text{d}\tau.$ (14)

The solution at the boundary exhibits the jump condition [30]

$\frac{C_{1}(L,t)}{C_{3}(L,t)}=\frac{D_{2}}{D_{1}}.$

####  2.2.4 Multiple growing domain solutions

We next focus on obtaining exact solutions for the case with multiple growing
domains. For $G=2$, the governing equations are

$\displaystyle\frac{\partial C_{1}(x,t)}{\partial t}$
$\displaystyle=D_{1}\frac{\partial}{\partial x}\bigg{(}\frac{\partial
C_{1}(x,t)}{\partial x}\bigg{)}-\frac{\partial}{\partial
x}\bigg{(}v_{1}(x,t)C_{1}(x,t)\bigg{)},\ -L_{1}(t)<x<L_{1}(t),$
$\displaystyle\frac{\partial C_{2}(x,t)}{\partial t}$
$\displaystyle=D_{2}\frac{\partial}{\partial x}\bigg{(}\frac{\partial
C_{2}(x,t)}{\partial x}\bigg{)}-\frac{\partial}{\partial
x}\bigg{(}v_{2}(x,t)C_{2}(x,t)\bigg{)},\ -L_{2}(t)-L_{1}(t)<x<-L_{1}(t),$
$\displaystyle\frac{\partial C_{3}(x,t)}{\partial t}$
$\displaystyle=D_{2}\frac{\partial}{\partial x}\bigg{(}\frac{\partial
C_{3}(x,t)}{\partial x}\bigg{)}-\frac{\partial}{\partial
x}\bigg{(}v_{3}(x,t)C_{3}(x,t)\bigg{)},\ L_{1}(t)<x<L_{1}(t)+L_{2}(t)$

where $L_{1}(t)$ and $L_{2}(t)$ are yet to be specified and
$v_{3}(x,t)=-v_{2}(-x,t)$ such that the domain is symmetric. We consider
transformations of $C_{1}(x,t)$, $C_{2}(x,t)$, $C_{3}(x,t)$, $x$, and $t$ to
$C_{1}(\xi_{1},T_{1})$, $C_{2}(\xi_{2},T_{2})$, $C_{3}(\xi_{3},T_{2})$,
$\xi_{1}$, $\xi_{2}$, $\xi_{3}$, $T_{1}$, and $T_{2}$. The three spatial
transformations are

$\xi_{1}=\bigg{(}\frac{x}{L_{1}(t)}\bigg{)}L_{1}(0),\qquad\xi_{2}=\bigg{(}\frac{x-L_{1}(t)}{L_{2}(t)}\bigg{)}L_{2}(0)+L_{1}(0),\qquad\xi_{3}=\bigg{(}\frac{x+L_{1}(t)}{L_{2}(t)}\bigg{)}L_{2}(0)-L_{1}(0).$

Here we have

$\displaystyle
v_{1}(x,t)=\bigg{(}\frac{x}{L_{1}(t)}\bigg{)}\frac{\text{d}L_{1}}{\text{d}t},$
$\displaystyle\qquad
v_{2}(x,t)=\bigg{(}\frac{x-L_{1}(t)}{L_{2}(t)}\bigg{)}\frac{\text{d}L_{2}}{\text{d}t}+\frac{\text{d}L_{1}}{\text{d}t},$
$\displaystyle\qquad v_{3}(x,t)$
$\displaystyle=\bigg{(}\frac{x+L_{1}(t)}{L_{2}(t)}\bigg{)}\frac{\text{d}L_{2}}{\text{d}t}-\frac{\text{d}L_{1}}{\text{d}t}.$

Applying the transformations, we obtain

$\displaystyle\frac{\partial C_{1}(\xi_{1},t)}{\partial t}$
$\displaystyle=D_{1}\bigg{(}\frac{L_{1}(0)}{L_{1}(t)}\bigg{)}^{2}\frac{\partial^{2}C_{1}(\xi_{1},t)}{\partial\xi_{1}^{2}}-\frac{1}{L_{1}(t)}\frac{\text{d}L_{1}}{\text{d}t}C_{1}(\xi_{1},t),\
-L_{1}(0)<\xi_{1}<L_{1}(0),$ $\displaystyle\frac{\partial
C_{2}(\xi_{2},t)}{\partial t}$
$\displaystyle=D_{2}\bigg{(}\frac{L_{2}(0)}{L_{2}(t)}\bigg{)}^{2}\frac{\partial^{2}C_{2}(\xi_{2},t)}{\partial\xi_{2}^{2}}-\frac{1}{L_{2}(t)}\frac{\text{d}L_{2}}{\text{d}t}C_{2}(\xi_{2},t),\
-L_{2}(0)-L_{1}(0)<\xi_{2}<-L_{1}(0),$ $\displaystyle\frac{\partial
C_{3}(\xi_{3},t)}{\partial t}$
$\displaystyle=D_{2}\bigg{(}\frac{L_{2}(0)}{L_{2}(t)}\bigg{)}^{2}\frac{\partial^{2}C_{3}(\xi_{3},t)}{\partial\xi_{3}^{2}}-\frac{1}{L_{2}(t)}\frac{\text{d}L_{2}}{\text{d}t}C_{3}(\xi_{3},t),\
L_{1}(0)<\xi_{3}<L_{1}(0)+L_{2}(0).$

The time transformations are

$T_{1}(t)=\int_{0}^{t}\bigg{(}\frac{L_{1}(0)}{L_{1}(s)}\bigg{)}^{2}\
\text{d}s,\qquad
T_{2}(t)=\int_{0}^{t}\bigg{(}\frac{L_{2}(0)}{L_{2}(s)}\bigg{)}^{2}\
\text{d}s,$

which give

$\displaystyle\frac{\partial C_{1}(\xi_{1},T_{1})}{\partial T_{1}}$
$\displaystyle=D_{1}\frac{\partial^{2}C_{1}(\xi_{1},T_{1})}{\partial\xi_{1}^{2}}-\frac{L_{1}(t)}{L_{1}(0)^{2}}\frac{\text{d}L_{1}}{\text{d}t}C_{1}(\xi_{1},T_{1}),\
-L_{1}(0)<\xi_{1}<L_{1}(0),$ $\displaystyle\frac{\partial
C_{2}(\xi_{2},T_{2})}{\partial T_{2}}$
$\displaystyle=D_{2}\frac{\partial^{2}C_{2}(\xi_{2},T_{2})}{\partial\xi_{2}^{2}}-\frac{L_{2}(t)}{L_{2}(0)^{2}}\frac{\text{d}L_{2}}{\text{d}t}C_{2}(\xi_{2},T_{2}),\
-L_{2}(0)-L_{1}(0)<\xi_{2}<-L_{1}(0),$ $\displaystyle\frac{\partial
C_{3}(\xi_{3},T_{2})}{\partial T_{2}}$
$\displaystyle=D_{2}\frac{\partial^{2}C_{3}(\xi_{3},T_{2})}{\partial\xi_{3}^{2}}-\frac{L_{2}(t)}{L_{2}(0)^{2}}\frac{\text{d}L_{2}}{\text{d}t}C_{3}(\xi_{3},T_{2}),\
L_{1}(0)<\xi_{3}<L_{1}(0)+L_{2}(0).$

If we define

$f_{1}(T_{1})=\frac{L_{1}(t)}{L_{1}(0)^{2}}\frac{\text{d}L_{1}}{\text{d}t},\qquad
f_{2}(T_{2})=\frac{L_{2}(t)}{L_{2}(0)^{2}}\frac{\text{d}L_{2}}{\text{d}t},$

and

$C_{1}(\xi_{1},T_{1})=U_{1}(\xi_{1},T_{1})\frac{L_{1}(0)}{L_{1}(t)},\qquad
C_{2}(\xi_{2},T_{2})=U_{2}(\xi_{2},T_{2})\frac{L_{2}(0)}{L_{2}(t)},\qquad
C_{3}(\xi_{3},T_{2})=U_{3}(\xi_{3},T_{2})\frac{L_{2}(0)}{L_{2}(t)},$

it can be seen that $U_{1}(\xi_{1},T_{1})$, $U_{2}(\xi_{2},T_{2})$ and
$U_{3}(\xi_{3},T_{2})$ satisfy

$\displaystyle\frac{\partial U_{1}(\xi_{1},T_{1})}{\partial T_{1}}$
$\displaystyle=D_{1}\frac{\partial^{2}U_{1}(\xi_{1},T_{1})}{\partial\xi_{1}^{2}},\
-L_{1}(0)<\xi_{1}<L_{1}(0),$ $\displaystyle\frac{\partial
U_{2}(\xi_{2},T_{2})}{\partial T_{2}}$
$\displaystyle=D_{2}\frac{\partial^{2}U_{2}(\xi_{2},T_{2})}{\partial\xi_{2}^{2}},\
-L_{2}(0)-L_{1}(0)<\xi_{2}<-L_{1}(0),$ $\displaystyle\frac{\partial
U_{3}(\xi_{3},T_{2})}{\partial T_{2}}$
$\displaystyle=D_{2}\frac{\partial^{2}U_{3}(\xi_{3},T_{2})}{\partial\xi_{3}^{2}},\
L_{1}(0)<\xi_{3}<L_{1}(0)+L_{2}(0),$

which is the linear diffusion equation on multiple fixed domains. As for the
previous multiple domain solutions, we impose the conservation of diffusive
flux across the internal boundaries, noting that the velocity fields at the
boundaries are continuous. We assume that there is a Dirac delta initial
condition at $x=x_{0}$ where $-L_{1}(0)<x_{0}<L_{1}(0)$ and hence $C_{1}(x,t)$
will have components arising from the initial mass and the re-introduced
boundary flux, as before, while $C_{2}(x,t)$ and $C_{3}(x,t)$ will only have
components from the internal boundary flux. Due to the growing domains, the
weights for the boundary flux terms now depend on the transformations of $t$
and $C_{i}(x,t)$. We again require the weights sum to one and that, if
$D_{1}=D_{2}$, then $C_{1}(L_{1},t)=C_{3}(L_{1},t)$ due to the continuity of
the velocity fields. The weights are

$\displaystyle w_{1}^{+}(t,\tau)$
$\displaystyle=\frac{\sqrt{D_{2}}\sqrt{T_{2}(t,\tau)}L_{2}(t)/L_{2}(\tau)}{\sqrt{D_{2}}\sqrt{T_{2}(t,\tau)}L_{2}(t)/L_{2}(\tau)+\sqrt{D_{1}}\sqrt{T_{1}(t,\tau)}L_{1}(t)/L_{1}(\tau)},$
$\displaystyle w_{3}^{-}(t,\tau)$
$\displaystyle=\frac{\sqrt{D_{1}}\sqrt{T_{1}(t,\tau)}L_{1}(t)/L_{1}(\tau)}{\sqrt{D_{2}}\sqrt{T_{2}(t,\tau)}L_{2}(t)/L_{2}(\tau)+\sqrt{D_{1}}\sqrt{T_{1}(t,\tau)}L_{1}(t)/L_{1}(\tau)},$

for the boundary at $x=L_{1}(t)$, and

$\displaystyle w_{1}^{-}(t,\tau)$
$\displaystyle=\frac{\sqrt{D_{2}}\sqrt{T_{2}(t,\tau)}L_{2}(t)/L_{2}(\tau)}{\sqrt{D_{2}}\sqrt{T_{2}(t,\tau)}L_{2}(t)/L_{2}(\tau)+\sqrt{D_{1}}\sqrt{T_{1}(t,\tau)}L_{1}(t)/L_{1}(\tau)},$
$\displaystyle w_{2}^{+}(t,\tau)$
$\displaystyle=\frac{\sqrt{D_{1}}\sqrt{T_{1}(t,\tau)}L_{1}(t)/L_{1}(\tau)}{\sqrt{D_{2}}\sqrt{T_{2}(t,\tau)}L_{2}(t)/L_{2}(\tau)+\sqrt{D_{1}}\sqrt{T_{1}(t,\tau)}L_{1}(t)/L_{1}(\tau)},$

for the boundary at $x=-L_{1}(t)$, where

$T_{1}(t,\tau)=\int_{\tau}^{t}\bigg{(}\frac{L_{1}(\tau)}{L_{1}(s)}\bigg{)}^{2}\
\text{d}s,\qquad
T_{2}(t,\tau)=\int_{\tau}^{t}\bigg{(}\frac{L_{2}(\tau)}{L_{2}(s)}\bigg{)}^{2}\
\text{d}s.$

The solutions are therefore

$C_{1}(x,t)=H(x,t)+\int_{0}^{t}w_{1}^{-}(t,\tau)F_{1}^{-}(\tau)P_{1}^{-}(x,-L_{1}(t),t-\tau)+w_{1}^{+}(t,\tau)F_{1}^{+}(\tau)P_{1}^{+}(x,L_{1}(t),t-\tau)\
\text{d}\tau,$ (15)

where

$\displaystyle H(x,t)=$ $\displaystyle\frac{1}{\sqrt{4\pi
D_{1}T_{1}(t)}}\frac{L_{1}(0)}{L_{1}(t)}\Bigg{[}\exp\bigg{(}-\frac{\Big{(}x\frac{L_{1}(0)}{L_{1}(t)}-x_{0}\Big{)}^{2}}{4D_{1}T_{1}(t)}\bigg{)}$
$\displaystyle+\Bigg{[}\sum_{j=1}^{\infty}(-1)^{j}\exp\bigg{(}-\frac{\Big{(}x\frac{L_{1}(0)}{L_{1}(t)}-(-2jL_{1}(0)+(-1)^{j}x_{0})\Big{)}^{2}}{4D_{1}T_{1}(t)}\bigg{)}$
$\displaystyle+(-1)^{j}\exp\bigg{(}-\frac{\Big{(}x\frac{L_{1}(0)}{L_{1}(t)}-(2jL_{1}(0)+(-1)^{j}x_{0})\Big{)}^{2}}{4D_{1}T_{1}(t)}\bigg{)}\Bigg{]},$
(16)

and

$F_{1}^{-}(t)=-D_{1}\frac{\partial H}{\partial x}\bigg{|}_{x=-L_{1}(t)},\qquad
F_{1}^{+}(t)=-D_{1}\frac{\partial H}{\partial x}\bigg{|}_{x=L_{1}(t)}.$

The solution for domains two and three are

$\displaystyle C_{2}(x,t)$
$\displaystyle=\int_{0}^{t}w_{2}^{+}(t,\tau)F_{1}^{-}(\tau)P_{2}^{+}(x,-L_{1}(t),t-\tau)\
\text{d}\tau,$ (17) $\displaystyle C_{3}(x,t)$
$\displaystyle=\int_{0}^{t}w_{3}^{-}(t,\tau)F_{1}^{+}(\tau)P_{3}^{-}(x,L_{1}(t),t-\tau)\
\text{d}\tau,$ (18)

where $P_{i}^{k}(x,x_{0},t)$ is defined as previously, albeit for the $i$th
growing domain, rather than a fixed domain.

###  2.3 Splitting and survival probabilities

The proportion of the individuals that have yet to cross either external
boundary, known as the survival probability, is a key metric for diffusive
processes on growing domains [24, 19, 20]. The proportion of the individuals
that have crossed either the left or right external boundary, known as the
left and right splitting probabilities, is similarly of interest. We can
leverage the derived exact solutions to calculate these metrics, both for the
single growing domain and multiple growing domain cases. Critically, we can
obtain expressions for the time-varying splitting and survival probabilities,
rather than just in the long-time limit.

For the case with a single growing domain, it is relatively straightforward to
calculate the splitting and survival probabilities. The time-varying splitting
probabilities for the boundaries at $x=-L(t)$ and $x=L(t)$, denoted $S^{-}(t)$
and $S^{+}(t)$, respectively, are simply the time integral of the appropriate
boundary flux. The fluxes for the two boundaries are

$\displaystyle F^{-}(t)=-D_{1}\frac{\partial C(x,t)}{\partial
x}\bigg{|}_{x=-L(t)}=$
$\displaystyle\sum_{j=0}^{\infty}(-1)^{j}\frac{(2j+1)L(0)+(-1)^{j}x_{0}}{\sqrt{4\pi
D_{1}T(t)^{3}}}\bigg{(}\frac{L(0)}{L(t)}\bigg{)}^{2}$
$\displaystyle\times\exp\bigg{(}-\frac{((2j+1)L(0)+(-1)^{j}x_{0})^{2}}{4D_{1}T(t)}\bigg{)},$
$\displaystyle F^{+}(t)=-D_{1}\frac{\partial C(x,t)}{\partial
x}\bigg{|}_{x=L(t)}=$
$\displaystyle\sum_{j=0}^{\infty}(-1)^{j}\frac{(2j+1)L(0)-(-1)^{j}x_{0}}{\sqrt{4\pi
D_{1}T(t)^{3}}}\bigg{(}\frac{L(0)}{L(t)}\bigg{)}^{2}$
$\displaystyle\times\exp\bigg{(}-\frac{((2j+1)L(0)-(-1)^{j}x_{0})^{2}}{4D_{1}T(t)}\bigg{)}.$

The time-varying splitting probabilities are therefore

$\displaystyle S^{-}(t)=\int_{0}^{t}F_{1}^{-}(\tau)\ \text{d}\tau$
$\displaystyle=\sum_{j=0}^{\infty}(-1)^{j}\text{erfc}\bigg{(}\frac{(2j+1)L(0)+(-1)^{j}x_{0}}{\sqrt{4D_{1}T(t)}}\bigg{)},$
(19) $\displaystyle S^{+}(t)=\int_{0}^{t}F_{1}^{+}(\tau)\ \text{d}\tau$
$\displaystyle=\sum_{j=0}^{\infty}(-1)^{j}\text{erfc}\bigg{(}\frac{(2j+1)L(0)-(-1)^{j}x_{0}}{\sqrt{4D_{1}T(t)}}\bigg{)}.$
(20)

The survival probability, $S(t)$, is simply the proportion of the population
that has not crossed either boundary

$\displaystyle S(t)=1-S^{-}(t)-S^{+}(t)=$ $\displaystyle
1-\sum_{j=0}^{\infty}\Bigg{[}(-1)^{j}\text{erfc}\bigg{(}\frac{(2j+1)L(0)+(-1)^{j}x_{0}}{\sqrt{4D_{1}T(t)}}\bigg{)}$
$\displaystyle+(-1)^{j}\text{erfc}\bigg{(}\frac{(2j+1)L(0)-(-1)^{j}x_{0}}{\sqrt{4D_{1}T(t)}}\bigg{)}\Bigg{]}.$
(21)

For the case with multiple growing domains, the process to calculate the
splitting and survival probabilities is more complicated. For an individual to
cross an external boundary, it must first cross an internal boundary. We
therefore first calculate the internal splitting probabilities, $S_{1}^{-}(t)$
and $S_{1}^{+}(t)$. Note that this is the probability that an individual
crosses an internal boundary and does not re-enter the internal domain. This
is calculated by weighting the internal boundary flux as before and taking the
integral over time

$S_{1}^{-}(t)=\int_{0}^{t}w_{2}^{+}(t,\tau)F_{1}^{-}(\tau)\text{d}\tau,\qquad
S_{1}^{+}(t)=\int_{0}^{t}w_{3}^{-}(t,\tau)F)1^{+}(\tau)\text{d}\tau.$

A closed-form expression for the internal splitting probabilities can be
obtained if
$\sqrt{T_{1}(t,\tau)}L_{1}(t)/L_{1}(\tau)=\sqrt{T_{2}(t,\tau)}L_{2}(t)/L_{2}(\tau)$
holds for all $t$ and $\tau$, that is, that the ratio of the new domain size
to the original domain size is the same for both the internal and external
domains. Where this holds, the weights reduce to a constant, and hence

$\displaystyle S_{1}^{-}(t)$
$\displaystyle=\frac{\sqrt{D_{1}}}{\sqrt{D_{1}}+\sqrt{D_{2}}}\sum_{j=0}^{\infty}(-1)^{j}\text{erfc}\bigg{(}\frac{(2j+1)L_{1}(0)+(-1)^{j}x_{0}}{\sqrt{4D_{1}T(t)}}\bigg{)},$
$\displaystyle S_{1}^{+}(t)$
$\displaystyle=\frac{\sqrt{D_{1}}}{\sqrt{D_{1}}+\sqrt{D_{2}}}\sum_{j=0}^{\infty}(-1)^{j}\text{erfc}\bigg{(}\frac{(2j+1)L_{1}(0)-(-1)^{j}x_{0}}{\sqrt{4D_{1}T(t)}}\bigg{)}.$

If this condition does not hold, we can calculate the splitting probabilities
via numerical integration. As before, the internal survival probability can be
calculated according to $S_{1}(t)=1-S_{1}^{-}(t)-S_{1}^{+}(t)$.

To calculate the external splitting probabilities, that is, the probability
that an individual crosses the boundary at either $x=-L_{1}(t)-L_{2}(t)$ or
$x=L_{1}(t)+L_{2}(t)$, we first require an expression for the flux through the
appropriate boundary. Recall that the solutions in domains two and three
correspond to an integral over solution kernels that represent mass initially
located at $x=\pm L_{1}(t)$. The flux through the boundary at
$x=L_{1}(t)+L_{2}(t)$ at time $t$ due to mass that is placed at $x=L_{1}(t)$
at $t=\tau$ is

$\displaystyle F_{3}^{+}(t,\tau)=$ $\displaystyle-D_{2}\frac{\partial
P_{3}(x,L_{1}(t),t-\tau)}{\partial x}\bigg{|}_{x=L_{1}(t)+L_{2}(t)}$
$\displaystyle=\sum_{j=0}^{\infty}(-1)^{j}\frac{(2j+1)(L_{1}(\tau)+L_{2}(\tau))-(-1)^{j}L_{1}(\tau)}{\sqrt{\pi
D_{2}T_{2}(t,\tau)^{3}}}\bigg{(}\frac{L_{2}(\tau)}{L_{2}(t)}\bigg{)}^{2}\times$
$\displaystyle\exp\bigg{(}-\frac{\Big{(}(2j+1)(L_{1}(\tau)+L_{2}(\tau))-(-1)^{j}L_{1}(\tau)\Big{)}^{2}}{4D_{2}T_{2}(t,\tau)}\bigg{)}.$

The survival probability at time $t$ in domain three for the mass that is
placed at $x=L_{1}(t)$ at $t=\tau$ is obtained through the relationship

$1-S_{3}(t,\tau)=\int_{\tau}^{t}F_{3}^{+}(s,\tau)\
\text{d}s=\sum_{j=0}^{\infty}(-1)^{j}\text{erfc}\bigg{(}\frac{(2j+1)(L_{1}(\tau)+L_{2}(\tau))-(-1)^{j}L_{1}(\tau)}{\sqrt{4D_{2}T_{2}(t,\tau)}}\bigg{)},$

noting that there is no flux from domain three into domain one, as can be seen
from the form of the solution kernel. While agents in the simulation do cross
the internal boundaries multiple times, these are accounted for in the kernel
via the weighting of the flux. These results are for a solution kernel
corresponding to a mass that is placed in a single instant; however the
solution in domain three involves an integral over the mass placed at each
instant. Therefore, the overall splitting probability for the boundary at
$x=L_{1}(t)+L_{2}(t)$ is

$S^{+}(t)=\int_{0}^{t}F_{1}^{+}(\phi)\int_{\phi}^{t}w_{3}^{-}(s,\phi)F_{3}^{+}(s,\phi)\
\text{d}s\ \text{d}\phi.$

The three components in the integral are the flux across the internal boundary
($F_{1}^{+}(\phi)$), the fraction of that flux that is placed into domain
three ($w_{3}^{-}(s,\phi)$) and the flux across the external boundary
($F_{3}^{+}(s,\phi)$), given that there was mass placed at $x=L_{1}(\phi)$. In
the case where
$\sqrt{T_{1}(t,\tau)}L_{1}(t)/L_{1}(\tau)=\sqrt{T_{2}(t,\tau)}L_{2}(t)/L_{2}(\tau)$
we can use the expressions for the internal boundary flux and the survival
probability for domain three to obtain a single integral expression for the
splitting probability for the boundary at $x=L_{1}(t)+L_{2}(t)$

$\displaystyle S^{+}(t)=$
$\displaystyle\frac{\sqrt{D_{1}}}{\sqrt{D_{1}}+\sqrt{D_{2}}}\int_{0}^{t}\bigg{[}\sum_{j=0}^{\infty}(-1)^{j}\frac{(2j+1)L_{1}(0)-(-1)^{j}x_{0}}{\sqrt{4\pi
D_{1}T_{1}(\phi)^{3}}}\bigg{(}\frac{L_{1}(0)}{L(\phi)}\bigg{)}^{2}$
$\displaystyle\times\exp\bigg{(}-\frac{((2j+1)L_{1}(0)-(-1)^{j}x_{0})^{2}}{4D_{1}T_{1}(\phi)}\bigg{)}\bigg{]}\times$
$\displaystyle\bigg{[}\sum_{j=0}^{\infty}(-1)^{j}\text{erfc}\bigg{(}\frac{(2j+1)(L_{1}(\phi)+L_{2}(\phi))-(-1)^{j}L_{1}(\phi)}{\sqrt{4D_{2}T_{2}(t,\phi)}}\bigg{)}\bigg{]}\
\text{d}\phi.$ (22)

As before, if this relationship does not hold, we can use numerical techniques
to evaluate the integral. Following similar arguments, we obtain the splitting
probability for the boundary at $x=-L_{1}(t)-L_{2}(t)$

$S^{-}(t)=\int_{0}^{t}F_{1}^{-}(\phi)\int_{\phi}^{t}w_{2}^{+}(s,\phi)F_{2}^{-}(s,\phi)\
\text{d}s\ \text{d}\phi,$

which, provided
$\sqrt{T_{1}(t,\tau)}L_{1}(t)/L_{1}(\tau)=\sqrt{T_{2}(t,\tau)}L_{2}(t)/L_{2}(\tau)$
is satisfied, can be expressed

$\displaystyle S^{-}(t)=$
$\displaystyle\frac{\sqrt{D_{1}}}{\sqrt{D_{1}}+\sqrt{D_{2}}}\int_{0}^{t}\bigg{[}\sum_{j=0}^{\infty}(-1)^{j}\frac{(2j+1)L_{1}(0)+(-1)^{j}x_{0}}{\sqrt{4\pi
D_{1}T_{1}(\phi)^{3}}}\bigg{(}\frac{L_{1}(0)}{L(\phi)}\bigg{)}^{2}$
$\displaystyle\times\exp\bigg{(}-\frac{((2j+1)L_{1}(0)+(-1)^{j}x_{0})^{2}}{4D_{1}T_{1}(\phi)}\bigg{)}\bigg{]}\times$
$\displaystyle\bigg{[}\sum_{j=0}^{\infty}(-1)^{j}\text{erfc}\bigg{(}\frac{(2j+1)(L_{1}(\phi)+L_{2}(\phi))-(-1)^{j}L_{1}(\phi)}{\sqrt{4D_{2}T_{2}(t,\phi)}}\bigg{)}\bigg{]}\
\text{d}\phi.$ (23)

The overall survival probability (i.e. that an agent remains in any domain) is
given by

$S(t)=1-S^{-}(t)-S^{+}(t).$ (24)

##  3 Results

For all results we choose unit lattice widths and timesteps such that
$\Delta=\delta t=1$. We consider three functional forms for domain growth:

  * 1.

Linear growth. Here an individual domain undergoes (positive) growth that is
linear in time according to

$L_{i}(t)=L_{i}(0)+\beta_{i}t,$ (25)

where $\beta_{i}>0$. The corresponding time transformation is

$T_{i}(t,\tau)=\bigg{(}\frac{L_{i}(0)+\beta_{i}\tau}{L_{i}(0)+\beta_{i}t}\bigg{)}(t-\tau).$
(26)

Note that the time transformation $T_{i}(t)$ is obtained by setting $\tau=0$.
It is possible to select $\beta_{i}<0$; however, this requires the restriction
that $t<-L_{i}(0)/\beta_{i}$ to ensure that the domain size remains positive
and the solution does not experience finite-time blow-up.

  * 2.

Exponential growth. Here an individual domain undergoes growth that is
exponential in time according to

$L_{i}(t)=L_{i,\text{min}}+(L_{i}(0)-L_{i,\text{min}})\exp(-\beta_{i}t).$ (27)

For exponential growth we do not have a restriction on the sign of $\beta_{i}$
as the domain will approach the finite size $L_{i,\text{min}}>0$ when
$\beta_{i}>0$, and hence we can study examples of both positive and negative
growth. The corresponding time transformation is

$\displaystyle T_{i}(t,\tau)=$
$\displaystyle\frac{\Big{[}\exp(-\beta_{i}\tau)\big{(}L_{i}(0)-L_{i,\text{min}}\big{)}+L_{i,\text{min}}\Big{]}^{2}}{\beta_{i}L_{i,\text{min}}^{2}}$
$\displaystyle\times$
$\displaystyle\bigg{[}\big{(}L_{i}(0)-L_{i,\text{min}}\big{)}\bigg{(}\frac{1}{L_{i,\text{min}}\big{[}1-\exp(\beta_{i}\tau)]-L_{i}(0)}+\frac{1}{L_{i}(0)+L_{i,\text{min}}\big{[}\exp(\beta_{i}t)-1\big{]}}\bigg{)}$
$\displaystyle-\log\Big{(}L_{i}(0)+L_{i,\text{min}}\big{[}\exp(\beta_{i}\tau)-1\big{]}\Big{)}+\log\Big{(}L_{i}(0)+L_{i,\text{min}}\big{[}\exp(\beta_{i}t)-1\big{]}\Big{)}\bigg{]}.$
(28)

  * 3.

Oscillatory growth. Here an individual domain undergoes growth that oscillates
in time according to

$L_{i}(t)=L_{i}(0)+(L_{i}(0)-L_{i,\text{min}})\sin(\beta_{i}t).$ (29)

Again, for oscillatory growth, there is no restriction on the sign of
$\beta_{i}$ as this simply dictates whether the domain initially experiences
positive or negative growth. The domain size is bounded according to
$L_{i,\text{min}}\leq L_{i}(t)\leq 2L_{i}(0)-L_{i,\text{min}}$, where
$0<L_{i,\text{min}}\leq L_{i}(0)$. The corresponding time transformation is

$\displaystyle T_{i}(t,\tau)=$
$\displaystyle\Bigg{(}\frac{L_{i}(0)+\big{(}L_{i}(0)-L_{i,\text{min}}\big{)}\sin(\beta_{i}\tau)}{L_{i}(0)+(L_{i}(0)-L_{i,\text{min}})\sin(\beta_{i}t)}\Bigg{)}$
$\displaystyle\times\Bigg{(}\frac{1}{\beta_{i}L_{i,\text{min}}\big{[}2L_{i}(0)-L_{i,\text{min}}\big{]}\big{[}L_{i}(0)^{2}-(L_{i}(0)-L_{i,\text{min}})^{2}\big{]}^{1/2}}\Bigg{)}$
$\displaystyle\times\Bigg{[}\big{[}L_{i}(0)-L_{i,\text{min}}\big{]}\big{[}L_{i}(0)^{2}-(L_{i}(0)-L_{i,\text{min}})^{2}\big{]}^{1/2}$
$\displaystyle\times\Big{(}\cos(\beta_{i}t)\big{[}L_{i}(0)+(L_{i}(0)-L_{i,\text{min}})\sin(\beta_{i}\tau)\big{]}$
$\displaystyle-\cos(\beta_{i}\tau)\big{[}L_{i}(0)+(L_{i}(0)-L_{i,\text{min}})\sin(\beta_{i}t)\big{]}\Big{)}$
$\displaystyle+2L_{i}(0)\big{[}L_{i}(0)+(L_{i}(0)-L_{i,\text{min}})\sin(\beta_{i}\tau)\big{]}\big{[}L_{i}(0)+(L_{i}(0)-L_{i,\text{min}})\sin(\beta_{i}t)\big{]}$
$\displaystyle\times\Bigg{[}\tan^{-1}\Bigg{(}\frac{L_{i}(0)\Big{[}\tan\big{(}\beta_{i}t/2)+1\Big{]}-L_{i,\text{min}}}{\big{[}L_{i}(0)^{2}-(L_{i}(0)-L_{i,\text{min}})^{2}\big{]}^{1/2}}\Bigg{)}$
$\displaystyle-\tan^{-1}\Bigg{(}\frac{L_{i}(0)\Big{[}\tan\big{(}\beta_{i}\tau/2)+1\Big{]}-L_{i,\text{min}}}{\big{[}L_{i}(0)^{2}-(L_{i}(0)-L_{i,\text{min}})^{2}\big{]}^{1/2}}\Bigg{)}\Bigg{]}\Bigg{]}$
$\displaystyle+\Big{[}L_{i}(0)+(L_{i}(0)-L_{i,\text{min}})\sin\big{(}\beta_{i}\tau\big{)}\Big{]}^{2}\Bigg{[}\frac{2\pi
L_{i}(0)}{\beta_{i}\big{[}L_{i}(0)^{2}-(L_{i}(0)-L_{i,\text{min}})^{2}\big{]}^{3/2}}\Bigg{]}$
$\displaystyle\times\Bigg{[}\bigg{\lfloor}\frac{\beta_{i}t+\pi}{2\pi}\bigg{\rfloor}-\bigg{\lfloor}\frac{\beta_{i}\tau+\pi}{2\pi}\bigg{\rfloor}\Bigg{]}.$
(30)

To ensure that our solutions are consistent with those derived in previous
investigations [19], we compare the average individual density in the random
walk against the exact solution for a single fixed domain and a single growing
domain (Appendix A, Figures 9 and 10). As in previous investigations [19], we
observe a close match between the density profiles obtained from repeated
realisations of the random walk and the exact solutions. We also examine
whether the derived exact solutions for multiple fixed domains are consistent
with the average behaviour in the random walk. We consider both the case where
$D_{1}=D_{2}$ and where $D_{1}\neq D_{2}$ (Appendix A, Figures 11 and 12). In
both cases we observe that the exact solutions match the average random walk
behaviour well, which suggests that the derived solutions are valid.

Figure 3: Comparison between the average behaviour in the lattice-based random
walk $\overline{C}(x,t)$ (black, dashed) and the exact solutions $C_{1}(x,t)$
(cyan), $C_{2}(x,t)$ (orange) and $C_{3}(x,t)$ (pink) as defined in Equations
(15)-(18) for multiple linearly growing domains (Equations (25)-(26)) and
domain-independent diffusivities. Parameters used are $D_{1}=D_{2}=0.5$,
$L_{1}(0)=L_{2}(0)=50$, $N=1000$, $x_{0}=0$, (a) $\beta_{1}=\beta_{2}=0.01$,
(b) $\beta_{1}=\beta_{2}=0.05$, (c) $\beta_{1}=\beta_{2}=0.1$. Solution
profiles are presented at $t=200$, $t=500$ and $t=1000$. The arrow indicates
the direction of increasing time. Dashed grey lines correspond to the position
of the boundary. Average random walk behaviour is obtained from 5000
identically-prepared realisations of the random walk.

We now impose domain growth and examine whether the derived exact solutions
are valid for the case with multiple growing domains. As before, we verify the
solutions via comparison against the average random walk behaviour. We first
consider examples with positive linear domain growth for three different rates
of growth and present the results in Figure 3. In each case we observe that
the exact solutions match the average random walk behaviour in each of the
three domains. The exact solutions for $C_{1}(x,t)$, $C_{2}(x,t)$ and
$C_{3}(x,t)$ are shown in cyan, orange and pink, respectively. The location of
the internal boundary is shown via the dashed grey line, which increases in
intensity as time increases. We observe that as we increase the rate of domain
growth a smaller proportion of individuals leaves the inner domain. We note
that as $D_{1}=D_{2}$ and the ratio
$L_{1}(t)/L_{1}(0)=L_{2}(t)/L_{2}(0)=L_{3}(t)/L_{3}(0)$ is consistent between
domains this problem is comparable to the problem of domain growth on a single
domain.

Figure 4: Comparison between the average behaviour in the lattice-based random
walk $\overline{C}(x,t)$ (black, dashed) and the exact solutions $C_{1}(x,t)$
(cyan), $C_{2}(x,t)$ (orange) and $C_{3}(x,t)$ (pink) as defined in Equations
(15)-(18) for multiple linearly growing domains (Equations (25)-(26)) and
domain-independent diffusivities. Parameters used are $D_{1}=D_{2}=0.5$,
$L_{1}(0)=L_{2}(0)=50$, $N=1000$, $x_{0}=0$, (a) $\beta_{1}=0.01$,
$\beta_{2}=0.02$, (b) $\beta_{1}=0.01$, $\beta_{2}=0.04$, (c)
$\beta_{1}=0.01$, $\beta_{2}=0.08$. Solution profiles are presented at
$t=400$, $t=1000$ and $t=2000$. The arrow indicates the direction of
increasing time. Dashed grey lines correspond to the position of the boundary.
Average random walk behaviour is obtained from 5000 identically-prepared
realisations of the random walk.

If the relative rates of domain growth are different between domains, we
expect different behaviour to the single domain problem. We consider three
examples where $D_{1}=D_{2}$ and $\beta_{1}\neq\beta_{2}$ and present solution
profiles obtained from repeated realisations of the random walk and
enumerating the exact solutions for three different domain growth rate pairs
in Figure 4. In each example $\beta_{1}=0.01$, while $\beta_{2}=0.02$ in
Figure 4(a), $\beta_{2}=0.04$ in Figure 4(b) and $\beta_{2}=0.08$ in Figure
4(c). The exact solutions match the average random walk behaviour well. As
$\beta_{2}$ increases, we see fewer individuals reaching the external
boundaries, though the solution in the inner domain is largely consistent. The
additional domain growth serves to stretch the solution profile (via
advection) in the outer domains. However, this additional advection is not
sufficient for the individuals to reach the outer boundary at higher
$\beta_{2}$ values.

Figure 5: Comparison between the average behaviour in the lattice-based random
walk $\overline{C}(x,t)$ (black, dashed) and the exact solutions $C_{1}(x,t)$
(cyan), $C_{2}(x,t)$ (orange) and $C_{3}(x,t)$ (pink) as defined in Equations
(15)-(18) for multiple linearly growing domains (Equations (25)-(26)) and
domain-dependent diffusivities. Parameters used are $L_{1}(0)=L_{2}(0)=50$,
$N=1000$, $x_{0}=0$, (a) $D_{1}=0.5$, $D_{2}=0.25$, $\beta_{1}=0.01$,
$\beta_{2}=0.02$, (b) $D_{1}=0.25$, $D_{2}=0.5$, $\beta_{1}=0.01$,
$\beta_{2}=0.01$, (c) $D_{1}=0.5$, $D_{2}=0.25$, $\beta_{1}=0.02$,
$\beta_{2}=0.01$. Solution profiles are presented at (a) $t=400$, $t=1000$ and
$t=2000$, (b),(c) $t=1000$, $t=2000$ and $t=5000$. The arrow indicates the
direction of increasing time. Dashed grey lines correspond to the position of
the boundary. Average random walk behaviour is obtained from 5000 identically-
prepared realisations of the random walk.

We now consider examples where $D_{1}\neq D_{2}$. The form of the exact
solution suggests that there will be a jump discontinuity at the internal
boundary that satisfies $D_{1}C_{1}(L_{1}(t),t)=D_{2}C_{3}(L_{1}(t),t)$. For
the three examples considered, presented in Figure 5, we observe this jump
discontinuity where if $D_{1}>D_{2}$ then
$C_{1}(L_{1}(t),t)<C_{3}(L_{1}(t),t)$ and if $D_{1}<D_{2}$ then
$C_{1}(L_{1}(t),t)>C_{3}(L_{1}(t),t)$, as expected. The jump discontinuity is
accurately captured in both the exact solution and the average random walk
behaviour, and we again see that the exact solution closely describes the
average behaviour in the random walk.

Figure 6: Comparison between the average behaviour in the lattice-based random
walk $\overline{C}(x,t)$ (black, dashed) and the exact solutions $C_{1}(x,t)$
(cyan), $C_{2}(x,t)$ (orange) and $C_{3}(x,t)$ (pink) as defined in Equations
(15)-(18) for multiple exponentially growing domains (Equations (27)-(28)).
Parameters used are $D_{1}=D_{2}=0.5$, $N=1000$, $x_{0}=0$, (a)
$L_{1}(0)=L_{2}(0)=60$, $L_{1,\text{min}}=L_{2,\text{min}}=30$,
$\beta_{1}=1\times 10^{-3}$, $\beta_{2}=2\times 10^{-3}$, (b) $L_{1}(0)=40$,
$L_{2}(0)=80$, $L_{1,\text{min}}=L_{2,\text{min}}=20$, $\beta_{1}=-1\times
10^{-3}$, $\beta_{2}=2\times 10^{-3}$, (c) $L_{1}(0)=80$, $L_{2}(0)=40$,
$L_{1,\text{min}}=L_{2,\text{min}}=20$, $\beta_{1}=2\times 10^{-3}$,
$\beta_{2}=-1\times 10^{-3}$. Solution profiles are presented at $t=200$,
$t=500$ and $t=1000$. The arrow indicates the direction of increasing time.
Dashed grey lines correspond to the position of the boundary. Average random
walk behaviour is obtained from 5000 identically-prepared realisations of the
random walk.

Thus far, we have only considered examples where the domains have experienced
linear growth. We now examine three cases of exponential growth, and present
the average behaviour in the random walk and the profiles obtained from the
exact solutions in Figure 6. In Figure 6(a), both domains decrease in size. In
Figure 6(b), the inner domain experiences positive growth while the outer
domains experience negative growth. Finally, in Figure 6(c), the inner domain
undergoes negative growth and the outer domain undergoes positive growth.
Critically, we select a positive value of $L_{i,\text{min}}$ so that even if
the domain shrinks in size we do not experience finite-time blow-up in the
solution due to the collision of the solution characteristics. In each case,
there is close agreement between the average random walk behaviour and the
exact solution, which confirms that our derived solutions are appropriate for
describing negative domain growth.

Figure 7: (a),(d),(g),(j) Left splitting probability, (b),(e),(h),(k) survival
probability and (c),(f),(j),(l) right splitting probability obtained from the
exact solutions defined in Equations (22)-(24) for multiple linearly growing
domains (Equations (25)-(26)) for a range of domain growth rates and initial
locations. Parameter used are $L_{1}(0)=L_{2}(0)=50$, (a)-(c)
$D_{1}=D_{2}=0.5$, (d)-(f) $D_{1}=0.25$, $D_{2}=0.5$, (g)-(i) $D_{1}=0.5$,
$D_{2}=0.25$, (j)-(l) $D_{1}=D_{2}=0.25$. Probabilities are presented at
$t=10^{4}$.

One benefit of the exact solutions is that (near) steady state values of
statistics of interest, such as the splitting and survival probabilities, can
be enumerated without the need to perform simulations over long time periods.
In Equations (19)-(24), we derived expressions for the time-varying left and
right splitting probabilities and survival probability. We verify that these
probabilities are correct for both the single and multiple growing domains by
comparing the derived expressions against the probabilities obtained from
repeated realisations of the random walk (Appendix A, Figure 13). We highlight
the usefulness of the expressions by exploring how domain growth rate, initial
location and domain-dependent diffusivity influence the splitting and survival
probabilities in Figure 7. Uniformly increasing the diffusivity increases both
the splitting probabilities while decreasing the survival probability.
Interestingly, comparing the results where $D_{1}=0.25$ and $D_{2}=0.5$
(Figures 7(d)-(f)) against the results where $D_{1}=0.5$ and $D_{2}=0.25$
(Figures 7(g)-(i)) suggests that, in terms of an individual reaching the
external boundary, an increased rate of movement in the outer domain is more
important than in the inner domain.

Figure 8: (a) Survival probability over time as defined in Equations (22)-(24)
for multiple oscillating domains (Equations (29)-(30)). The grey line
corresponds to no oscillation ($\beta_{i}=0$). The solid lines correspond to
$L_{i,\text{min}}=50$. The dashed lines correspond to $L_{i,\text{min}}=25$.
Oscillation rates are $\beta_{i}=-8\pi\times 10^{-4}$ (cyan), $-4\pi\times
10^{-4}$ (black), $4\pi\times 10^{-4}$ (orange), $8\pi\times 10^{-4}$ (pink).
(b),(c) Survival probability at (b) $t=8\times 10^{3}$ and (c) $t=1\times
10^{4}$ for a range of oscillation rates and minimum domain sizes. (d)
Survival probability at $t=8\times 10^{3}$ (black) and $t=1\times 10^{4}$
(pink) for $L_{i,\text{min}}=35$. The black line in (b) and the pink line in
(c) correspond to the profiles in (d). Parameters used are
$L_{1}(0)=L_{2}(0)=75$, $D_{1}=D_{2}=0.5$, and $x_{0}=0$.

Finally, we consider the question of survival probabilities on oscillatory
domains. That is, we consider domains that oscillate between a minimum (but
positive) and maximum size as defined in Equation (29). Contractile cells,
such as smooth muscle cells and cardiomyocytes, exhibit this behaviour [1].
Cardiomyocytes contract and relax with a regular rhythm as the heart beats,
while various molecular species undergo diffusion within the cell [1]. We
examine how the minimum domain size ($L_{i,\text{min}}$) and the rate of
oscillation ($\beta_{i}$) impact the survival probability, and present the
results in Figure 8. We see that in all cases, after one oscillation has been
completed, an individual is more likely to have reached the external boundary
in comparison to the case where there is no oscillation (Figure 8(a)). This is
interesting as the average length of the domain is not altered by the presence
of oscillation. However, as the size of individuals’ movement events is
independent of the overall size of the domain, a single movement event on a
domain that is at its minimum size is proportionally larger than if at the
domain is at its initial size. We see that for smaller minimum domain sizes
the survival probability decreases more rapidly. We explore how changing the
rates of oscillation and minimum domain sizes affects the survival probability
at $t=8\times 10^{3}$ and $t=1\times 10^{4}$ (Figures 8(b),(c)). The times and
oscillation rates are chosen such that at $t=1\times 10^{4}$ all domains have
completed an integer number of oscillations and are at the initial domain
size. Interestingly, in Figure 8(c), we see that all oscillation rates except
for $\beta_{i}=0$ have the same survival probability at $t=1\times 10^{4}$.
This suggests that the survival probability is insensitive to the rate of
oscillation; rather it is impacted by the current phase of oscillatory growth
and the magnitude of the oscillation. If we consider a fixed minimum domain
size and vary the rate of oscillation, as in Figure 8(d), we can see this
impact clearly. There is an increased sensitivity to the phase of growth if
there have been few oscillations, which occurs for $\beta_{i}$ values with a
low magnitude. This sensitivity decreases as the magnitude of the oscillation
rate increases. Again, when the oscillating domains align in size (at
$t=1\times 10^{4}$) we see that the rate of oscillation does not affect the
survival probability, provided there is any oscillation. In the context of
contracting cardiomyocytes, this suggests that the traversal of chemical
species across the cell may not be impacted by transient changes to heart
rate.

##  4 Discussion and conclusions

Diffusive processes occur on expanding and contracting domains, from the
smallest of biological scales [1, 8] to the largest of cosmological scales
[3]. Mathematical investigations have yielded detailed insight into the
dynamics of diffusive processes on both single uniformly growing domains and
multiple non-growing domains [21, 22, 23, 10, 11, 12, 13, 14, 15, 16, 17, 18,
19, 20]. However, certain diffusive processes of interest, such as drug
delivery and cell migration, occur on domains that can grow in a spatially
non-uniform manner [5, 8].

Here we have presented exact solutions to a mathematical model of a diffusive
process on multiple growing domains. Each domain exhibits spatially-uniform
growth and hence the overall domain can exhibit spatially non-uniform growth.
We compare the enumerated exact solution profiles against profiles obtained
from repeated realisations of a corresponding lattice-based random walk to
ensure that the exact solutions are valid. In all cases, we observe that the
derived exact solutions accurately describe the average random walk behaviour.
From our exact solutions, we derive expressions for relevant time-varying
statistics such as the survival probability and splitting probabilities. We
reveal how domain-specific model parameters influence these statistics,
including how the interplay between domain growth rates, domain-specific
diffusivities and initial location drive long-term survival. Finally, we show
how oscillating domain growth enhances diffusion. Intriguingly, the rate of
oscillation (provided it is non-zero) does not influence the survival
probability, provided the comparison is made at a point of time such that the
domains have completed an integer number of oscillations. Instead, the
magnitude of the oscillation is the key factor influencing the survival
probability. This result has interesting implications for biological phenomena
that exhibit oscillatory growth. For example, in the context of
cardiomyocytes, this result suggests that the internal diffusion of chemical
species should be robust to changes in heart rate, such as those due to
exercise or stress.

The work presented here may be extended in several directions. We have only
considered linearly diffusive processes (i.e. random motion of the agents).
However, many processes exhibit characteristics attributable to nonlinear
diffusion, such as compact support and a finite rate of spread [31, 32, 33].
These characteristics are particularly relevant for growing domains, as the
relationship between the rate of spread and the rate of domain growth dictates
whether it is possible for agents to reach the external domain boundary. It is
instructive to determine which nonlinear diffusivity functions permit exact
solutions following the transformation techniques presented here. Similarly,
many biological processes require the inclusion of a reaction term to
accurately capture the underlying population dynamics [31]. It is unclear
which classes of reaction terms would admit exact solutions under the
transformations examined here; however, it is possible that non-classical
transformation techniques may yield analytical progress [34].

## Code availability

The code used to generate the results presented here can be found on Github
at:
https://github.com/DrStuartJohnston/heterogeneous-growing-domains.

## Acknowledgements

The authors would like to acknowledge the organisers of the MATRIX workshop
“The mathematics of tissue dynamics,” where this research project commenced.

## Funding

This work was in part supported by the Australian Research Council (STJ:
DE200100988, MJS: DP200100177).

## References

  * [1] D. M. Bers, Cardiac excitation–contraction coupling, Nature 415 (6868) (2002) 198–205.
  * [2] K. A. Landman, M. J. Simpson, D. F. Newgreen, Mathematical and experimental insights into the development of the enteric nervous system and hirschsprung’s disease, Development, Growth & Differentiation 49 (4) (2007) 277–286.
  * [3] A. G. Riess, A. V. Filippenko, P. Challis, A. Clocchiatti, A. Diercks, P. M. Garnavich, R. L. Gilliland, C. J. Hogan, S. Jha, R. P. Kirshner, et al., Observational evidence from supernovae for an accelerating universe and a cosmological constant, The Astronomical Journal 116 (3) (1998) 1009.
  * [4] B. I. Shraiman, Mechanical feedback as a possible regulator of tissue growth, Proceedings of the National Academy of Sciences 102 (9) (2005) 3318–3323.
  * [5] B. J. Binder, K. A. Landman, M. J. Simpson, M. Mariani, D. F. Newgreen, Modeling proliferative tissue growth: a general approach and an avian case study, Physical Review E 78 (3) (2008) 031912.
  * [6] E. J. Crampin, W. W. Hackborn, P. K. Maini, Pattern formation in reaction-diffusion models with nonuniform domain growth, Bulletin of Mathematical Biology 64 (2002) 747–769.
  * [7] C. A. Yates, R. E. Baker, R. Erban, P. K. Maini, Going from microscopic to macroscopic on nonuniform growing domains, Physical Review E 86 (2) (2012) 021921\.
  * [8] O. Zavgorodnya, C. A. Carmona-Moran, V. Kozlovskaya, F. Liu, T. M. Wick, E. Kharlampieva, Temperature-responsive nanogel multilayers of poly (n-vinylcaprolactam) for topical drug delivery, Journal of Colloid and Interface Science 506 (2017) 589–602.
  * [9] J. Crank, The Mathematics of Diffusion, Oxford University Press, 1979.
  * [10] E. J. Crampin, E. A. Gaffney, P. K. Maini, Reaction and diffusion on growing domains: scenarios for robust pattern formation, Bulletin of Mathematical Biology 61 (6) (1999) 1093–1120.
  * [11] E. J. Crampin, E. A. Gaffney, P. K. Maini, Mode-doubling and tripling in reaction-diffusion patterns on growing domains: A piecewise linear model, Journal of Mathematical Biology 44 (2) (2002) 107–128.
  * [12] K. A. Landman, G. J. Pettet, D. F. Newgreen, Mathematical models of cell colonization of uniformly growing domains, Bulletin of Mathematical Biology 65 (2) (2003) 235–262.
  * [13] F. Le Vot, S. Yuste, Continuous-time random walks and fokker-planck equation in expanding media, Physical Review E 98 (4) (2018) 042117.
  * [14] V. Klika, E. A. Gaffney, History dependence and the continuum approximation breakdown: the impact of domain growth on turing’s instability, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 473 (2199) (2017) 20160744.
  * [15] D. R. McLean, B. P. Graham, Mathematical formulation and analysis of a continuum model for tubulin-driven neurite elongation, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 460 (2048) (2004) 2437–2456.
  * [16] R. J. Ross, R. E. Baker, C. A. Yates, How domain growth is implemented determines the long-term behavior of a cell population through its effect on spatial correlations, Physical Review E 94 (1) (2016) 012408.
  * [17] A. Ryabov, E. Berestneva, V. Holubec, Brownian motion in time-dependent logarithmic potential: exact results for dynamics and first-passage properties, The Journal of Chemical Physics 143 (11) (2015) 114117.
  * [18] M. J. Simpson, J. A. Sharp, L. C. Morrow, R. E. Baker, Exact solutions of coupled multispecies linear reaction–diffusion equations on a uniformly growing domain, PloS One 10 (9) (2015) e0138894.
  * [19] M. J. Simpson, J. A. Sharp, R. E. Baker, Survival probability for a diffusive process on a growing domain, Physical Review E 91 (4) (2015) 042701.
  * [20] S. B. Yuste, E. Abad, C. Escudero, Diffusion in an expanding medium: Fokker-planck equation, green’s function, and first-passage properties, Physical Review E 94 (3) (2016) 032118.
  * [21] E. Carr, I. Turner, A semi-analytical solution for multilayer diffusion in a composite medium consisting of a large number of layers, Applied Mathematical Modelling 40 (15-16) (2016) 7034–7050.
  * [22] E. J. Carr, G. Pontrelli, Modelling mass diffusion for a multi-layer sphere immersed in a semi-infinite medium: application to drug delivery, Mathematical Biosciences 303 (2018) 1–9.
  * [23] E. J. Carr, M. J. Simpson, New homogenization approaches for stochastic transport through heterogeneous media, The Journal of Chemical Physics 150 (4) (2019) 044104.
  * [24] M. J. Simpson, R. E. Baker, Exact calculations of survival probability for diffusion on growing lines, disks, and spheres: the role of dimension, The Journal of Chemical Physics 143 (9) (2015) 094109.
  * [25] I. Barrass, E. J. Crampin, P. K. Maini, Mode transitions in a model reaction–diffusion system driven by domain growth and noise, Bulletin of Mathematical Biology 68 (2006) 981–995.
  * [26] B. J. Binder, K. A. Landman, Exclusion processes on a growing domain, Journal of Theoretical Biology 259 (3) (2009) 541–551.
  * [27] C. A. Yates, Discrete and continuous models for tissue growth and shrinkage, Journal of Theoretical Biology 350 (2014) 37–48.
  * [28] M. N. Nilsson, The moving-eigenvalue method: hitting time for ito processes and moving boundaries, Journal of Physics A: Mathematical and Theoretical 53 (40) (2020) 405201.
  * [29] R. J. Ross, C. A. Yates, R. E. Baker, The effect of domain growth on spatial correlations, Physica A: Statistical Mechanics and its Applications 466 (2017) 334–345.
  * [30] K. Gudnason, S. Sigurdsson, B. S. Snorradottir, M. Masson, F. Jonsdottir, A numerical framework for drug transport in a multi-layer system with discontinuous interlayer condition, Mathematical Biosciences 295 (2018) 11–23.
  * [31] S. T. Johnston, R. E. Baker, D. McElwain, M. J. Simpson, Co-operation, competition and crowding: a discrete framework linking allee kinetics, nonlinear diffusion, shocks and sharp-fronted travelling waves, Scientific Reports 7 (1) (2017) 1–19.
  * [32] Y. Li, P. R. Buenzli, M. J. Simpson, Interpreting how nonlinear diffusion affects the fate of bistable populations using a discrete modelling framework, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 478 (2262) (2022) 20220013.
  * [33] J. L. Vázquez, The porous medium equation: mathematical theory, Oxford University Press, 2007.
  * [34] B. H. Bradshaw-Hajek, Nonclassical symmetry solutions for non-autonomous reaction-diffusion equations, Symmetry 11 (2) (2019) 208.

## Appendix A. Additional solutions

Figure 9: Comparison between the average behaviour in the lattice-based random
walk $\overline{C}(x,t)$ (black, dashed) and the exact solution $C(x,t)$ as
defined in Equation (10) (cyan) for a single fixed domain. Parameters used are
$D=0.5$, $L=100$, $N=1000$. Solution profiles are presented at $t=100$,
$t=250$ and $t=500$ The arrow indicates the direction of increasing time.
Average random walk behaviour is obtained from 5000 identically-prepared
realisations of the random walk. Figure 10: Comparison between the average
behaviour in the lattice-based random walk $\overline{C}(x,t)$ (black, dashed)
and the exact solution $C(x,t)$ as defined in Equation (10) (cyan) for a
single linearly growing domain (Equations (25)-(26)). Parameters used are
$D=0.5$, $L=50$, $N=1000$, $\beta=0.05$. Solution profiles are presented at
$t=200$, $t=500$ and $t=1000$. The arrow indicates the direction of increasing
time. Average random walk behaviour is obtained from 5000 identically-prepared
realisations of the random walk. Figure 11: Comparison between the average
behaviour in the lattice-based random walk $\overline{C}(x,t)$ (black, dashed)
and the exact solutions $C_{1}(x,t)$ (cyan), $C_{2}(x,t)$ (orange) and
$C_{3}(x,t)$ (pink) as defined in Equations (15)-(18) for multiple fixed
domains with domain-independent diffusivities. Parameters used are
$D_{1}=D_{2}=0.5$, $L_{1}=L_{2}=50$, $N=1000$, (a)$x_{0}=0$, (b) $x_{0}=20$.
Solution profiles are presented at $t=200$, $t=500$ and $t=1000$. The arrow
indicates the direction of increasing time. Average random walk behaviour is
obtained from 5000 identically-prepared realisations of the random walk.
Figure 12: Comparison between the average behaviour in the lattice-based
random walk $\overline{C}(x,t)$ (black, dashed) and the exact solutions
$C_{1}(x,t)$ (cyan), $C_{2}(x,t)$ (orange) and $C_{3}(x,t)$ (pink) as defined
in Equations (15)-(18) for multiple fixed domains and domain-independent
diffusivities. Parameters used are $L_{1}=L_{2}=50$, $N=1000$, $x_{0}=0$, (a)
$D_{1}=0.5$, $D_{2}=0.25$, (b) $D_{1}=0.25$, $D_{2}=0.5$, (c) $D_{1}=0.5$,
$D_{2}=0.05$. Solution profiles are presented at $t=200$, $t=500$ and
$t=1000$. The arrow indicates the direction of increasing time. Average random
walk behaviour is obtained from 5000 identically-prepared realisations of the
random walk. Figure 13: Comparison between the survival (black, cyan), left
splitting (dark grey, orange) and right splitting (light grey, pink)
probabilities obtained from the lattice-based random walk (solid) and the
exact solutions defined in (a) Equations (19)-(21), (b) Equations (22)-(24)
(dashed) for (a) one linearly growing domain and (b) two exponentially growing
domains (Equations (25)-(26)). Parameters used are $N=1000$, (a)
$L_{1}(0)=50$, $D_{1}=0.5$ and $\beta_{1}=0.05$ and (b)
$L_{1}(0)=L_{2}(0)=50$, $D_{1}=D_{2}=0.5$ and $\beta_{1}=\beta_{2}=0.05$.
Probabilities are presented at $t=10^{4}$. Average random walk behaviour is
obtained from 100 identically-prepared realisations of the random walk.