4pt 4pt

11institutetext: A. Sargood 22institutetext: E. A. Gaffney 33institutetext: A.
L. Krause 44institutetext: Wolfson Centre for Mathematical Biology,
Mathematical Institute, University of Oxford, Andrew Wiles Building, Radcliffe
Observatory Quarter, Woodstock Road, Oxford, OX2 6GG, United Kingdom
A. L. Krause 55institutetext: Mathematical Sciences Department, Durham
University, Upper Mountjoy Campus, Stockton Rd, Durham DH1 3LE, United Kingdom

# Fixed and Distributed Gene Expression Time Delays in Reaction-Diffusion
Systems

Alec Sargood  Eamonn A. Gaffney  Andrew L. Krause

(Received: date / Accepted: date)

###### Abstract

Time delays, modelling the process of intracellular gene expression, have been
shown to have important impacts on the dynamics of pattern formation in
reaction-diffusion systems. In particular, past work has shown that such time
delays can shrink the Turing space, thereby inhibiting patterns from forming
across large ranges of parameters. Such delays can also increase the time
taken for pattern formation even when Turing instabilities occur. Here we
consider reaction-diffusion models incorporating fixed or distributed time
delays, modelling the underlying stochastic nature of gene expression
dynamics, and analyze these through a systematic linear instability analysis
and numerical simulations for several sets of different reaction kinetics. We
find that even complicated distribution kernels (skewed Gaussian probability
density functions) have little impact on the reaction-diffusion dynamics
compared to fixed delays with the same mean delay. We show that the location
of the delay terms in the model can lead to changes in the size of the Turing
space (increasing or decreasing) as the mean time delay, $\tau$, is increased.
We show that the time to pattern formation from a perturbation of the
homogeneous steady state scales linearly with $\tau$, and conjecture that this
is a general impact of time delay on reaction-diffusion dynamics, independent
of the form of the kinetics or location of the delayed terms. Finally we show
that while initial and boundary conditions can influence these dynamics,
particularly the time-to-pattern, the effects of delay appear robust under
variations of initial and boundary data. Overall our results help clarify the
role of gene expression time delays in reaction-diffusion patterning, and
suggest clear directions for further work in studying more realistic models of
pattern formation.

###### Keywords:

Pattern formation time delay linear instability analysis

##  1 Introduction

The self-organisation of cells into an apparent order appears across many
different fields within biology. For example, the distribution of cells during
the developmental process of an embryo, the growth of cancerous tissue morph
(2), vertebrate limb development miura1 (30, 18, 29), and pattern formation on
animal coats (e.g. spots on a jaguar painter (33), feathers on birds bailleul
(1)). Wolpert wolpert (49) presented the idea that, underpinning the
development of shape and form (morphogenesis) is a cell’s ability to
differentiate according to its position in space and time. Furthermore, the
concentration of signalling molecules (morphogens), or the concentration
gradients of certain morphogens across a spatial domain of cells, affects the
cell differentiation mechanism, and thus cells adopt a state relative to the
concentration of a specific morphogen that they are exposed to.

The mechanism allowing cells to adopt an appropriate state is known as
differential gene expression, and depends crucially on the communication
between cells, achieved through cellular signalling gaffmonk (14). Typical
reaction-diffusion systems modelling such signalling implicitly assume a
negligible timescale on which the internal cellular signalling and gene
expression processes occur. The gene expression process, however, is extremely
complex and proceeds through several stages gaffmonk (14), including gene
transcription and gene translation. These sub-processes can take large amounts
of time, and it has been experimentally shown that these time delays are
typically on the order of minutes, but in some cases can be as large as
several hours gaffmonk (14, 47). Such time delays can therefore approach the
timescale of pattern formation process itself. For example, the basic body
plan of a zebrafish is established in less than 24 hours kimmel (21). Hence it
seems important to assess the impact of such time delays on pattern formation
in developmental biology.

In 1952, Alan Turing proposed that the pattern formation process could be
mathematically modelled on a purely chemical basis via the interaction of
morphogens, whose evolutions are described by a system of coupled reaction-
diffusion equations turing (48). Turing showed that a stable steady state,
robust to small perturbations in the spatially homogeneous setting (i.e.
without diffusion), could become unstable and sensitive to small perturbations
with the introduction of diffusion, leading to spatially inhomogeneous
patterns. Cell fate decisions are then based on these morphogen
concentrations, where regions of high morphogen or low morphogen concentration
can lead to different cell fate decisions. Turing’s model is therefore one of
pre-patterning, where the morphogen pattern concentrations across a spatial
domain are modelled, which in turn lead to cell fate decisions at a later
stage. Typical reaction-diffusion systems in the context of Turing pattern
formation consist of two partial differential equations describing the
interaction and evolution of two morphogens, the activator and inhibitor.
Empirical evidence suggests that Turing instabilities are present in real
biological systems, and can be used to explain complex biological phenomena
yigaffneyli (52, 17, 31, 29, 46). However, whether Turing patterns may be
found experimentally in biological systems with simple two-species systems is
still very much an active field of research bespoke (51).

Time delays have been investigated in the context of Turing patterns, both
numerically and analytically, through the incorporation of constant fixed time
delays. One of the canonical reaction-diffusion mechanisms that exhibits
Turing instabilities is the Schnakenberg model schnakenberg (44), also known
as the activator-depleted model. This system has been extensively studied in
the context of Turing pattern formation with incorporated gene expression time
delays. Two biologically motivated variants incorporating time delays are the
ligand-internalisation (LI) and reverse ligand-binding (RLB) models leegaffney
(25). The LI model places gene expression time delay in purely the activator
dynamics, whereas the RLB model contains time delay in both the activator and
inhibitor dynamics.

Numerical results in gaffmonk (14) on the LI model showed that the time taken
until pattern formation occurs drastically increases as the time delay in the
model is increased. In particular, small delays on the order of minutes can
cause a large increase in time-to-pattern, on the order of several hours,
compared to a model with no time delay. This highlights the impact of
including gene expression time delays, especially when considering patterning
events that occur on a fast timescale.

The two papers jiang (20, 52) consider both the LI and RLB variants of the
Schnakenberg model. Using linear analysis and numerical simulations, the
results in both papers suggest that the RLB model can exhibit spatially
inhomogeneous temporal oscillations, as well as de-stabilisation of spatially
inhomogeneous steady states, inhibiting pattern formation via Turing
instabilities. The results in yigaffneyli (52) indicate that extensive ligand
internalization may inhibit the formation of patterns within some regions of
the parameter space. Results in both leegaffney (25, 26) suggest that time
delay causes a significant effect on the time taken until pattern formation
occurs. A final observation from gaffmonk (14, 26), for both LI and RLB
models, is that an increasing time delay may also increase the sensitivity of
the final pattern formed to variations in initial conditions.

In contrast, analysis of one-dimensional spike solutions of the Gierer-
Meinhardt (GM) model in fadai1 (12, 11) shows that the placement of delay
terms in the model can affect the size of the parameter regimes for which the
spike solution is linearly stable. It was found that, depending on the
positioning of time delayed terms, an increasing time delay can have a
stabilising or de-stabilising effect, enlarging or shrinking the stable
parameter region of the spike. Further details of spike solutions of the GM
model and their stability analysis can also be found in spike (19). This
analysis highlighted the importance of time delay positioning in the GM model
for the stability of spike solutions.

The key insights from the literature on delayed reaction-diffusion models
above are that fixed time delays can (a) change the space of parameters in
which spatially inhomogeneous steady states are observed (the Turing space),
and (b) delay the time for patterns to form. These results have implications
on the use of such models as mechanistic explanations of pattern formation in
development. Namely, in cases where the Turing space shrinks for realistic
sizes of delay, this can challenge the robustness of Turing pattern formation,
which is already an open problem in using non-delayed reaction-diffusion
models to explain robust biological phenomena maini2012turing (27, 45).
Additionally, the increasing time needed for patterns to form can become an
obstacle in quantitatively relating these models to developmental perspectives
for patterning events that occur on relatively fast timescales leegaffmonk
(26). Thus our first objective will be to explore these observations in
variants of Schnakenberg and GM systems in more detail by explicitly
constructing Turing spaces via a linear stability analysis, and exploring the
time taken for patterns to form as a function of the delay time.

The current literature on Turing pattern formation in development has
primarily considered gene expression dynamics via fixed time delays in
reaction-diffusion processes. On a cellular level however, the biological
processes responsible for gene expression are inherently stochastic raj (37,
10, 28, 35). Distributed time delays have been incorporated to model
biological phenomena such as hematopoiesis, and lactose operon dynamics
newdist (4), Wnt/$\beta$-catenin signalling pathways signal (5), and Oncolytic
virotherapy treatments for cancer cancer (9), among other applications. A
distributed delay can be thought of as a more ‘general and realistic’ cancer
(9) approach to modelling, on a larger scale, a process which may possess an
intrinsic stochasticity on a small scale. Within the context of Turing pattern
formation, introducing a fixed time delay into the reaction-diffusion
mechanism is a simplification of the underlying biological process on a
microscopic level. This leads us to our second objective: to model distributed
gene expression time delays at the macroscopic level to ascertain if the
results obtained in the case of fixed time delays are sensitive to the shape
of a distribution of delays, which is in principle a more realistic model for
gene expression dynamics bratsun (3, 22). Of particular interest within this
objective is to assess whether the distribution of time delays induces
substantial changes compared to a Dirac-delta kernel for the time delay,
noting the latter is equivalent to a fixed time delay.

In Section 2 we outline a general formulation of reaction-diffusion systems
with distributed time delay, as well as some specific reaction kinetics that
we will examine. We briefly review the underlying theory of Turing pattern
formation and carry out a linear stability analysis of our general distributed
delay model in Section 3. In Section 4, we apply this instability analysis to
explore how Turing spaces change for increasing time delays in fixed and
distributed scenarios. We compare these predictions from the linear theory
against numerical solutions in Section 5, as well as exploring questions of
the time-to-pattern as a function of delay, boundary, and initial conditions.
We close with a discussion of our results in Section 6, further highlighting
the importance of the placement of time delayed terms within the models, and
noting the surprising fact that distributed time delays have almost no
influence on pattern formation dynamics compared to fixed delay models with
the same mean delay.

##  2 Models of Gene Expression Time Delay

We consider two morphogens, $u$ and $v$, obeying the following non-dimensional
reaction-diffusion equations:

$\begin{split}\frac{\partial u}{\partial
t}&=\frac{\epsilon^{2}}{L^{2}}\frac{\partial^{2}u}{\partial
x^{2}}+f(u,v,\mathcal{F}(u,v)),\\\ \frac{\partial v}{\partial
t}&=\frac{1}{L^{2}}\frac{\partial^{2}v}{\partial
x^{2}}+g(u,v,\mathcal{G}(u,v)),\end{split}$ (1)

where $x\in\Omega=[0,1]$, $t>0$ and

