Fixed and Distributed Gene Expression Time Delays in Reaction–Diffusion Systems

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 analyse 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, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau $$\end{document}τ, is increased. We show that the time to pattern formation from a perturbation of the homogeneous steady state scales linearly with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau $$\end{document}τ, 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.


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 (Bard 1992), vertebrate limb development (Miura and Shiota 2000a;Hentschel et al. 2004;Miura and Maini 2004), and pattern formation on animal coats (e.g. spots on a jaguar (Painter et al. 2000), feathers on birds (Bailleul et al. 2020)). Wolpert (1969) 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 (Gaffney and Monk 2006). 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 (Gaffney and Monk 2006), 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 (Gaffney and Monk 2006;Tennyson et al. 1995). 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 h (Kimmel et al. 1995). 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 1952). 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 (Yi et al. 2017;Harris et al. 2005;Miura and Shiota 2000b;Miura and Maini 2004;Sick et al. 2006). However, whether Turing patterns may be found experimentally in bio-logical systems with simple two-species systems is still very much an active field of research (Woolley et al. 2021).
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 1979), 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 (Seirin Lee and ). 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 Gaffney and Monk (2006) 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 et al. 2019;Yi et al. 2017) 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 Yi et al. (2017) indicate that extensive ligand internalisation may inhibit the formation of patterns within some regions of the parameter space. Results in both  suggest that time delay causes a significant effect on the time taken until pattern formation occurs. A final observation from Monk (2006), Seirin Lee et al. (2010), 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 Fadai et al. (2017Fadai et al. ( , 2018 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 Iron and Ward (2000). 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 (Maini et al. 2012;Scholes et al. 2019). 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 . 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 and Van Oudenaarden 2008;Elowitz et al. 2002;McAdams and Arkin 1997;Paulsson 2005). Distributed time delays have been incorporated to model biological phenomena such as hematopoiesis, and lactose operon dynamics (Cassidy 2021), Wnt/β-catenin signalling pathways (Cavallo et al. 2020), and Oncolytic virotherapy treatments for cancer (Elaiw and Al Agha 2020), among other applications. A distributed delay can be thought of as a more 'general and realistic' (Elaiw and Al Agha 2020) 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 et al. 2005;. 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 Sect. 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 Sect. 3. In Sect. 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 Sect. 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 Sect. 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.

Models of Gene Expression Time Delay
We consider two morphogens, u and v, obeying the following non-dimensional reaction-diffusion equations: where x ∈ = [0, 1], t > 0 and the function evaluations of u and v at some time delay s, which is an integration variable. We also have that τ 1 and τ 2 are the minimum and maximum time delays of the distributions, L 2 is the non-dimensional scaling of the domain length, and 2 the diffusion ratio between the activator u and inhibitor v. Unless otherwise stated, we use the same value of 2 = 0.001 as in Gaffney and Monk (2006), 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 Gaffney and Monk (2006), 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 et al. 1999). We primarily use no flux (homogeneous Neumann) boundary conditions on the boundary of the spatial domain, namely though we will consider an alternative set of boundary conditions in Sect. 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 , v ). Unless otherwise stated, initial conditions (u 0 , v 0 ) are chosen as a small random Gaussian perturbation from this homogeneous steady state given by where r u (x), r v (x) are normally distributed random variables for each x ∈ [0, 1] with zero mean and standard deviation σ I C . Unless otherwise mentioned, we take σ I C = 0.01 throughout. The stochastic nature of gene expression delays leads us to consider a mean-field approach to modelling the time delay (Bratsun et al. 2005;. 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; p), with s the integral variable s ∈ [τ 1 , τ 2 ], 0 < τ 1 < τ 2 and p the vector of distribution parameters. We will also consider non-Normal distributions by investigating skewed Gaussian distributions as well.