$\mathcal{F}(u,v)=\int_{\tau_{1}}^{\tau_{2}}K(s;\bm{p})F(u_{s},v_{s})\
\text{d}s,\quad\mathcal{G}(u,v)=\int_{\tau_{1}}^{\tau_{2}}K(s;\bm{p})G(u_{s},v_{s})\
\text{d}s,$

with $u_{s}=u(x,t-s)$ and $v_{s}=v(x,t-s)$ the function evaluations of $u$ and
$v$ at some time delay $s$, which is an integration variable. We also have
that $\tau_{1}$ and $\tau_{2}$ are the minimum and maximum time delays of the
distributions, $L^{2}$ is the non-dimensional scaling of the domain length,
and $\epsilon^{2}$ the diffusion ratio between the activator $u$ and inhibitor
$v$. Unless otherwise stated, we use the same value of $\epsilon^{2}=0.001$ as
in gaffmonk (14), and a non-dimensional domain size, $L^{2}=4.5$. The latter
allows larger numbers of stripes to form should the system pattern compared to
the domain length scales used in gaffmonk (14), enabling the prospect of
greater sensitivity of the final pattern when comparing different
distributions of time delay, since pattern mode selection is particularly
sensitive on larger domains crampin (6). We primarily use no flux (homogeneous
Neumann) boundary conditions on the boundary of the spatial domain, namely

$\frac{\partial u}{\partial x}=\frac{\partial v}{\partial x}=0,\quad\quad
x=0,1,$ (2)

though we will consider an alternative set of boundary conditions in Section
5.3.

The functions $f$ and $g$ are the reaction kinetics, which we will vary in
subsequent sections. We will consider kinetics that have a unique positive
homogeneous steady state, $(u_{\star},v_{\star})$. Unless otherwise stated,
initial conditions $(u_{0},v_{0})$ are chosen as a small random Gaussian
perturbation from this homogeneous steady state given by

$\begin{pmatrix}u_{0}\\\
v_{0}\end{pmatrix}=\begin{pmatrix}u_{\star}(1+r_{u}(x))\\\
v_{\star}(1+r_{v}(x))\end{pmatrix},$ (3)

where $r_{u}(x)$, $r_{v}(x)$ are Normally distributed random variables for
each $x\in[0,1]$ with zero mean and standard deviation $\sigma_{IC}$. Unless
otherwise mentioned, we take $\sigma_{IC}=0.01$ throughout.

The stochastic nature of gene expression delays leads us to consider a mean-
field approach to modelling the time delay bratsun (3, 22). The function $K$
is the kernel of our distributed delay, representing the distribution of delay
times from the underlying process. The functions $F$ and $G$ model reaction
steps which incorporate this delay. By assuming each individual mechanism
within the gene expression process occurs independently and identically, we
use the central limit theorem to model the delay as a (truncated) Gaussian
distribution, $K(s;\bm{p})$, with $s$ the integral variable
$s\in[\tau_{1},\tau_{2}]$, $0<\tau_{1}<\tau_{2}$ and $\bm{p}$ the vector of
distribution parameters. We will also consider non-Normal distributions by
investigating skewed Gaussian distributions as well.

###  2.1 Symmetric Truncated Gaussian Formulation

Using a symmetric truncated Gaussian kernel, we have $\bm{p}=(\tau,\sigma)$,
for some mean delay $\tau$ and standard deviation $\sigma$, with integration
limits chosen as $\tau_{1}=\tau-n\sigma$ and $\tau_{2}=\tau+n\sigma$ for some
$n\in\mathbb{N}$, such that $\tau_{1}=\tau-n\sigma>0$ to ensure positive time
delays only. The distribution is then given by

$K(s;\bm{p})=K_{N}(s;\tau,\sigma)=\Phi_{c}\frac{1}{\sigma\sqrt{2\pi}}\exp\left(-\frac{1}{2}\left(\frac{s-\tau}{\sigma}\right)^{2}\right),$
(4)

where $\Phi_{c}$ denotes the truncation scaling constant, which ensures that
$K_{N}(s;\tau,\sigma)$ integrates to $1$ over the given integration domain
$[\tau_{1},\tau_{2}]$, and is computed as

$\Phi_{c}=\frac{1}{\phi\left(\frac{\tau_{2}-\tau}{\sigma}\right)-\phi\left(\frac{\tau_{1}-\tau}{\sigma}\right)},$
(5)

with $\phi(x)$ the cdf of the standard Gaussian distribution,111The error
function is given by
$\text{erf}(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-z^{2}}\ dz$.

$\phi(x)=\frac{1}{2}\left(1+\text{erf}\left(\frac{x}{\sqrt{2}}\right)\right).$
(6)

Throughout this paper, we use $n=3$ so that the integration limits are
$\tau_{1}=\tau-3\sigma$ and $\tau_{2}=\tau+3\sigma$. This was chosen so that a
relatively large $\sigma$ value could be used for each $\tau$ while
maintaining $\tau_{1}>0$. For each $\tau$ value, a maximum $\sigma$ value can
be computed such that $\tau_{1}=\tau-3\sigma\geq 0$ as
$\sigma_{\max}=\frac{\tau}{3}$. By setting $\sigma<\sigma_{\max}$, we ensure
that the integration domain strictly considers positive time delays only.

As $\sigma\to 0$, we have that $K_{N}(s;\tau,\sigma)\to\delta(s-\tau)$. Using
the sifting property of the delta function, we remark that the continuous time
delay problem formulation in (1) reduces in the limit of $\sigma\to 0$ to a
discrete fixed delay case with delay $\tau$, where

$\mathcal{F}(u,v)=F(u_{\tau},v_{\tau}),\quad\quad\mathcal{G}(u,v)=G(u_{\tau},v_{\tau}).$

Hence, the continuously distributed time delay model is a generalisation of
the fixed time delay model.

###  2.2 Skewed Gaussian Formulation

We will also study the role that the shape of a distribution has by
considering a skewed Gaussian kernel, using the parameters
$\bm{p}=(\mu,\omega,\rho)$. Such a distribution exhibits asymmetry but is
otherwise a relatively simple generalization of the Normal distribution. The
probability density function of the skewed truncated Gaussian distribution is
given by

$K(s;\bm{p})=K_{S}(s;\mu,\omega,\rho)=\frac{\Psi_{c}}{\omega}\sqrt{\frac{2}{\pi}}\exp\left(-\frac{1}{2}\left(\frac{s-\mu}{\omega}\right)^{2}\right)\phi\left(\rho\frac{s-\mu}{\omega}\right),$
(7)

where $\phi(x)$ is the same as in (6). The parameter $\rho$ denotes the skew
factor. The distribution is negatively skewed for $\rho<0$ and positively
skewed for $\rho>0$. Finally, we have that $\Psi_{c}$ is the truncation
scaling constant, given as

$\Psi_{c}=\frac{1}{\Psi\left(\frac{\tau_{2}-\mu}{\omega},\rho\right)-\Psi\left(\frac{\tau_{1}-\mu}{\omega},\rho\right)},$
(8)

with $\Psi(x,\rho)$ the cdf of a skewed Gaussian distribution, described by

$\Psi(x,\rho)=\phi(x)-2T(x,\rho).$ (9)

The function $T(x,\rho)$ denotes the Owen’s T function owenst (34) and is
written as an integral in the form

$T(x,\rho)=\frac{1}{2\pi}\int_{0}^{\rho}\frac{e^{-\frac{1}{2}x^{2}(1+s^{2})}}{1+s^{2}}\
\text{d}s\quad-\infty<x,\rho<\infty.$ (10)

In the computational implementation of the skewed truncated Gaussian pdf, the
integral $T(x,\rho)$ is resolved numerically using the composite Simpson’s
rule with $100,000$ discretisation points.

The parameters $\mu$ and $\omega$ no longer denote the mean and standard
deviation of the distribution, but instead the location of the maximum and a
scaling factor. To compare how the skewed distribution affects the onset of
patterning compared to that of the fixed delay case, we consider the mean of
the skewed distribution, $\tau$, which is given by

$\tau=\int_{\tau_{1}}^{\tau_{2}}sK_{S}(s;\mu,\omega,\rho)\ \ \text{d}s.$ (11)

Following skewed (13), the mean of the skewed truncated Gaussian distribution
is computed as

$\tau=\mu+\omega\Psi_{c}\left[K_{S}(\tau_{1};\mu,\omega,\rho)-K_{S}(\tau_{2};\mu,\omega,\rho)+\frac{2\rho}{\hat{\rho}\sqrt{2\pi}}\left(\phi\left(\hat{\rho}\frac{\tau_{2}-\mu}{\omega}\right)-\phi\left(\hat{\rho}\frac{\tau_{1}-\mu}{\omega}\right)\right)\right],$
(12)

with $\hat{\rho}=\left(1+\rho^{2}\right)^{1/2}$. We note that, for a given
$\rho$, all of the terms on the right-hand side of (12), namely $\omega$,
$\Psi_{c}$, and $K_{S}(s;\mu,\omega,\rho)$, can be written explicitly in terms
of $\mu$. Equation (12) can therefore be solved implicitly for $\mu(\tau)$,
for a given $\tau$.

The integration limits were set to $\tau_{1}=\mu-3\omega$,
$\tau_{2}=\mu+3\omega$, where $\omega$ was chosen such that
$\omega<\omega_{\max}$, with $\omega_{\max}=\frac{\mu}{3}$ to ensure only
positive time delays were considered. In Figure 1 we give examples of these
Skew Normal distributions compared to a standard symmetric Normal
distribution, for different values of the standard deviation/scaling in each
panel.

(a) $\sigma=0.05\sigma_{\textrm{max}}$, $\omega=0.05\omega_{\textrm{max}}$

(b) $\sigma=0.2\sigma_{\textrm{max}}$, $\omega=0.2\omega_{\textrm{max}}$

(c) $\sigma=0.99\sigma_{\textrm{max}}$, $\omega=0.99\omega_{\textrm{max}}$

Figure 1: Probability density functions of the (symmetric) Normal, and Skew
Normal distributions with a fixed mean of $\tau=1$ in all cases, and different
standard deviations. The positively skewed distributions have $\rho=10$, and
the negatively skewed $\rho=-10$. Note that the axes labels cover very
different ranges, with small values of $\sigma$ and $\omega$ leading to very
narrow distributions. For large values of $\omega$, the negative and positive
skew cases are not symmetric about the mean due to the truncation used. For
this value of $\tau$, we have that $\sigma_{\textrm{max}}\approx 0.3333$ and
$\omega_{\textrm{max}}\approx 0.3354$.

###  2.3 Reaction Kinetics

We consider the Schnakenberg schnakenberg (44) and Gierer-Meinhardt (GM) gm
(15) models as two canonical reaction-diffusion systems that exhibit Turing
pattern formation. These models have also been studied extensively in the
context of fixed time delays as described above.

The stoichiometry of the Schnakenberg kinetics baker (50) is given by

$A\xrightleftharpoons[c_{-1}]{c_{1}}U,\quad B\xrightarrow{c_{2}}V,\quad
2U+V\xrightarrow{c_{3}}3U,$ (13)

where the $c_{i}$ represent reaction rates. The quantities $A$ and $B$ are
reservoirs whose evolution is not considered, and we assume a constant supply
rate. We use $u$, $v$, $a$, and $b$ to denote the (nondimensional)
concentrations of substances $U$, $V$, $A$ and $B$ respectively, with $a$ and
$b$ taken as constants. The form of the model we consider with time delay is
the Ligand Internalistion (LI) variant of the standard Schnakenberg model. The
LI model assumes that, a reaction at the cell surface is followed by
internalisation of a morphogen, before the gene expression process can
continue and morphogen production can occur leegaffney (25, 52). The time
delayed terms in the LI model only appear in the activator’s dynamics. This is
based on the assumption that the gene expression process, and thus the source
of the time delay, is responsible for autocatalysis of the activator in the
reaction-diffusion mechanism gaffmonk (14). Considering mass-action kinetics
and adding standard diffusion terms murray (32), we nondimensionalize the
resulting system to find that the LI model with a distributed time delay is
given by

$\begin{split}\frac{\partial u}{\partial
t}&=\frac{\epsilon^{2}}{L^{2}}\frac{\partial^{2}u}{\partial
x^{2}}+a-u-2u^{2}v+3\int_{\tau_{1}}^{\tau_{2}}K(s;\textbf{p})u_{s}^{2}v_{s}\
\text{d}s,\\\ \frac{\partial v}{\partial
t}&=\frac{1}{L^{2}}\frac{\partial^{2}v}{\partial x^{2}}+b-u^{2}v.\end{split}$
(14)

The steady state of the Schnakenberg model is given by
$(u_{\star},v_{\star})=\left(a+b,\frac{b}{(a+b)^{2}}\right)$. We do not
consider the RLB variant proposed in leegaffney (25), as it was shown to be
able to exhibit negative concentrations for positive values of the feed rates
$a$ and $b$ william (7). We will also consider how this model compares to a
fixed time delay LI model which is obtained in the limit $\sigma\to 0$.

Additionally, we consider delayed forms of the GM model in the fixed-delay
case in order to explore the role of different kinetics and delay terms on our
results. We will focus on two well-studied fixed time delay variants for
simplicity, leaving an analysis of distributed delay in GM models to future
work. A chemical interpretation of the kinetic reactions for the GM Model can
be found in leegaffmonk (26). The two non-dimensionalised model descriptions
we consider, with kinetic reactions taken from murray (32), and time delayed
terms motivated by fadai1 (12) and fadai2 (11), are given by (15) and (16),
and labelled GM1 and GM2 respectively, 222We note that the papers fadai1 (12,
11) label $v$ as the activator and $u$ as the inhibitor.

$\left.\begin{split}\frac{\partial u}{\partial
t}&=\frac{\epsilon^{2}}{L^{2}}\frac{\partial^{2}u}{\partial
x^{2}}+a-bu+\frac{u_{\tau}^{2}}{v_{\tau}},\\\ \frac{\partial v}{\partial
t}&=\frac{1}{L^{2}}\frac{\partial^{2}v}{\partial
x^{2}}+u_{\tau}^{2}-v,\end{split}\right\\}\text{GM}_{1}$ (15)
$\left.\begin{split}\frac{\partial u}{\partial
t}&=\frac{\epsilon^{2}}{L^{2}}\frac{\partial^{2}u}{\partial
x^{2}}+a-bu+\frac{u_{\tau}^{2}}{v},\\\ \frac{\partial v}{\partial
t}&=\frac{1}{L^{2}}\frac{\partial^{2}v}{\partial
x^{2}}+u_{\tau}^{2}-v,\end{split}\right\\}\text{GM}_{2}$ (16)

with $u_{\tau}=u(x,t-\tau)$ and $v_{\tau}=v(x,t-\tau)$. The key difference
between these models is the inhibitor term in the activator’s kinetics being
delayed or not. The homogeneous steady state of the GM model is given by

$(u_{\star},v_{\star})=\left(\frac{a+1}{b},\left(\frac{a+1}{b}\right)^{2}\right).$

##  3 General Linear Instability Analysis

We now consider the linear stability of homogeneous equilibria of the system
(1). Denoting the steady state as $(u_{\star},v_{\star})$, we consider a small
perturbation, $u(x,t)=u_{\star}+\delta\xi(x,t)$,
$v(x,t)=v_{\star}+\delta\eta(x,t)$, where $|\delta|\ll 1$. We have that

$f(u_{\star},v_{\star},\mathcal{F}(u_{\star},v_{\star}))=g(u_{\star},v_{\star},\mathcal{G}(u_{\star},v_{\star}))=0.$

Truncating at $O(\delta)$, we thus find that perturbations evolve according to

$\begin{split}\frac{\partial\xi}{\partial
t}&=\frac{\epsilon^{2}}{L^{2}}\frac{\partial^{2}\xi}{\partial
x^{2}}+\xi\frac{\partial f}{\partial u}+\eta\frac{\partial f}{\partial
v}+\frac{\partial
f}{\partial\mathcal{F}}\int_{\tau_{1}}^{\tau_{2}}K(s;\bm{p})\left[\xi_{s}\frac{\partial
F}{\partial u}+\eta_{s}\frac{\partial F}{\partial v}\right]\text{d}s,\\\
\frac{\partial\eta}{\partial
t}&=\frac{1}{L^{2}}\frac{\partial^{2}\eta}{\partial x^{2}}+\xi\frac{\partial
g}{\partial u}+\eta\frac{\partial g}{\partial v}+\frac{\partial
g}{\partial\mathcal{G}}\int_{\tau_{1}}^{\tau_{2}}K(s;\bm{p})\left[\xi_{s}\frac{\partial
G}{\partial u}+\eta_{s}\frac{\partial G}{\partial
v}\right]\text{d}s,\end{split}$ (17)

where all partial derivatives are evaluated at the steady state, and
$\xi_{s}=\xi(x,t-s)$, $\eta_{s}=\eta(x,t-s)$. Substituting into (17) an ansatz
of the form $(\xi,\eta)^{T}=e^{\lambda_{k}t}\cos(k\pi
x)(\xi_{0},\eta_{0})^{T}$ yields a homogeneous linear system for
$(\xi_{0},\eta_{0})^{T}$ given by

$\underbrace{\begin{pmatrix}\lambda_{k}+\frac{\epsilon^{2}}{L^{2}}(k\pi)^{2}-\frac{\partial
f}{\partial u}-\frac{\partial f}{\partial\mathcal{F}}\frac{\partial
F}{\partial u}E_{k}(\lambda_{k})&-\frac{\partial f}{\partial v}-\frac{\partial
f}{\partial F}\frac{\partial F}{\partial v}E_{k}(\lambda_{k})\\\
-\frac{\partial g}{\partial u}-\frac{\partial g}{\partial G}\frac{\partial
G}{\partial
u}E_{k}(\lambda_{k})&\lambda_{k}+\frac{1}{L^{2}}(k\pi)^{2}-\frac{\partial
g}{\partial v}-\frac{\partial g}{\partial\mathcal{G}}\frac{\partial
G}{\partial
v}E_{k}(\lambda_{k})\end{pmatrix}}_{\bm{M}}\begin{pmatrix}\xi_{0}\\\
\eta_{0}\end{pmatrix}=\begin{pmatrix}0\\\ 0\end{pmatrix},$ (18)

with

$E_{k}(\lambda_{k})=\int_{\tau_{1}}^{\tau_{2}}K(s;\bm{p})e^{-\lambda_{k}s}\
\text{d}s.$

Looking for non-trivial solutions of this system, we set $\det(\bm{M})=0$ to
compute values of $\lambda_{k}$. This leads to a characteristic equation (or
dispersion relation) of the form,

$\mathcal{D}_{k}=\lambda_{k}^{2}+\alpha_{k}\lambda_{k}+\beta_{k}+(\gamma_{k}\lambda_{k}+\delta_{k})E_{k}(\lambda_{k})+\chi_{k}E^{2}_{k}(\lambda_{k})=0,$
(19)

with the coefficients given as

$\begin{split}\alpha_{k}&=\left(\frac{\epsilon^{2}}{L^{2}}+\frac{1}{L^{2}}\right)(k\pi)^{2}-\frac{\partial
f}{\partial u}-\frac{\partial g}{\partial v},\\\
\beta_{k}&=\left(\frac{\epsilon^{2}}{L^{2}}(k\pi)^{2}-\frac{\partial
f}{\partial
u}\right)\left(\frac{1}{L^{2}}(k\pi)^{2}-\frac{\mathrm{d}g}{\mathrm{d}v}\right)-\frac{\partial
f}{\partial v}\frac{\partial g}{\partial u},\\\ \gamma_{k}&=-\frac{\partial
g}{\partial\mathcal{G}}\frac{\partial\mathcal{G}}{\partial v}-\frac{\partial
f}{\partial\mathcal{F}}\frac{\partial\mathcal{F}}{\partial u},\\\
\delta_{k}&=-\frac{\partial
g}{\partial\mathcal{G}}\frac{\partial\mathcal{G}}{\partial
v}\left(\frac{\epsilon^{2}}{L^{2}}(k\pi)^{2}-\frac{\partial f}{\partial
u}\right)-\frac{\partial
f}{\partial\mathcal{F}}\frac{\partial\mathcal{F}}{\partial
u}\left(\frac{1}{L^{2}}(k\pi)^{2}-\frac{\partial g}{\partial
v}\right)-\frac{\partial f}{\partial v}\frac{\partial g}{\partial
G}\frac{\partial G}{\partial u}-\frac{\partial g}{\partial u}\frac{\partial
f}{\partial F}\frac{\partial F}{\partial v},\\\ \chi_{k}&=\frac{\partial
f}{\partial\mathcal{F}}\frac{\partial\mathcal{F}}{\partial u}\frac{\partial
g}{\partial\mathcal{G}}\frac{\partial\mathcal{G}}{\partial v}-\frac{\partial
f}{\partial F}\frac{\partial F}{\partial v}\frac{\partial g}{\partial
G}\frac{\partial G}{\partial u}.\end{split}$ (20)

The coefficients $\alpha_{k}$, $\beta_{k}$, $\gamma_{k}$, $\delta_{k}$, and
$\chi_{k}$ depend only on the intrinsic model parameters, which are
independent of the time delay distribution and presented for the LI and the
two GM models, GM1 and GM2, in table 1.

Hence, one may note that the only term in the dispersion relation that varies
with the distribution is $E_{k}(\lambda_{k})$, which can be computed for a
skewed Gaussian distribution via:

$E_{k}=\int_{\tau_{1}}^{\tau_{2}}K_{S}(s;\mu,\omega,\rho)e^{-\lambda_{k}s}\
\text{d}s=\frac{\Psi_{c}}{\omega\sqrt{2\pi}}\bigintsss_{\tau_{1}}^{\tau_{2}}\left(1+\text{erf}\left(\rho\frac{s-\mu}{\omega\sqrt{2}}\right)\right)\exp\left(-\frac{1}{2}\left(\frac{s-\mu}{\omega}\right)^{2}-\lambda_{k}s\right)\text{d}s.$
(21)

Setting $\mu=\tau$, $\omega=\sigma$, and $\rho=0$ to consider a symmetric
distribution yields:

$E_{k}=\int_{\tau_{1}}^{\tau_{2}}K_{N}(s;\tau,\sigma)e^{-\lambda_{k}s}\
\text{d}s=\frac{\Phi_{c}}{2}\left[\exp\left(\frac{\lambda_{k}(\lambda_{k}\sigma^{2}-2\tau)}{2}\right)\text{erf}\left(\frac{\lambda_{k}\sigma^{2}+s-\tau}{\sqrt{2}\sigma}\right)\right]\Bigg{|}_{\tau_{1}}^{\tau_{2}}.$
(22)

Finally, taking $\sigma\to 0$ results in

$E_{k}=e^{-\lambda_{k}\tau},$ (23)

as expected, corresponding to the fixed delay case with time delay $\tau$.

The $E_{k}$ term represents the main impact of delay on linear stability. Even
in the relatively simple fixed delay case, we note that the (19) is not a
quadratic equation for $\lambda_{k}$, but a transcendental one. The more
complicated forms of $E_{k}$ in the distributed cases can be evaluated
numerically, which is in general necessary to analyze these transcendental
dispersion relations. We also remark that while (21) seemingly contains many
parameters, all of these can be rewritten in terms of $\tau$ and $\sigma$, and
hence compared directly with the symmetric case given in (22) (though the
expressions for $\rho\neq 0$ provide no obvious insight, and again one must
resort to numerical computation of $\lambda_{k}$).

The characteristic equation (19) can be used to determine the parameter sets
$(a,b,\epsilon^{2},L^{2},\tau,\sigma)$ in which a Turing instability occurs
(the ‘Turing space’) and hence where we expect pattern formation. If there
exists a $k\neq 0$ for a given set of parameters such that
$\max_{k}(\Re(\lambda_{k}))>0$, then we expect pattern formation. The largest
value of $\Re(\lambda_{k})$ also gives some indication of how quickly a
perturbation grows away from the homogeneous steady state, and hence gives a
heuristic estimate of the time taken for pattern formation to start.

| LI | GM1 | GM2
---|---|---|---
Parameter |  |  |
$\alpha_{k}$ | $\left(\frac{\epsilon^{2}}{L^{2}}+\frac{1}{L^{2}}\right)k^{2}\pi^{2}+u_{\star}^{2}+4u_{\star}v_{\star}+1$ | $\left(\frac{\epsilon^{2}}{L^{2}}+\frac{1}{L^{2}}\right)k^{2}\pi^{2}+b+1$ | $\left(\frac{\epsilon^{2}}{L^{2}}+\frac{1}{L^{2}}\right)k^{2}\pi^{2}+b+1$
$\beta_{k}$ | $\left(\frac{1}{L^{2}}\pi^{2}k^{2}+u_{\star}^{2}\right)\left(\frac{\epsilon^{2}}{L^{2}}\pi^{2}k^{2}+4u_{\star}v_{\star}+1\right)-4u_{\star}^{3}v_{\star}$ | $\left(\frac{\epsilon^{2}}{L^{2}}\pi^{2}k^{2}+b\right)\left(\frac{1}{L^{2}}\pi^{2}k^{2}+1\right)$ | $\left(\frac{\epsilon^{2}}{L^{2}}\pi^{2}k^{2}+b\right)\left(\frac{1}{L^{2}}\pi^{2}k^{2}+1\right)$
$\gamma_{k}$ | $-6u_{\star}v_{\star}$ | $-2\frac{u_{\star}}{v_{\star}}$ | $-2\frac{u_{\star}}{v_{\star}}$
$\delta_{k}$ | $-\frac{6}{L^{2}}u_{\star}v_{\star}k^{2}\pi^{2}$ | $-2\frac{u_{\star}}{v_{\star}}\left(\frac{1}{L^{2}}k^{2}\pi^{2}+1\right)$ | $-2\frac{u_{\star}}{v_{\star}}\left(\frac{1}{L^{2}}k^{2}\pi^{2}+1\right)+2\frac{u_{\star}^{3}}{v_{\star}^{2}}$
$\chi_{k}$ | 0 | $2\frac{u_{\star}^{3}}{v_{\star}^{2}}$ | 0
Table 1: A table of the coefficients of the characteristic equation (19) for
the LI, GM1, and GM2 models. We note that these coefficients agree with those
computed in yigaffneyli (52) for the LI model, up to a small typographical
error in the coefficient of $\beta_{k}$.

##  4 Turing Spaces Under Delay

In this Section we apply the results of the linear instability analysis to
produce bifurcation diagrams of Turing pattern formation as the form of delay
varies. Our results highlight the importance of the placement of time delay
terms in the reaction kinetics, and show that this can alter the effect that
delay has on the Turing space. We finally examine the effects that a
continuously distributed delay has on the dominant eigenvalue of perturbation
growth compared with that of a fixed delay with the same mean delay.

###  4.1 Fixed Time Delay

(a) Turing space for Schnakenberg model.

(b) Turing space for GM model.

Figure 2: Turing spaces for the Schnakenberg and GM models without delay, with
$\epsilon^{2}=0.001$ and $L^{2}=4.5$. The dashed black arc (leftmost boundary)
corresponds to parameters in which $\Re(\lambda_{0})=0$, and the solid black
outer arc (rightmost boundary) corresponds to parameters satisfying
$\max_{k}(\Re(\lambda_{k}))=0$ for $k\neq 0$.

In Figure 2 we show Turing spaces of the Schnakenberg and GM models, which are
bifurcation diagrams indicating regions of Turing instability in the absence
of delay. We note two separate curves which separate the parameter space into
its distinct regions. These will be referred to as the stability lines, and
are highlighted for the models in red and green. The inner green arc
corresponds to the values of $(a,b)$ such that $\Re(\lambda_{0})=0$ for the
spatially homogeneous characteristic equation. The outer red boundary is
comprised of the points $(a,b)$ such that $\max_{k}(\Re(\lambda_{k}))=0$,
$k\neq 0$. We are interested in how these Turing spaces change as time delay
is varied, and are also interested in the quantitative effects that changing
delay has on the value of $\max_{k}(\Re(\lambda_{k}))$, as this can be thought
of as a proxy for the time taken for patterns to form.

The roots of the characteristic equation were found using the _roots_ command
of the MATLAB package Chebfun chebfun (8). In order to compute
$\max_{k}(\Re(\lambda_{k}))$, we varied $k\in\mathbb{Z}$ over $[0,50]$ for a
given $\tau$ and then took the maximum over $k$. We do not consider $k>50$ as
full numerical solutions for the parameter values used tended towards patterns
with four ‘spikes’, so do not expect much larger wavenumbers to be excited
(and full numerical simulations were used to check this assumption throughout,
as well as to check the predictions of pattern formation).

####  4.1.1 The LI Model

(a) $\tau=0$.

(b) $\tau=1.5$.

Figure 3: $\max_{k}(\Re(\lambda_{k}))$ computed over $(a,b)$ parameter space
by solving (19) for the fixed time delay LI model (14), with
$\epsilon^{2}=0.001$, $L^{2}=4.5$. As $\tau$ increases,
$|\max_{k}(\Re(\lambda_{k}))|$ decreases. Stability lines for
$\Re(\lambda_{0})=0$ and $\max_{k}(\Re(\lambda_{k}))=0$, $k\neq 0$, are
overlaid, indicating the Turing space between them.

Considering first the LI model with fixed time delay, in Figure 3 we plot a
heatmap of $\max_{k}(\Re(\lambda_{k}))$ over the $(a,b)$ parameter space for
the two fixed time delays $\tau=0$, $\tau=1.5$. We overlay contour lines
corresponding to where $\Re(\lambda_{0})=0$ and
$\max_{k}(\Re(\lambda_{k}))=0$, $k\neq 0$, highlighting the Turing instability
region. As $\tau$ increases, the region of homogeneous instability (leftmost
curve in Figure 3 corresponding to $\Re(\lambda_{0})=0$) decreases, increasing
the Turing space. It can be seen that the absolute value
$|\max_{k}(\Re(\lambda_{k}))|$ also decreases for $\tau=1.5$. This suggests
pattern formation will take longer to occur, but this will happen over a
larger Turing instability region of the $(a,b)$ parameter space. It similarly
suggests that for $(a,b)$ such that $\max_{k}(\Re(\lambda_{k}))<0$, it will
take a longer time for the eigenfunctions with modes $k\neq 0$ to decay to a
spatially homogeneous steady state. These results have been verified through
full numerical solutions. Furthermore the outer curve, corresponding to
$\max_{k}(\Re(\lambda_{k}))=0$ for $k\neq 0$, does not move at the resolution
of the plotting $\tau$ changes, while analogous results (not shown) were found
for other values of $\tau$ and $\epsilon$ sargood2022gene (41).

These results were confirmed through full numerical solutions, and pattern
formation was found where suggested by linear theory. We note that care ought
to be taken when numerically simulating these models. The parameter region in
the bottom left of the parameter space is a delicate region that can exhibit
both Turing and (homogeneous) Hopf bifurcations, leading to complex spatio-
temporal behaviours. This type of dynamics in reaction-diffusion systems has
been studied more extensively in krausefixed (40, 20) and many others.

####  4.1.2 The GM Models

The results in fadai1 (12) showed that an increasing time delay in GM1 had an
antagonistic effect, shrinking the parameter space exhibiting stable spike
solutions. In contrast, an increasing time delay for GM2 caused an expansion
of the stable spike solution parameter regime fadai2 (11). Here, we use our
linear analysis of the spatially homogeneous steady states to examine how an
increasing time delay will affect the Turing space for each of these variants.

Results in Figure 4(b) show a shrinking Turing space for increasing $\tau$ for
GM1, consistent with the analysis conducted in fadai1 (12), which showed a de-
stabilisation of the stable spike solution parameter space with an increasing
$\tau$. Similarly, results in fadai2 (11) showed a stabilising effect of
increasing $\tau$ on the stable spike solution parameter space for model GM2
in (16), and we find an analogous result for the Turing space, shown in Figure
4(c). In both cases we also observe a decrease in
$|\max_{k}(\Re(\lambda_{k}))|$ as the delay is increased. We remark that the
choice of delay times $\tau$ for these and other plots using the dispersion
relation (19) was sufficiently large to demonstrate key trends for increasing
$\tau$, but not so large that there were convergence issues in solving the
transcendental dispersion relation (due to the multiple-scale nature of
polynomial and exponential rootfinding).

(a) $\tau=0$

(b) $\tau=0.5$ GM1 (15)