Symmetric Truncated Gaussian Formulation
Using a symmetric truncated Gaussian kernel, we have p = (τ, σ ), for some mean delay τ and standard deviation σ , with integration limits chosen as τ 1 = τ − nσ and τ 2 = τ + nσ for some n ∈ N, such that τ 1 = τ − nσ > 0 to ensure positive time delays only. The distribution is then given by where c denotes the truncation scaling constant, which ensures that K N (s; τ, σ ) integrates to 1 over the given integration domain [τ 1 , τ 2 ], and is computed as with φ(x) the cdf of the standard Gaussian distribution, 1 Throughout this paper, we use n = 3 so that the integration limits are τ 1 = τ − 3σ and τ 2 = τ + 3σ . This was chosen so that a relatively large σ value could be used for each τ while maintaining τ 1 > 0. For each τ value, a maximum σ value can be computed such that τ 1 = τ − 3σ ≥ 0 as σ max = τ 3 . By setting σ < σ max , we ensure that the integration domain strictly considers positive time delays only.
As σ → 0, we have that K N (s; τ, σ ) → δ(s − τ ). Using the sifting property of the delta function, we remark that the continuous time delay problem formulation in (1) reduces in the limit of σ → 0 to a discrete fixed delay case with delay τ , where Hence, the continuously distributed time delay model is a generalisation of the fixed time delay model. 1 The error function is given by erf(x) = 2 √ π x 0 e −z 2 dz.

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 p = (μ, ω, ρ). Such a distribution exhibits asymmetry but is otherwise a relatively simple generalisation of the Normal distribution. The probability density function of the skewed truncated Gaussian distribution is given by where φ(x) is the same as in (6). The parameter ρ denotes the skew factor. The distribution is negatively skewed for ρ < 0 and positively skewed for ρ > 0. Finally, we have that c is the truncation scaling constant, given as with (x, ρ) the cdf of a skewed Gaussian distribution, described by The function T (x, ρ) denotes the Owen's T function (Patefield and Tandy 2000) and is written as an integral in the form In the computational implementation of the skewed truncated Gaussian pdf, the integral T (x, ρ) is resolved numerically using the composite Simpson's rule with 100,000 discretisation points. The parameters μ and ω 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, τ , which is given by Following Flecher et al. (2010), the mean of the skewed truncated Gaussian distribution is computed as Fig. 1 Probability density functions of the (symmetric) normal, and skew normal distributions with a fixed mean of τ = 1 in all cases, and different standard deviations. The positively skewed distributions have ρ = 10, and the negatively skewed ρ = −10. Note that the axes labels cover very different ranges, with small values of σ and ω leading to very narrow distributions. For large values of ω, the negative and positive skew cases are not symmetric about the mean due to the truncation used. For this value of τ , we have that σ max ≈ 0.3333 and ω max ≈ 0.3354 withρ = 1 + ρ 2 1/2 . We note that, for a given ρ, all of the terms on the right-hand side of (12), namely ω, c , and K S (s; μ, ω, ρ), can be written explicitly in terms of μ. Equation (12) can therefore be solved implicitly for μ(τ ), for a given τ .
The integration limits were set to τ 1 = μ − 3ω, τ 2 = μ + 3ω, where ω was chosen such that ω < ω max , with ω max = μ 3 to ensure only positive time delays were considered. In Fig. 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.