(c) $\tau=0.8$ GM2 (16)

Figure 4: The maximum growth rate, $\max_{k}(\Re(\lambda_{k}))$, corresponding
to solutions of (19) for models (15) in (b) and (16) in (c) with a fixed time
delay, plotted for $(a,b)\in[0,1]\times[0,4]$, with varying $\tau$. Note that
these models are equivalent for $\tau=0$ in (a). Parameters
$\epsilon^{2}=0.001$ and $L^{2}=4.5$ used. Stability lines for
$\Re(\lambda_{0})=0$ and $\max_{k}(\Re(\lambda_{k}))=0$, $k\neq 0$, are
overlaid, indicating the Turing space between them.

These results were also checked through full numerical solutions, and we found
that the linear theory in all cases successfully predicted where Turing
patterns would or would not form. These results indicate the importance of the
positioning of time delayed terms within the kinetic reactions, as these
differences in positioning of time delayed terms shows both expansion and
contraction of the associated the Turing spaces. We finally note that in all
the models considered, the ligand internalisation (LI) model and the two
Gierer-Meinhardt (GM) models, increasing $\tau$ changes the size of the Turing
space only via the curve associated with the homogeneous characteristic
equation ($k=0$), suggesting a specific mechanism by which gene expression
time delays impact Turing space size.

###  4.2 A Symmetric Gaussian Distributed Delay

Next we consider the LI model with a symmetric Gaussian distribution for the
time delay associated with the modelling of gene expression. The parameter
$\sigma$ corresponds to a measure of the width of this distribution, with the
limit $\sigma\to 0$ capturing the fixed time delay case. In particular
$\sigma_{\max}$ has been defined to be the largest feasible value of $\sigma$
to avoid negative delays, as described in subsection 2.1. We consider the
absolute difference of $\max_{k}(\Re(\lambda_{k}))$ for varying $\sigma$
values as a fraction of $\sigma_{\max}$, for multiple $\tau$ and $\epsilon$
values. For each $(\tau,\epsilon^{2})$, bifurcation plots analogous to those
of the previous subsection have been computed for the distributed delay case
with varying $\sigma\in\\{\sigma_{\max}\times 0.99,\sigma_{\max}\times
0.2,\sigma_{\max}\times 0.1\\}$. However, the differences are below plotting
resolution when compared to the fixed delay case, as may be found in
sargood2022gene (41); hence for each $(\tau,\epsilon^{2})$ we instead tabulate
the absolute difference of $\max_{k}(\Re(\lambda_{k}))$ between each
distributed delay case and the fixed delay case, across the $(a,b)$ parameter
space as summarised in Table 2.

The largest absolute difference in $\max_{k}(\Re(\lambda_{k}))$ in Table 2 for
all $\sigma$, $\tau$ and $\epsilon^{2}$ considered across the parameter space
$(a,b)\in[0,1.4]\times[0,2]$ is $O(10^{-3})$. We therefore expect that for all
$(a,b)\in[0,1.4]\times[0,2]$, using a symmetric Gaussian distribution centred
at some mean $\tau$ (for small $\tau$) will not significantly affect Turing
instabilities compared to the fixed delay case, independent of the standard
deviation $\sigma$ of the distribution. These results are a fairly robust
indication that, at the level of linear instability, the distributed delay
does not appreciably change the impact of time delay on Turing instability
analysis.

Parameters Used | $\sigma_{\max}\times 0.99$ | $\sigma_{\max}\times 0.2\ $ | $\sigma_{\max}\times 0.1\ $
---|---|---|---
$\epsilon^{2}=0.001$ | $\tau=0.2$ | $0.0010$ | $4.2\times 10^{-5}$ | $1.1\times 10^{-5}$
$\epsilon^{2}=0.001$ | $\tau=1.0$ | $0.0078$ | $3.3\times 10^{-4}$ | $8.2\times 10^{-5}$
$\epsilon^{2}=0.01$ | $\tau=0.2$ | $0.0025$ | $9.4\times 10^{-5}$ | $2.3\times 10^{-5}$
$\epsilon^{2}=0.01$ | $\tau=0.5$ | $0.0076$ | $2.6\times 10^{-4}$ | $6.4\times 10^{-5}$
Table 2: Table showing the maximum value over the parameter spaces $(a,b)$ of
the absolute difference of $\max_{k}(\Re(\lambda_{k}))$ between distributed
delay cases and the fixed delay case, across the $(a,b)\in[0,1.4]\times[0,2]$
parameter space, for multiple $\tau$ and $\epsilon^{2}$ values. Throughout we
set $L^{2}=4.5$. Results are displayed to $2$ significant figures.

###  4.3 A Skewed Gaussian Distributed Delay

Evaluating the dispersion relation (19) in the skewed distribution case is
substantially more costly computationally, especially for larger values of
$\tau$ where root-finding becomes numerically difficult. Instead, we plot
dispersion relations for certain fixed $(a,b)$ parameter values close to the
Turing space boundaries of Figure 2(a), and compare these plots to the fixed
delay case. We show that for a small mean $\tau$, the skew, positive or
negative, does not significantly affect the value of
$\max_{k}(\Re(\lambda_{k}))$. We highlight here again that equation (12) can
be solved implicitly for $\mu(\tau)$, for a given $\tau$, using the fzero
command in MATLAB. For a given $\rho$, and each found $\mu$, we compute
$\max_{k}(\Re(\lambda_{k}))$ by solving for roots of the characteristic
equation (19). In Figure 5, we plot $\max_{k}(\Re(\lambda_{k}))$ against
$\tau\in[0,0.8]$ for skew parameter values of $\rho=-10,10$, and
$\omega=\omega_{\max}\times 0.99$, with two different $(a,b)$ parameter sets.
The value $\omega_{\max}$ is defined analogously to that of $\sigma_{\max}$. A
plot of $\max_{k}(\Re(\lambda_{k}))$ for the fixed delay case is also added
for comparison in each case.

(a) $(a,b)=(0.1,0.9)$.

(b) $(a,b)=(0.4,0.4)$.

Figure 5: Comparison of $\max_{k}(\Re(\lambda_{k}))$ computed from (19)
plotted against $\tau\in[0,0.8]$ for $\rho=-10,10$ against the fixed delay
case, with the distribution spread given by $\omega=\omega_{\max}\times 0.99$.
Parameter values $\epsilon^{2}=0.001$ and $L^{2}=4.5$ used. In generating the
plot, the time delay scale, $\tau$, was varied at regular intervals of $0.05$
with $k\in\mathbb{Z}$ ranging over $k\in[0,50]$ before a maximum was taken.

From Figure 5 we see that the curves differ slightly for small $\tau$, with
$\rho=-10$ having a slightly higher value of the maximum growth rate, and
$\rho=10$ a slightly lower value. The overall effect is very small despite the
large skew implemented in the distribution. Noting the overall unit scale of
$|\max_{k}(\Re(\lambda_{k}))|$ in Figs. 3, 4(b), 4(c), we anticipate that
these effects are small enough not to have a significant impact on the
timescale for patterning onset, as confirmed numerically via simulations of
the full model with distributed delay near the boundary of these
instabilities.

##  5 Numerical Exploration of the Models

To verify our linear instability results, as well as explore the dynamics of
pattern formation beyond the linear regime, full numerical solutions were
conducted. The spatial derivatives were discretized via the method-of-lines
using the standard three-point stencil for the Laplacian, with $m=500$ equally
spaced points on the domain $x\in\Omega=[0,1]$. This discretization results in
$m=500$ delay differential equations (DDEs) in time, which are solved via
built-in time-stepping solvers in Julia, typically _Rodas5_ , a 5-th order
A-stable solver, from the family of Rosenbrock methods rodas (36, 39). The
distributed delay terms were approximated via a composite Simpson’s rule using
$50$ quadrature points; see sargood2022gene (41) for a more thorough
discussion of the numerical procedures used. The Julia code used to generate
all numerical solutions throughout this paper can be found at the open source
repository git (42).

With a maximum delay of $\tau_{2}=\tau+3\sigma$ (or $\mu+3\omega$ with a
skewed-distribution), we note that to solve DDEs, a history function is
required to define the solution for $t\in[-\tau,0)$ for the fixed delay case,
$t\in[-\tau-3\sigma,0)$ for the symmetric distributed case, and
$t\in[-\mu-3\omega,0)$ for the skewed distributed case. Unless otherwise
stated, a constant history function equal to the initial conditions is used.
Finally, as the LI model has cross reaction kinetics, we have that when the
concentration of the activator $u$ is high, the concentration of the inhibitor
$v$ is low, and vice-versa murray (32). The concentration gradients of the two
morphogens $u$ and $v$ are thus effectively ‘out of phase’, and so it is
sufficient to show just the numerical solution of the activator $u$. Hence,
only the numerical solution of the activator $u$ is plotted.

###  5.1 A Relationship between Fixed Delay and Onset of Patterning

Here we show that, for small $\tau$ and $L^{2}$, the linear theory provides a
good approximation to the time taken until pattern formation occurs, and in
fact, the relationship between $\tau$ and time-to-pattern under these
conditions is linear for the fixed delay case. We also show that, through full
numerical solutions, the relationship between $\tau$ and time-to-pattern on a
longer timescale for larger $\tau$ is also linear. We first consider the
former.

(a) Dispersion curve plotted with domain size $L^{2}=0.2$. Curve produced by
varying $k\in[0,5]$ at regular intervals of $1$.

(b) Dispersion curve plotted with domain size $L^{2}=4.5$. Curve produced by
varying $k\in[0,60]$ at regular intervals of $1$.

Figure 6: Dispersion curves of the characteristic equation given in (19) for
the LI model plotted for $(a,b,\tau)=(0.4,1.8,0.2)$ and $\epsilon^{2}=0.001$.
Discrete values of $k$ overlaid as scatter points. A larger $L^{2}$ results in
more unstable modes $\lambda_{k}$ such that $\Re(\lambda_{k})>0$.

We restrict the domain size to $L^{2}=0.2$, as a smaller domain results in
fewer unstable modes and thus less competition for the dominant mode,
resulting in a better approximation of the linear theory. This finite size
effect can be seen in Figure 6, where $\Re(\lambda_{k})$ is plotted against
$k$ for two different domain sizes, for a given $(a,b,\tau)$. We also use an
initial small perturbation of the form of Eqn. (3) for $t\in[-\tau,0]$, but
vary the standard deviation, $\sigma_{IC}$, of these perturbations.

The linear theory suggests that the perturbation will be of the form
$\xi(x,t)\sim A_{k}(t)\cos(k\pi x)$, where $k$ is the dominant mode and
$A_{k}(t)$ denotes the corresponding Fourier coefficient at time $t$. For a
given parameter set $(a,b,\epsilon^{2},\tau,L^{2})$, we can solve the
characteristic equation (19), and plot $\Re(\lambda_{k})$ against $k$, to
determine the dominant mode $k$ and the corresponding growth rate,
$\lambda_{k}$. When the perturbation $\xi$ has grown sufficiently, in absolute
value, beyond a threshold where pattern formation is considered, we call this
time $t=T^{\star}$. More specifically, the time $T^{\star}$ is the first time
such that $\max_{x}|u(T^{\star},x)-u_{\star}|>u_{T^{\star}}$, where
$u_{T^{\star}}$ is a given threshold. This then gives the first time such that
any solution point across the whole spatial domain is large enough, in
absolute difference, from the steady state. We refer to this value $T^{\star}$
as the ‘true’ time-to-pattern. Finally, using the relation
$A_{k}(T^{\star})\approx A_{k}(0)e^{\lambda_{k}T^{\star}}$, we can rearrange
for $T^{\star}$ and thus compute a linear theory approximation for the ‘true’
time-to-pattern, via

$T^{\star}\approx\frac{1}{\lambda_{k}}\ln\left(\frac{A_{k}(T^{\star})}{A_{k}(0)}\right)=:T.$
(24)

The approximate time-to-pattern, $T$, can thus be determined via
$\lambda_{k}$, from (19), as well as $A_{k}(0)$ and $A_{k}(T^{\star})$,
computed via the Fast Fourier Transform.

(a) $(a,b)=(0.4,1.8)$.

(b) $(a,b)=(0.1,0.9)$.

(c) $(a,b)=(0.2,1.3)$.

Figure 7: Predicted vs ‘true’ time-to-pattern for numerical solutions of the
LI model with fixed time delay and boundary conditions given by (2). Initial
perturbations from the steady state are taken with
$\sigma_{\text{IC}}=10^{-5}$ and a threshold of $u_{T^{*}}=0.1$ is used to
compute the ‘true’ time-to-pattern, $T$. The predicted time-to-pattern is
computed using (24) with $\lambda_{k}$ computed via (19), for three different
parameter sets, with $L^{2}=0.2$, $\epsilon^{2}=0.001$, and $\tau\in[0,1.6]$,
varied at regular intervals of $0.2$.

We repeat the evaluation of both the true and predicted times to pattern for
varying $(a,b,\tau)$, and compare the evaluations in Figure 7, which shows the
comparison for three different parameter sets for $(a,b)$, with the time delay
varied over $\tau\in[0,1.6]$ at intervals of $0.2$. Other parameter values,
and other choices of the threshold $u_{T^{*}}$, gave qualitatively similar
results, essentially indicating that the computed value of $\lambda_{k}$ for
the dominant mode accurately predicts the time to pattern, which is observed
to vary linearly in $\tau$ to a very good approximationsargood2022gene (41).
We remark that computing $\lambda_{k}$ for $\tau>1.6$ is numerically
difficult; thus full numerical solutions were used to plot ‘true’ time-to-
pattern against $\tau$ on a longer timescale, for both the LI model and GM
models, to verify the linear relationship between the time delay and time to
pattern. In Figure 8 we consider the time delay $\tau\in[1,16]$ at unit
intervals for the LI model, with two different parameter sets of $(a,b)$ and
plot $T$, the time taken for a perturbation to grow up to a threshold value
$u_{T^{*}}=0.1$ from an initial perturbation of size
$\sigma_{\text{IC}}=10^{-5}$, which once more reveals a linearly increasing
relationship between $T$ and $\tau$.

For the GM1 model (15), as a result of the shrinking Turing space, we only
consider small $\tau\in[0.1,1]$, varied at regular intervals of $0.1$. For the
GM2 model (16), we consider both small $\tau\in[0.1,2]$ and larger
$\tau\in[1,16]$, varied at regular intervals of $0.1$ and $1$ respectively. We
set an initial perturbation from the steady state as
$\sigma_{\text{IC}}=0.001$, and a threshold value of $10$ (these larger values
are chosen for computational reasons, and qualitatively identical results are
found for other choices). ‘True’ time-to-pattern results, analogous to the LI
model, are presented in Figure 9, for both Gierer-Meinhardt models, Eqn. (15)
and Eqn. (16), where a linear relationship may be observed for both models.
Finally, we note that even though a linear relation is consistently observed
between the time delay and time to pattern, there is nonetheless a difference
in line slope between the different models, illustrating that kinetics can
affect the sensitivity in the timing of patterning onset with the time delay.
Nonetheless, in all cases shown, the slopes are generally steep, with delays
of $O(1)$ time units leading to an order-of-magnitude increase in the time-to-
pattern.

(a) $(a,b)=(0.1,0.9)$

(b) $(a,b)=(0.4,0.8)$

Figure 8: ‘True’ time-to-pattern for full numerical solutions of the fixed
time delay LI model plotted against $\tau\in[1,16]$ for
$\sigma_{\text{IC}}=10^{-5}$ and threshold $u_{T^{*}}=0.1$. Parameters used
are $\epsilon^{2}=0.001$ and $L^{2}=4.5$.

(a) ‘True’ time-to-pattern vs $\tau$ for $\tau\in[0.1,2]$ varied at intervals
of $0.1$. The blue line shows results for the model in (15), and the red line
for (16).

(b) ‘True’ time-to-pattern vs $\tau$ for $\tau\in[1,16]$ varied at intervals
of $1$. Model given in (16).

Figure 9: ‘True’ time-to-pattern results for the two fixed time delay GM
variants given in (15) and (16). Initial random perturbations given with
$\sigma_{\text{IC}}=0.001$, and threshold value as $u_{T^{*}}=10$. Parameters
$(a,b)=(0.75,0.5)$, $\epsilon^{2}=0.001$ and $L^{2}=4.5$ are used. Boundary
conditions set as in (2).

###  5.2 Onset of Patterning in Models with Distributed Delay

We have shown the linear stability analysis predicts that a continuously
distributed time delay with a symmetric or skewed Gaussian kernel has only a
small (and often negligible) influence on the value of
$\max_{k}(\Re(\lambda_{k}))$, regardless of the distribution properties,
compared to a fixed delay model with the same mean delay. We thus look to
confirm this result directly via full numerical simulations of the distributed
delay cases.

We first consider Figures 10 and 11, which show numerical solutions for the LI
model with a symmetric Gaussian distributed delay, using $(a,b)=(0.1,0.9)$ for
$\tau\in\\{1,16\\}$ and varying $\sigma$, and compared with the appropriate
fixed delay case given an initial small perturbation of the form of Eqn. (3)
for $t\in[-\tau,0]$. The results illustrate that the onset of patterning, and
the emergent pattern do not depend on the value of $\sigma$ used at the
resolution of the plotting.

(a) Fixed delay LI model.

(b) Distributed delay LI model with $\sigma=\sigma_{\max}\times 0.99$.

(c) Distributed delay LI model with $\sigma=\sigma_{\max}\times 0.1$.

Figure 10: Numerical simulations of (14) showing a comparison of the fixed vs
symmetric distributed delay case for $\tau=1$. Boundary conditions are given
by (2) and initial conditions by (3), with parameters $(a,b)=(0.1,0.9)$,
$\epsilon^{2}=0.001$, and $L^{2}=4.5$.

(a) Fixed delay LI model.

(b) Distributed delay LI model with $\sigma=\sigma_{\max}\times 0.99$.

(c) Distributed delay LI model with $\sigma=\sigma_{\max}\times 0.1$.

Figure 11: Numerical simulations of (14) showing a comparison of the fixed vs
symmetric distributed delay case for $\tau=16$. Boundary conditions are given
by (2) and initial conditions by (3), with parameters $(a,b)=(0.1,0.9)$,
$\epsilon^{2}=0.001$, and $L^{2}=4.5$.

We next present numerical results to show that even the skew of an asymmetric
distribution for both a small and large mean delay $\tau$, does not have a
significant effect on observable patterning results, as highlighted for the LI
model by Figure 12 for $\tau=0.1$ and Figure 13 for $\tau=16$. Despite a range
of skews under consideration (see panel (a) in both Figures), the simulations
are indistinguishable by eye to one another and to the fixed delay simulations
given in panel (b) in each case. We finally remark that additional simulations
in other parameter regimes, and for the Gierer-Meinhardt models, illustrate
the same independence of the skew, with further examples can be found in
sargood2022gene (41).

(a) Probability density functions of skewed truncated Gaussian distributions,
with $\rho=-10,10$ and mean $\tau=0.1$.

(b) Fixed delay case with $\tau=0.1$.

(c) Numerical simulation with skewed distribution of $\rho=-10$. Distribution
parameters are $\mu=0.124(3s.f.)$ and $\omega=0.0408(3s.f.)$.

(d) Numerical simulation with skewed distribution of $\rho=10$. Distribution
parameters are $\mu=0.0863(3s.f.)$ and $\omega=0.0285(3s.f.)$.

Figure 12: Simulations of (14) with fixed and distributed delay with
parameters $(a,b)=(0.1,0.9)$, $\epsilon^{2}=0.01$ and $L^{2}=4.5$,
$\rho=-10,10$ and $\tau=0.1$. Initial conditions are given by (3) and boundary
conditions by (2).

(a) Probability density functions of skewed truncated Gaussian distributions,
with $\rho=-10,10$, and mean $\tau=16$.

(b) Fixed delay case with $\tau=16$.

(c) Numerical simulation with distribution of $\rho=-10$. Distribution
parameters are $\mu=21.7(3s.f.)$ and $\omega=7.16(3s.f.)$.

(d) Numerical simulation with distribution of $\rho=10$. Distribution
parameters are $\mu=12.7(3s.f.)$ and $\omega=4.18(3s.f.)$.

Figure 13: Simulations of (14) with fixed and distributed delay with
parameters $(a,b)=(0.1,0.9)$, $\epsilon^{2}=0.01$ and $L^{2}=4.5$,
$\rho=-10,10$ and $\tau=16$. Initial conditions are given by (3) and boundary
conditions by (2).

###  5.3 Boundary and Initial Conditions

Figure 14: The functions $u_{1}(x)$, $v_{1}(x)$ used to generate the initial
conditions IC1; algebraic expressions for these functions are given in the
Appendices of gaffmonk (14).

We next consider the robustness of our results under variation of the initial
and boundary conditions, focusing on the LI model with a fixed delay for
simplicity. We first consider the sensitivity of pattern formation in the
context of a fixed time delay to varying initial conditions. Fixing parameters
by

$(a,b,\epsilon^{2},L^{2})=(0.1,0.9,0.001,4.5)$

we proceed to consider three different sets of initial conditions. We denote
by $\text{IC}_{1}$ initial conditions based on those used in gaffmonk (14).
With $u_{1}(x)$ and $v_{1}(x)$ the functions depicted in Fig. 14, these
conditions are given by

$u(x,t)=u_{1}(x),\leavevmode\nobreak\ \leavevmode\nobreak\ v(x,t)=v_{1}(x),$
(25)