Reaction Kinetics
We consider the Schnakenberg (1979) and Gierer and Meinhardt (1972) 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 (Woolley et al. 2012) is given by 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 Internalisation (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 (Seirin Lee and Gaffney 2010; Yi et al. 2017). 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 (Gaffney and Monk 2006). Considering mass-action kinetics and adding standard diffusion terms (Murray 2001), we nondimensionalise the resulting system to find that the LI model with a distributed time delay is given by The steady state of the Schnakenberg model is given by (u , v ) = a + b, b (a+b) 2 . We do not consider the RLB variant proposed in Seirin Lee and Gaffney (2010), as it was shown to be able to exhibit negative concentrations for positive values of the feed rates a and b (Dash 2020). We will also consider how this model compares to a fixed time delay LI model which is obtained in the limit σ → 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 Seirin . The two non-dimensionalised model descriptions we consider, with kinetic reactions taken from Murray (2001), and time-delayed terms motivated by Fadai et al. (2017) and Fadai et al. (2018), are given by (15) and (16), and labelled GM 1 and GM 2 , respectively, 2 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

General Linear Instability Analysis
We now consider the linear stability of homogeneous equilibria of the system (1). Denoting the steady state as (u , v ), we consider a small perturbation, u( Truncating at O(δ), we thus find that perturbations evolve according to where all partial derivatives are evaluated at the steady state, and (17) an ansatz of the form (ξ, η) T = e λ k t cos(kπ x)(ξ 0 , η 0 ) T yields a homogeneous linear system for (ξ 0 , η 0 ) T given by with Looking for non-trivial solutions of this system, we set det(M) = 0 to compute values of λ k . This leads to a characteristic equation (or dispersion relation) of the form, with the coefficients given as The coefficients α k , β k , γ k , δ k , and χ 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, GM 1 and GM 2 , in Table 1. Hence, one may note that the only term in the dispersion relation that varies with the distribution is E k (λ k ), which can be computed for a skewed Gaussian distribution via: (21) Setting μ = τ , ω = σ , and ρ = 0 to consider a symmetric distribution yields: Finally, taking σ → 0 results in as expected, corresponding to the fixed delay case with time delay τ . 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 λ 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 analyse these transcendental dispersion relations. We also remark that while (21) seemingly contains many parameters, all of these can be rewritten in terms of τ and σ , and hence compared directly with the symmetric case given in (22) (though the expressions for ρ = 0 provide no obvious insight, and again one must resort to numerical computation of λ k ).
The characteristic Eq. (19) can be used to determine the parameter sets (a, b, 2 , L 2 , τ, σ ) in which a Turing instability occurs (the 'Turing space') and hence where we expect pattern formation. If there exists a k = 0 for a given set of parameters such that max k ( (λ k )) > 0, then we expect pattern formation. The largest value of (λ 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.  Fig. 2 Turing spaces for the Schnakenberg and GM models without delay, with 2 = 0.001 and L 2 = 4.5. The dashed black arc (leftmost boundary) corresponds to parameters in which (λ 0 ) = 0, and the solid black outer arc (rightmost boundary) corresponds to parameters satisfying max k ( (λ k )) = 0 for k = 0

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.

Fixed Time Delay
In Fig. 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 (λ 0 ) = 0 for the spatially homogeneous characteristic equation. The outer red boundary is comprised of the points (a, b) such that max k ( (λ k )) = 0, k = 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 ( (λ 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 (Driscoll et al. 2014). In order to compute max k ( (λ k )), we varied k ∈ Z over [0, 50] for a given τ 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  (14), with 2 = 0.001, L 2 = 4.5. As τ increases, | max k ( (λ k ))| decreases. Stability lines for (λ 0 ) = 0 and max k ( (λ k )) = 0, k = 0, are overlaid, indicating the Turing space between them excited (and full numerical simulations were used to check this assumption throughout, as well as to check the predictions of pattern formation).

The LI Model
Considering first the LI model with fixed time delay, in Fig. 3, we plot a heatmap of max k ( (λ k )) over the (a, b) parameter space for the two fixed time delays τ = 0, τ = 1.5. We overlay contour lines corresponding to where (λ 0 ) = 0 and max k ( (λ k )) = 0, k = 0, highlighting the Turing instability region. As τ increases, the region of homogeneous instability (leftmost curve in Fig. 3 corresponding to (λ 0 ) = 0) decreases, increasing the Turing space. It can be seen that the absolute value | max k ( (λ k ))| also decreases for τ = 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 ( (λ k )) < 0, it will take a longer time for the eigenfunctions with modes k = 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 ( (λ k )) = 0 for k = 0, does not move at the resolution of the plotting τ changes, while analogous results (not shown) were found for other values of τ and (Sargood 2022a).
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 Sanchez-Garduno et al.

The GM Models
The results in Fadai et al. (2017) showed that an increasing time delay in GM 1 had an antagonistic effect, shrinking the parameter space exhibiting stable spike solutions. In Fig. 4 Maximum growth rate, max k ( (λ k )), corresponding to solutions of (19) for models (15) in (b) and (16) in (c) with a fixed time delay, plotted for a, b ∈ [0, 1] × [0, 4], with varying τ . Note that these models are equivalent for τ = 0 in a. Parameters 2 = 0.001 and L 2 = 4.5 used. Stability lines for (λ 0 ) = 0 and max k ( (λ k )) = 0, k = 0, are overlaid, indicating the Turing space between them contrast, an increasing time delay for GM 2 caused an expansion of the stable spike solution parameter regime (Fadai et al. 2018). 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 Fig. 4b show a shrinking Turing space for increasing τ for GM 1 , consistent with the analysis conducted in Fadai et al. (2017), which showed a de-stabilisation of the stable spike solution parameter space with an increasing τ . Similarly, results in Fadai et al. (2018) showed a stabilising effect of increasing τ on the stable spike solution parameter space for model GM 2 in (16), and we find an analogous result for the Turing space, shown in Fig. 4c. In both cases, we also observe a decrease in | max k ( (λ k ))| as the delay is increased. We remark that the choice of delay times τ for these, and other plots using the dispersion relation (19) was sufficiently large to demonstrate key trends for increasing τ , 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).
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 timedelayed 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 τ 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.

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 σ corresponds to a measure of the width of this distribution, with the limit σ → 0 capturing the fixed time delay case. In particular, σ max has been defined to be the largest feasible value of σ to avoid negative delays, as described in Sect. 2.1. We consider the absolute difference of max k ( (λ k )) for varying σ values as a fraction of σ max , for multiple τ and values. For each (τ, 2 ), bifurcation plots analogous to those of the previous subsection have been computed for the distributed delay case with varying σ ∈ {σ max × 0.99, σ max × 0.2, σ max × 0.1}. However, the differences are below plotting resolution when compared to the fixed delay case, as may be found in Sargood (2022a); hence for each (τ, 2 ), we instead tabulate the absolute difference of max k ( (λ 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 ( (λ k )) in Table 2 for all σ , τ and 2 considered across the parameter space (a, b) ∈ [0, 1.4] × [0, 2] is O(10 −3 ). We therefore expect that for all (a, b) ∈ [0, 1.4] × [0, 2], using a symmetric Gaussian distribution centred at some mean τ (for small τ ) will not significantly affect Turing instabilities compared to the fixed delay case, independent of the standard deviation σ 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.

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 τ 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 Fig. 2a, and compare these plots to the fixed delay case. We show that for a small mean τ , the skew, positive or negative, does not significantly affect the value of max k ( (λ k )). We highlight here again that Eq. (12) can be solved implicitly for μ(τ ), for a given τ , using the fzero command in MATLAB. For a given ρ, and each found μ, we compute max k ( (λ k )) by solving for roots of the characteristic Eq. (19). In Fig. 5, we plot max k ( (λ k )) against τ ∈ [0, 0.8] for skew parameter values of ρ = −10, 10, and ω = ω max × 0.99, with two different (a, b) parameter sets. The value ω max is defined analogously to that of σ max . A plot of max k ( (λ k )) for the fixed delay case is also added for comparison in each case.  (19) plotted against τ ∈ [0, 0.8] for ρ = −10, 10 against the fixed delay case, with the distribution spread given by ω = ω max × 0.99. Parameter values 2 = 0.001 and L 2 = 4.5 used. In generating the plot, the time delay scale, τ , was varied at regular intervals of 0.05 with k ∈ Z ranging over k ∈ [0, 50] before a maximum was taken From Fig. 5, we see that the curves differ slightly for small τ , with ρ = −10 having a slightly higher value of the maximum growth rate, and ρ = 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 ( (λ k ))| in Figs. 3, 4b, 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.

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 discretised via the method-of-lines using the standard threepoint stencil for the Laplacian, with m = 500 equally spaced points on the domain x ∈ = [0, 1]. This discretisation 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 (Rackauckas and Nie 2017;Rosenbrock 1963). The distributed delay terms were approximated via a composite Simpson's rule using 50 quadrature points; see Sargood (2022a) 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 (Sargood 2022b).
With a maximum delay of τ 2 = τ + 3σ (or μ + 3ω with a skewed-distribution), we note that to solve DDEs, a history function is required to define the solution for t ∈ [−τ, 0) for the fixed delay case, t ∈ [−τ − 3σ, 0) for the symmetric distributed case, and t ∈ [−μ − 3ω, 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 2001). The concentration gradients of the two morphogens u and v are thus effectively 'out

A Relationship Between Fixed Delay and Onset of Patterning
Here we show that for small τ and L 2 , the linear theory provides a good approximation to the time taken until pattern formation occurs, and in fact, the relationship between τ 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 τ and time-topattern on a longer timescale for larger τ is also linear. We first consider the former.
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 Fig. 6, where (λ k ) is plotted against k for two different domain sizes, for a given (a, b, τ ). We also use an initial small perturbation of the form of Eq. (3) for t ∈ [−τ, 0], but vary the standard deviation, σ I C , of these perturbations.
The linear theory suggests that the perturbation will be of the form ξ(x, t) ∼ A k (t) cos(kπ 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, 2 , τ, L 2 ), we can solve the characteristic Eq. (19), and plot (λ k ) against k, to determine the dominant mode k and the corresponding growth rate, λ k . When the perturbation ξ has grown sufficiently, in absolute value, beyond a threshold where pattern formation is considered, we call this time t = T . More specifically, the time T is the first time such that max x |u(T , x)− u | > u T , where u T 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 as the 'true' time-to-pattern. Finally, using the relation A k (T ) ≈ A k (0)e λ k T , we can rearrange for T and thus compute a linear theory approximation for the 'true' time-to-pattern, via The approximate time-to-pattern, T , can thus be determined via λ k , from (19), as well as A k (0) and A k (T ), computed via the Fast Fourier Transform. We repeat the evaluation of both the true and predicted times to pattern for varying (a, b, τ ), and compare the evaluations in Fig. 7, which shows the comparison for three different parameter sets for (a, b), with the time delay varied over τ ∈ [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 λ k for the dominant mode accurately predicts the time to pattern, which is observed to vary linearly in τ to a very good approximation (Sargood 2022a). We remark that computing λ k for τ > 1.6 is numerically difficult; thus full numerical solutions were used to plot 'true' time-to-pattern against τ 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 Fig. 8, we consider the time delay τ ∈ [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 σ IC = 10 −5 , which once more reveals a linearly increasing relationship between T and τ .
For the GM 1 model (15), as a result of the shrinking Turing space, we only consider small τ ∈ [0.1, 1], varied at regular intervals of 0.1. For the GM 2 model (16), we consider both small τ ∈ [0.1, 2] and larger τ ∈ [1, 16], varied at regular intervals of 0.1 and 1, respectively. We set an initial perturbation from the steady state as σ 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 Fig. 9, for both Gierer-Meinhardt models, Eqs. (15) and (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-ofmagnitude increase in the time-to-pattern.

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 ( (λ 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 Figs. 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 τ ∈ {1, 16} and varying σ , and compared with the appropriate fixed delay case given an initial small perturbation of the form of Eq. (3) for t ∈ [−τ, 0]. The results illustrate that the onset of patterning, and the emergent pattern do not depend on the value of σ used at the resolution of the plotting. We next present numerical results to show that even the skew of an asymmetric distribution for both a small and large mean delay τ , does not have a significant effect on observable patterning results, as highlighted for the LI model by Fig. 12 for τ = 0.1 and Fig. 13 for τ = 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 Sargood (2022a).

Boundary and Initial Conditions
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 we proceed to consider three different sets of initial conditions. We denote by IC 1 initial conditions based on those used in Gaffney and Monk (2006). With u 1 (x) and v 1 (x), the functions depicted in Fig. 14, these conditions are given by for t ∈ [−τ, 0] and represent a specific perturbation of the homogeneous steady state. We denote by IC 2 the initial conditions defined in Eq. (3) with σ I C = 0.01, and IC 3 denote the same initial conditions, but with σ I C = 0.1. As previously, both initial conditions were specified for t ∈ [−τ, 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 Eq. (3) are the same in each simulation. Figures 15 and 16 show simulations for each of these initial conditions for varying fixed time delay τ ∈ {1, 16}. The final pattern is somewhat sensitive to the choice of initial conditions. As one might expect, the larger σ IC used in IC 3 , compared to that of 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 , v 1 (x) used to generate the initial conditions IC 1 ; algebraic expressions for these functions are given in Appendices of Gaffney and Monk (2006) 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 Sargood (2022a) for examples.
Finally, we consider the effect of varying boundary conditions. Motivated by the analysis in , homogeneous Dirichlet boundary conditions are implemented for the activator term, and homogeneous Neumann boundary conditions for the inhibitor term. Hence, we set The results in Fig. 17 were generated using IC 2 , with time delays of τ = 1 and τ = 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  as 'isolated' patterns, we do not observe a noticeable difference in the time taken for patterns to form.

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 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 Gaffney and Monk (2006). 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 Fadai et al. (2017Fadai et al. ( , 2018, 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 and Monk (2006) that a Hopf bifurcation cannot occur for a spatially inhomogeneous mode with nonzero wavenumber for this system. Thus, the growth rate λ k is real at a bifurcation for k = 0, and hence at the boundary of the Turing space dictated by the spatially inhomogeneous modes, we have λ k = 0. Given the time delay only occurs in the dispersion relation, Eq. (19), via terms of the form exp(−λ k τ ), these boundaries thus must be independent of τ. 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 behaviour 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 reactiondiffusion equations which can have approximately the same linearisation as delayed reaction-diffusion systems for small delay (Ritchie et al. 2022).
Finally, driven by the inherent stochasticity of the molecular processes underpinning gene expression (Raj and Van Oudenaarden 2008;Elowitz et al. 2002;McAdams and Arkin 1997;Paulsson 2005), 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 (Gaffney and Monk 2006). 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 Satnoianu et al. 2000;Scholes et al. 2019). 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 and Scott (1986) or Thomas (Murray 2001) 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 (Gaffney and Monk 2006;Sanchez-Garduno et al. 2019). 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 Krause et al. (2020) which have similar numerical difficulties due to the presence of exponential functions.