for $t\in[-\tau,0]$ and represent a specific perturbation of the homogeneous
steady state. We denote by $\text{IC}_{2}$ the initial conditions defined in
Eqn. (3) with $\sigma_{IC}=0.01$, and $\text{IC}_{3}$ denote the same initial
conditions, but with $\sigma_{IC}=0.1$. As previously, both initial conditions
were specified for $t\in[-\tau,0]$ and a fixed random seed was set for the
random variable, $r$, and hence the perturbations $r_{u}(x)$ and $r_{v}(x)$ of
Eqn. (3) are the same in each simulation.

(a) $\text{IC}_{1}$ given by Eqn. (25).

(b) $\text{IC}_{2}$ given by Eqn. (3), with $\sigma_{IC}=0.01$.

(c) $\text{IC}_{3}$ given by Eqn. (3), with $\sigma_{IC}=0.1$.

Figure 15: Simulations of LI model (14) for varying ICs and $\tau=1$. Boundary
conditions given by (2) with parameters $(a,b)=(0.1,0.9)$,
$\epsilon^{2}=0.001$, and $L^{2}=4.5$.

(a) $\text{IC}_{1}$ given by Eqn. (25).

(b) $\text{IC}_{2}$ given by equation (3), with $\sigma_{IC}=0.01$.

(c) $\text{IC}_{3}$ given by Eqn. (3), with $\sigma_{IC}=0.1$.

Figure 16: Simulations of LI model (14) for varying ICs and $\tau=16$.
Boundary conditions given by (2) with parameters $(a,b)=(0.1,0.9)$,
$\epsilon^{2}=0.001$, and $L^{2}=4.5$.

Figures 15 and 16 show simulations for each of these initial conditions for
varying fixed time delay $\tau\in\\{1,16\\}$. The final pattern is somewhat
sensitive to the choice of initial conditions. As one might expect, the larger
$\sigma_{\text{IC}}$ used in $\text{IC}_{3}$, compared to that of
$\text{IC}_{2}$, results in a faster onset of pattern formation, as does using
a random perturbation rather than one along a single mode. Although the time
taken until onset of patterning varies with different initial conditions, the
increase in time-to-pattern with an increasing time delay is consistent
independent of the initial conditions chosen. By considering the varying
$x$-axis, we also note that in each case, this relationship appears to be
linear, as predicted earlier.

The effects of time-dependent histories were also considered. Taking initial
data to be the homogeneous steady states multiplied by $1+r\sin(wt)$, with
varying real $w$ and $r$ Normally distributed at each spatial point, had only
transient effects with no noticeable impact on the long-time pattern structure
or time to pattern results. See sargood2022gene (41) for examples.

Finally, we consider the effect of varying boundary conditions. Motivated by
the analysis in krausemixed (24), homogeneous Dirichlet boundary conditions
are implemented for the activator term, and homogeneous Neumann boundary
conditions for the inhibitor term. Hence we set

$u=\frac{\partial v}{\partial x}=0\quad x=0,1.$ (26)

The results in Figures 17 were generated using $\text{IC}_{2}$, with time
delays of $\tau=1$ and $\tau=16$, respectively. These show simulations
generated with homogeneous Neumann conditions for both $u$ and $v$, as in the
boundary conditions (2), and those generated with the mixed conditions (26).
While these mixed conditions changed the kind and structure of pattern
observed, as described in krausemixed (24) as ‘isolated’ patterns, we do not
observe a noticeable difference in the time taken for patterns to form.

(a) $\tau=1$, Homogeneous Neumann boundary conditions, as given by Eqn. (2).

(b) $\tau=1$, Mixed boundary conditions, as given by Eqn. (26).

(c) $\tau=16$, Homogeneous Neumann boundary conditions, as given by Eqn. (2).

(d) $\tau=16$, Mixed boundary conditions, as given by Eqn. (26).

Figure 17: Comparison of varying BCs for simulations of the LI model (14) with
fixed delay and parameters $(a,b)=(0.1,0.9)$, $\epsilon^{2}=0.001$,
$L^{2}=4.5$, with the initial conditions specified by Eqn. (3). In the top row
the time delay is fixed at $\tau=1$, and in the bottom row $\tau=16.$

##  6 Discussion

We have studied reaction-diffusion systems gene expression delays modelled as
both discrete and continuous distributions of delay. We used linear stability
theory and systematic numerical simulation to explore the impact of varying
delay distribution, as well as initial and boundary data, on the pattern
formation process. Our results contribute to the clarification of aspects of
the impact of such delays from previous work, and shed light on which model
features are important for consideration in biological applications.

We have provided evidence, via both linear stability and full numerical
simulations, that the time until pattern onset linearly increases with
increasing mean delay. This was observed across different reaction kinetics
and the form of the delay, as well as initial and boundary data considered.
This suggests that the impact of time delay on slowing pattern formation
processes in reaction-diffusion systems is a general phenomenon which scales
linearly with time delay. Although the type of pattern seen can change with
these variations (particularly for changes in boundary conditions), the
increase in lag until onset of patterning as a result of time delay appears to
be robust. Importantly, for the models studied here, this increase in time-to-
pattern was substantial for even small and moderate delays, confirming earlier
predictions in gaffmonk (14). An important avenue of future work will be to
understand what leads to these particular slopes, and hence refine our
understanding of the plausibility of these models in the presence of gene
expression time delays for pattern forming systems.

For the Schnakenberg kinetics, where fixed gene expression delays were
motivated by ligand internalisation models, we have noted that increased time
delays act to expand the Turing space. This effect was displayed here through
the use of bifurcation diagrams, and was confirmed via numerical solutions.
Motivated by the stability analysis of spike solutions of the GM model in
fadai1 (12, 11), we demonstrated the importance of the positioning of time
delayed terms within a reaction-diffusion mechanism.

We further found increasing the time delay can either expand or contract the
Turing space and this was solely dependent on the boundary of the Turing space
provided by the stability of the spatially homogeneous mode in the presence of
delay. For the Schnakenberg ligand internalisation (LI) model this observation
may be explained by first noting from Appendix A.3.1 of Gaffney $\&$ Monk
gaffmonk (14) that a Hopf bifurcation cannot occur for a spatially
inhomogeneous mode with non-zero wavenumber for this system. Thus the growth
rate $\lambda_{k}$ is real at a bifurcation for $k\neq 0$ and hence at the
boundary of the Turing space dictated by the spatially inhomogeneous modes, we
have $\lambda_{k}=0.$ Given the time delay only occurs in the dispersion
relation, Eqn. (19), via terms of the form $\exp(-\lambda_{k}\tau)$ these
boundaries thus must be independent of $\tau.$ Thus if the Turing space is to
change for the LI model with changes in the time delay, it must be via the
boundary of the Turing space dictated by the spatially homogeneous mode, with
$k=0$ and a Hopf bifurcation on the boundary, which is generally possible, as
consistent with our observations. This further yields the question as to when
the change in the Turing space with increasing time delay is solely due to the
behavior of the spatially homogeneous mode, as also observed numerically for
the Gierer-Meinhardt (GM) models, in turn requiring a study of when a Hopf
bifurcation cannot occur for spatially inhomogeneous modes given a time delay.
It is an open question whether delay can induce Hopf bifurcations for
inhomogeneous modes in models beyond the ones studied here. In particular,
such bifurcations can occur in hyperbolic reaction-diffusion equations which
can have approximately the same linearization as delayed reaction-diffusion
systems for small delay ritchie_hyperbolic_2020 (38).

Finally, driven by the inherent stochasticity of the molecular processes
underpinning gene expression raj (37, 10, 28, 35), gene expression time delays
were modelled as both a symmetric and skewed Gaussian distribution. Through
linear analysis, and verified by numerical simulations, it was shown that the
distribution of delay has no qualitative (and negligible quantitative) impact
on solutions to these reaction-diffusion systems compared to a fixed delay
model with the same mean delay. Namely, solutions of the Schnakenberg model
seem to be dependent on the mean delay of the distribution used, but not on
the standard deviation or skew, and thus can effectively be modelled as purely
a fixed delay.

Our findings, that a distributed representation of time delay does not
alleviate the increased timescales of patterning events caused with a fixed
delay, reinforces earlier work indicating that such results are at odds with
rapid developmental biology patterning events gaffmonk (14). Of course, there
are still important limitations of our simple two-species models. It is
becoming increasingly clear that going beyond two-species models is crucial
for representing biological reality krausenew (22, 43, 45). Hence, one
potentially important avenue for further research would be to investigate the
effect of time delay on Turing mechanisms encapsulating larger systems. This
would aid in improving the possible applicability and similarity of Turing’s
models to more intricate biological dynamics.

Throughout this paper, we considered two different types of kinetics, namely
those based on Schnakenberg and GM models. Our results indicate some general
attributes that are common to both sets of kinetics. The first being a
linearly increasing relationship between incorporated time delay and time
until onset of patterning. The second being that the effect of time delay on
the Turing space is only dependent on the stability of the homogeneous
equilibrium in the absence of diffusion. A clear extension to these
observations would be to explore these effects for different reaction-
diffusion systems that can exhibit Turing patterns. Typical systems that could
be examined include the Gray-Scott grayscott (16) or Thomas murray (32)
models.

The use of other forms of distribution, such as the gamma or exponential
distributions could be considered, in order to verify our findings that, when
onset of patterning is being considered, the only relevant modelling parameter
required is the mean delay. We also only considered the problem on a one-
dimensional stationary spatial domain. Previous research has been conducted on
higher-dimensional spatial domains, and growing domains, with fixed delay
gaffmonk (14, 40). Although we hypothesise that our results with a distributed
delay will be consistent across variations in the spatial domain considered,
there is room to explore these possibilities. Finally, we note that due to
numerical limitations when using Chebfun to find roots of the transcendental
characteristic equations, the linear theory could only be tested for small-
time delays. We found that, for these small-time delays, the linear theory
generally provided a good approximation to the time-to-pattern, and all
conclusions from the linear theory could be verified using full numerical
solutions. Further work could therefore explore solving the characteristic
equations derived in this paper for larger time delay values, and examine
whether the linear theory still provides good approximations to the model
behaviour. The development of such techniques has applicability for other non-
classical dispersion relations, such as those derived in krause2020turing (23)
which have similar numerical difficulties due to the presence of exponential
functions.

###### Acknowledgements.

In the final stages of document preparation and review, A. S. was supported by
the UK Medical Research Council (G116768, SUAG/090) and The Wellcome Trust
Institutional Strategic Support Fund (RG89305, RG89529).

Data Statement In compliance with UKRI Research Council initiatives, all code
and data used to generate these Figures can be

found at http://dx.doi.org/xxx/xxx as well as at the open access repository
git (42).

## References

  * (1) Richard Bailleul, Marie Manceau and Jonathan Touboul  “A “Numerical Evo-Devo” Synthesis for the Identification of Pattern-Forming Factors”  In _Cells_ 9.8 Multidisciplinary Digital Publishing Institute, 2020, pp. 1840
  * (2) Jonathan BL Bard  “Morphogenesis: the cellular and molecular processes of developmental anatomy”  Cambridge University Press, 1992
  * (3) Dmitri Bratsun et al.  “Delay-induced stochastic oscillations in gene regulation”  In _Proceedings of the National Academy of Sciences_ 102.41 National Acad Sciences, 2005, pp. 14593–14598
  * (4) Tyler Cassidy  “Distributed delay differential equation representations of cyclic differential equations”  In _SIAM Journal on Applied Mathematics_ 81.4 SIAM, 2021, pp. 1742–1766
  * (5) James C Cavallo, Steffen Scholpp and Mark B Flegg  “Delay-driven oscillations via Axin2 feedback in the Wnt/$\beta$-catenin signalling pathway”  In _Journal of Theoretical Biology_ 507 Elsevier, 2020, pp. 110458
  * (6) Edmund J Crampin, Eamonn A Gaffney and Philip K Maini  “Reaction and Diffusion on Growing Domains: Scenarios for Robust Pattern Formation”  In _Bull. Math. Bio._ 61, 1999, pp. 1093–1120
  * (7) William Dash  “Distributed Delay in Reaction-Diffusion systems”  In _Mathematical Institute, University of Oxford_ , 2020
  * (8) Tobin A Driscoll, Nicholas Hale and Lloyd N Trefethen  “Chebfun guide”  Pafnuty Publications, Oxford, 2014
  * (9) AM Elaiw and AD Al Agha  “A reaction–diffusion model for oncolytic M1 virotherapy with distributed delays”  In _The European Physical Journal Plus_ 135.1 Springer, 2020, pp. 1–40
  * (10) Michael B Elowitz et al.  “Stochastic gene expression in a single cell”  In _Science_ 297.5584 American Association for the Advancement of Science, 2002, pp. 1183–1186
  * (11) Nabil T Fadai, Michael J Ward and Juncheng Wei  “A time-delay in the activator kinetics enhances the stability of a spike solution to the gierer-meinhardt model”  In _Discrete & Continuous Dynamical Systems-B_ 23.4 American Institute of Mathematical Sciences, 2018, pp. 1431
  * (12) Nabil T Fadai, Michael J Ward and Juncheng Wei  “Delayed Reaction Kinetics and the Stability of Spikes in the Gierer–Meinhardt Model”  In _SIAM Journal on Applied Mathematics_ 77.2 SIAM, 2017, pp. 664–696
  * (13) Cedric Flecher, Denis Allard and Philippe Naveau  “Truncated skew-normal distributions: moments, estimation by weighted moments and application to climatic data”  In _Metron_ 68.3 Springer, 2010, pp. 331–345
  * (14) EA Gaffney and NAM Monk  “Gene expression time delays and Turing pattern formation systems”  In _Bulletin of mathematical biology_ 68.1 Springer, 2006, pp. 99–130
  * (15) Alfred Gierer and Hans Meinhardt  “A theory of biological pattern formation”  In _Kybernetik_ 12.1 Springer, 1972, pp. 30–39
  * (16) Peter Gray and Stephen K Scott  “A new model for oscillatory behaviour in closed systems: the autocatalator”  In _Berichte der Bunsengesellschaft für physikalische Chemie_ 90.11 Wiley Online Library, 1986, pp. 985–996
  * (17) Matthew P Harris et al.  “Molecular evidence for an activator–inhibitor mechanism in development of embryonic feather branching”  In _Proceedings of the National Academy of Sciences_ 102.33 National Acad Sciences, 2005, pp. 11734–11739
  * (18) HGE Hentschel et al.  “Dynamical mechanisms for skeletal pattern formation in the vertebrate limb”  In _Proceedings of the Royal Society of London. Series B: Biological Sciences_ 271.1549 The Royal Society, 2004, pp. 1713–1722
  * (19) David Iron and Michael J Ward  “A metastable spike solution for a nonlocal reaction-diffusion model”  In _SIAM Journal on Applied Mathematics_ 60.3 SIAM, 2000, pp. 778–802
  * (20) Weihua Jiang, Hongbin Wang and Xun Cao  “Turing instability and Turing–Hopf bifurcation in diffusive Schnakenberg systems with gene expression time delay”  In _Journal of Dynamics and Differential Equations_ 31.4 Springer, 2019, pp. 2223–2247
  * (21) Charles B Kimmel et al.  “Stages of embryonic development of the zebrafish”  In _Developmental dynamics_ 203.3 Wiley Online Library, 1995, pp. 253–310
  * (22) Andrew L Krause et al.  “Modern Perspectives on Near-Equilibrium Analysis of Turing Systems”  In _arXiv preprint arXiv:2106.08375_ , 2021
  * (23) Andrew L Krause et al.  “Turing patterning in stratified domains”  In _Bulletin of Mathematical Biology_ 82.10 Springer, 2020, pp. 1–37
  * (24) Andrew L Krause et al.  “Isolating Patterns in Open Reaction–Diffusion Systems”  In _Bulletin of Mathematical Biology_ 83.7 Springer, 2021, pp. 1–35
  * (25) S Seirin Lee and EA Gaffney  “Aberrant behaviours of reaction diffusion self-organisation models on growing domains in the presence of gene expression time delays”  In _Bulletin of mathematical biology_ 72.8 Springer, 2010, pp. 2161–2179
  * (26) S Seirin Lee, EA Gaffney and NAM Monk  “The influence of gene expression time delays on Gierer–Meinhardt pattern formation systems”  In _Bulletin of mathematical biology_ 72.8 Springer, 2010, pp. 2139–2160
  * (27) Philip K Maini et al.  “Turing’s model for biological pattern formation and the robustness problem”  In _Interface focus_ 2.4 The Royal Society, 2012, pp. 487–496
  * (28) Harley H McAdams and Adam Arkin  “Stochastic mechanisms in gene expression”  In _Proceedings of the National Academy of Sciences_ 94.3 National Acad Sciences, 1997, pp. 814–819
  * (29) Takashi Miura and Philip K Maini  “Speed of pattern appearance in reaction-diffusion models: implications in the pattern formation of limb bud mesenchyme cells”  In _Bulletin of mathematical biology_ 66.4 Elsevier, 2004, pp. 627–649
  * (30) Takashi Miura and Kohei Shiota  “Extracellular matrix environment influences chondrogenic pattern formation in limb bud micromass culture: experimental verification of theoretical models”  In _The Anatomical Record: An Official Publication of the American Association of Anatomists_ 258.1 Wiley Online Library, 2000, pp. 100–107
  * (31) Takashi Miura and Kohei Shiota  “TGF$\beta$2 acts as an “Activator” molecule in reaction-diffusion model and is involved in cell sorting phenomenon in mouse limb micromass culture”  In _Developmental dynamics: an official publication of the American Association of Anatomists_ 217.3 Wiley Online Library, 2000, pp. 241–249
  * (32) James D Murray  “Mathematical biology II: spatial models and biomedical applications”  Springer-Verlag, 2001
  * (33) Kevin J Painter, Philip K Maini and Hans G Othmer  “Complex spatial patterns in a hybrid chemotaxis reaction-diffusion model”  In _J. Math. Biol_ 41.4, 2000, pp. 285–314
  * (34) Mike Patefield and David Tandy  “Fast and accurate calculation of Owen’s T function”  In _Journal of Statistical Software_ 5.5 American Statistical Association, 2000, pp. 1–25
  * (35) Johan Paulsson  “Models of stochastic gene expression”  In _Physics of life reviews_ 2.2 Elsevier, 2005, pp. 157–175
  * (36) Christopher Rackauckas and Qing Nie  “Differentialequations.jl–a performant and feature-rich ecosystem for solving differential equations in Julia”  In _Journal of Open Research Software_ 5.1 Ubiquity Press, 2017
  * (37) Arjun Raj and Alexander Van Oudenaarden  “Nature, nurture, or chance: stochastic gene expression and its consequences”  In _Cell_ 135.2 Elsevier, 2008, pp. 216–226
  * (38) J.. Ritchie, A.. Krause and R.. Van Gorder  “Turing and wave instabilities in hyperbolic reaction-diffusion systems: The role of second-order time derivatives and cross-diffusion terms on pattern formation” arxiv.org/abs/2204.13820  In _arXiv preprint_ , 2022
  * (39) HH Rosenbrock  “Some general implicit processes for the numerical solution of differential equations”  In _The Computer Journal_ 5.4 Oxford University Press, 1963, pp. 329–330
  * (40) Faustino Sánchez-Garduno et al.  “Turing–Hopf patterns on growing domains: the torus and the sphere”  In _Journal of theoretical biology_ 481 Elsevier, 2019, pp. 136–150
  * (41) Alec Sargood  “Gene Expression Time Delays in Reaction-Diffusion Systems” MSc Dissertation  In _arXiv preprint arXiv:2202.11544_ , 2022
  * (42) Alec Sargood  “GitHub Repository”  URL: https://github.com/alecsargood/Dissertation_Julia
  * (43) Razvan A Satnoianu, Michael Menzinger and Philip K Maini  “Turing instabilities in general systems”  In _Journal of mathematical biology_ 41.6 Springer, 2000, pp. 493–512
  * (44) J Schnakenberg  “Simple chemical reaction systems with limit cycle behaviour”  In _Journal of theoretical biology_ 81.3 Elsevier, 1979, pp. 389–400
  * (45) Natalie S Scholes et al.  “A comprehensive network atlas reveals that turing patterns are common but not robust”  In _Cell systems_ 9.3 Elsevier, 2019, pp. 243–257
  * (46) Stefanie Sick et al.  “WNT and DKK determine hair follicle spacing through a reaction-diffusion mechanism”  In _Science_ 314.5804 American Association for the Advancement of Science, 2006, pp. 1447–1450
  * (47) Christine N Tennyson, Henry J Klamut and Ronald G Worton  “The human dystrophin gene requires 16 hours to be transcribed and is cotranscriptionally spliced”  In _Nature genetics_ 9.2 Nature Publishing Group, 1995, pp. 184–190
  * (48) Alan M Turing  “The Chemical Basis of Morphogenesis”  In _Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences_ 237.641, 1952, pp. 37–72
  * (49) Lewis Wolpert  “Positional information and the spatial pattern of cellular differentiation”  In _Journal of theoretical biology_ 25.1 Elsevier, 1969, pp. 1–47
  * (50) Thomas E Woolley et al.  “Effects of intrinsic stochasticity on delayed reaction-diffusion patterning systems”  In _Physical Review E_ 85.5 APS, 2012, pp. 051914
  * (51) Thomas E Woolley, Andrew L Krause and Eamonn A Gaffney  “Bespoke Turing Systems”  In _Bulletin of Mathematical Biology_ 83.5 Springer, 2021, pp. 1–32
  * (52) Fengqi Yi, Eamonn A Gaffney and Sungrim Seirin-Lee  “The bifurcation analysis of Turing pattern formation induced by delay and diffusion in the Schnakenberg system”  In _Discrete & Continuous Dynamical Systems-B_ 22.2 American Institute of Mathematical Sciences, 2017, pp. 647