THE BIFURCATION ANALYSIS OF TURING PATTERN FORMATION INDUCED BY DELAY AND DIFFUSION IN THE SCHNAKENBERG SYSTEM

. A delayed reaction-diﬀusion Schnakenberg system with Neumann boundary conditions is considered in the context of long range biological self- organisation dynamics incorporating gene expression delays. We perform a detailed stability and Hopf bifurcation analysis and derive conditions for deter- mining the direction of bifurcation and the stability of the bifurcating periodic solution. The delay-diﬀusion driven instability of the unique spatially homoge- neous steady state solution and the diﬀusion-driven instability of the spatially homogeneous periodic solution are investigated, with limited simulations to support our theoretical analysis. These studies analytically demonstrate that the modelling of gene expression time delays in Turing systems can eliminate or disrupt the formation of a stationary heterogeneous pattern in the Schnakenberg system.

1. Introduction. A core question in developmental self-organisation concerns resolving the mechanisms by which long range spatial patterning emerges from an essentially homogeneous system. In 1952, Turing [42] suggested that interacting chemicals at a homogeneous steady state can be destabilized by diffusion, as explored via a system of two coupled reaction-diffusion equations. This class of instability is usually referred to as the Turing instability or the diffusively-driven instability. This concept has been confirmed in chemical systems [16,33] but awaits validation in biological systems, though numerous studies are indicative [5,18,19,20,22,23,24,34,32,39].
Over the years, Turing's idea has attracted the attention of a great number of investigators and has been explored theoretically. As well as being applied in biological pattern formation via morphogenic systems, Turing's idea has extended substantial influence over several fields including chemistry, economics, and semiconductor physics. Nevertheless, in the context of developmental biology, further complications arise when considering the fact that signaling morphogens induce each other's production via signaling pathways, which ultimately requires intracellular dynamics, including gene expression. This, in turn, requires DNA transcription and mRNA translation, which takes time. Empirical estimates for the magnitudes of transcriptional and translational timescales are estimated to be in the range of 10 minutes to several hours [21,41]. Recently, the impact of this delayed feedback due to the process of gene expression for the resulting morphogenic dynamics and its temporal behaviour has been intensively explored for Turing systems by modelling studies pursued mostly via computational simulations ( [15,36,37,38]), with the observation that that delayed dynamics can critically influence the prospect of robust patterns.
In particular, incorporating gene expression time delays into Turing models is most relevant for the self-enhancement production of a morphogen. Firstly, protein production by a cell requires gene expression, as observed in the context of morphogens by in-situ hybridisation which records mRNA levels, as illustrated in developmental self-organisation via Nodal and Lefty gene products in zebrafish mesendodermal induction [32]. Secondly, the extra-cellular regulation of such processes by an extracellular protein such as a morphogen requires signal transduction from the extra-cellular environment to the intra-cellular compartment and, for selfregulation, this signalling is also driven by the morphogen itself. Thus for Turing's mechanism in the context of cell biology, the extra-cellular morphogen binds a receptor on the cell membrane, consequently inducing a target gene expression [1], which results in the production of another new morphogen protein. Hence feedback delays appear in the production terms of the reaction kinetics, which are directly formulated via the law of kinetic mass action including time delays; in addition, the diffusion terms are unchanged since the delays are induced by intracellular dynamics, which is not coupled to the extracellular morphogen transport [15].
Furthermore, bound ligands, such as morphogens in the current context, have the potential to be endocytosed by cells and this mechanism is contemplated to be a means of regulation, by attenuating receptor signalling [40]. Thus, in addition to exploring reversible ligand binding (RLB) systems where the bound morphogen is released back into the extracellular environment, we also consider models where the bound ligand, here morphogen protein, is internalised via endocytosis. This results in an altered set of delayed kinetics, as detailed and derived by Seirin and Gaffney [36], and are referred to as ligand internalisation (LI) models below. In particular we will study both representations of morphogen signalling dynamics to assess the impact of delays on the behaviour of Turing systems, predominantly exploiting mathematical tools in distinct contrast to previous simulation based studies.
1.1. Main Equations and Objectives. Throughout this paper we consider Schnakenberg kinetics [27,35], firstly in the context of a reference RLB (reversible ligand binding) delayed model, whereby morphogen-regulated-morphogen-production proceeds via signal transduction driven by the reversible binding of ligand to the cell surface. The detailed model derivation can be found elsewhere [36] and, in summary, this model is given by Here x ∈ Ω ⊂ R is a one dimensional open bounded domain, t > 0 with the initial conditions and Neumann boundary conditions respectively given by In particular, u(x, t) and v(x, t) are concentrations of the activator and the inhibitor at (x, t) respectively, and a, b, d 1 , d 2 , τ are all positive constants. We also consider an alternative delayed model, based on ligand internalisation, and denoted as the LI model, which assumes that interactions proceed via morphogen internalisation to the cell interior, prior to the induction of the required signal transduction and downstream gene expression. As derived in [36], the LI model is given by together with the initial conditions and boundary conditions of equations (2), (3). The model systems (1) and (4) have been formulated based on the incorporation of morphogen-induced-morphogen production delays in the Schnakenberg system [27,35], which is recovered by setting τ = 0 in the system (1) or (4). Clearly, all these systems, (1)-(5), have a unique positive constant equilibrium solution E * = (u * , v * ), where Reaction-diffusion systems including time delays have been explored in many papers [6,7,8,9,13,14,17,25,28,30], but there has not been an investigation at the analytic level of how delays impact Schnakenberg Turing pattern formations. In particular, it is important to assess how delays influence the symmetry breaking dynamics from a homogenous steady state, as well as investigating the interaction of delay, kinetics and diffusion in self-organisation dynamics. Hence in studying the above delayed Schnakenberg systems, our first objective is to prove the existence (or absence) of spatially homogeneous and non-homogenous periodic solutions resulting from Hopf bifurcations, by using Hopf bifurcation theory, center manifold theory and normal form methods. With the help of analytical observations, we also confirm the presence of temporally periodic solutions of both spatially homogeneous and heterogeneous solutions via numerical simulations when relevant. Our final objective is to show that the interplay of delay and diffusion can induce the instability of the homogeneous solution leading to a stable stationary pattern and also to highlight cases where such an instability is not present.
1.2. Definitions and Notations. For clarity, we define the following instability properties. A solution of the system (1) or (4) is unstable via a diffusion driven instability if and only if the solution is stable in the system (1) or (4) when d 1 = d 2 = 0 but unstable for suitably chosen d 1 , d 2 > 0. A solution of the system (1) or (4) is unstable via a delay-diffusion driven instability if and only if the solution is stable when both d 1 = d 2 = 0, τ ≥ 0 and d 1 , d 2 > 0, τ = 0 but unstable for suitably chosen d 1 , d 2 > 0 and τ > 0. Trivially, the delay-diffusion driven instability is a special case of diffusion driven instability.
Throughout the paper, we denote by N the set of all the positive integers, and N 0 = N ∪ {0}. We respectively use FDEs and PFDEs to denote Functional Differential Equations (no diffusion but a delay is present) and Partial Functional Differential Equations (which possess both diffusion and delay terms).
The remaining parts of the paper are structured as follows. In Section 2, a detailed stability and Hopf bifurcation analysis for FDEs is presented. In section 3, we perform a stability and Hopf bifurcation analysis of PFDEs; in particular, we consider the delay-diffusion driven instability of the unique constant equilibrium solution and the diffusion-driven instability of the spatially homogeneous periodic solution of PFDEs. In section 4, the conclusions of the study are presented.

Dynamics of FDEs.
2.1. RLB model. To examine the dynamics of the homogeneous steady state E * = (u * , v * ) for the system (1), we first consider the following FDEs of the main equations (1): The linearized equations of the system (7) evaluated at (u * , v * ) are given by and thus the characteristic equation of system (2.1) is given by Since D(0, τ ) = u 2 * = 0, λ = 0 can never be a root of (8), so that the Bogdanov-Takens singularity will never arise in system (7) for all τ > 0.
Lemma 2.2. Suppose that τ j , with j ∈ N 0 , are defined as in (14). Then, Proof. Since the eigenvalue problem (9) is a special case of problem 4.1 of [2], we use Theorem 4.1 of [2] to prove this Lemma. Following [2], we define, Here in the context of our eigenvalue problem (9), we have, which implies that (θ(τ ) + 2jπ)/ω(τ ) = τ j , independent of τ . Then, by Theorem 4.1 in [2], we have, This completes the proof of the lemma.
Consequently, we obtain the following theorem directly from Lemma 2.1(2(b)) and Lemma 2.2. Theorem 2.3. Assume that (H 1 ) holds and let τ j be defined as in (14). Then, In what follows, we investigate the stability of the periodic solutions and consider the bifurcation direction of the local Hopf bifurcation associated with Theorem 2.3.
Let u(t) = u(τ t) − u * , and v(t) = v(τ t) − u * . Then, system (7) can be reduced to the following equations where we drop the tildes of u and v for the sake of simplicity if this is unambiguous. (15) is equivalent to the following abstract form in where , and, f By the Rieze Representation Theorem, there exists a 2×2 matrix function η(θ, τ ), with θ ∈ [−1, 0], whose elements are of bounded variation such that, for any φ ∈ Actually, we can choose where δ(θ) is the Heaviside function. We Taylor expand F (φ, τ ) as We next introduce a new parameter µ = τ − τ j , and then system (16) can be rewritten as where , we define the following bilinear inner product Let A be the infinitesimal generator of dU dt = L(τ j )U t , A and A * are adjoint operators, and ±iρ j (ρ j := τ j ω 0 ) are eigenvalues of A and A * . We let Λ = {−iρ j , iρ j }. Then, by applying the formal adjoint theory for functional differential equations in [11], C := C([−1, 0], R 2 ) can be decomposed by Λ as C = P ⊕ Q, where P is the generalized eigenspace associated with Λ, and its dimension is two. We extend the Let Φ and Ψ be the bases for P and P * associated with the eigenvalues ±iρ j of the adjoint equations, respectively. Note thatΦ = ΦB, where B = diag(iρ j , −iρ j ). We choose

FENGQI YI, EAMONN A. GAFFNEY AND SUNGRIM SEIRIN LEE
Clearly, (Ψ, Φ) = I, the 2 × 2 identity. As in [11], we enlarge the phase space C by introducing BC, denoted to be the space of the functions from [−1, 0] to C 2 which are uniformly continuous on [−1, 0) and which have a jump discontinuity at 0. The projection ϕ| is now replaced by π : BC → P , such that π(ϕ+X 0 α) = Φ(Ψ, ϕ)+Ψ(0)α. By using the decomposition U t = Φz(t) + y t , z ∈ C 2 , y t ∈ Q 1 (see [11,12]), we decompose system (20) as, with F 0 defined in (20), and X 0 = X 0 (θ) is given by In reducing to normal form, as in [12], we need to consider the following term: where f 1 j (z, y, µ), with j = 1, 2, are the homogeneous polynomials in (z, y, µ) of degree j with coefficients in C 2 and h.o.t denotes higher order terms. Then, in a finite dimensional locally invariant manifold tangent to the invariant subspace P of the linearized system (15) at z = 0, µ = 0, the normal form of system (15) is given byż where g 1 2 and g 1 3 are the second and the cubic terms in (z, µ). As in [11,12], we define, From (21), we have (19). A direct calculation shows that where Since by [12], we have Next, we need to consider the cubic terms g 1 3 (z, 0, µ). Since we concentrate on the properties of the local Hopf bifurcation, we only need to calculate the coefficient of z 2 1 z 2 instead of computing the coefficients of the terms o(|z|µ 2 ) which are irrelevant for the Hopf bifurcation.
From [3], it follows that, the sign of c 1 c 2 determines the direction of the bifurcation and the sign of c 2 determines the stability of the nontrivial periodic solution bifurcating from Hopf bifurcation on the center manifold.

Remark 1.
In order to confirm the existence of parameter space satisfying Theorem 2.4, we show the result numerically in Figure 1.

LI model.
To examine the dynamics of the homogeneous steady state E * = (u * , v * ) for the system (4), we now consider the following FDEs of the main equations (4): BIFURCATION ANALYSIS OF TURING PATTERN FORMATION 13 The linearized system of (37) at E * is simply given by and thus the characteristic equation is Theorem 2.5. Assume that (H 1 ) and Then E * is always asymptotically stable for all τ ≥ 0 so that there is never a Hopf bifurcation in system (37).
Then, we have the following results.
(1): If 0 <ε < 1 and then both the Turing instability and Turing bifurcation do not occur for all τ ≥ 0, whilst the Hopf bifurcation can only occur with τ = τ j , where the latter is precisely defined in (14).
Clearly, r k + q k > 0 for either ε ≥ 1, or Since now we know the sign of r k + q k , it will be sufficient to determine the sign of r k − q k in order to specify the sign of r 2 k − q 2 k . From (H 2 ), we have 2u * v * > u 2 * − 1 and thus Thus, given (46), we have r k − q k > 0 for k ≥ 1, so that (R 3 ) is satisfied. Note that the condition (46) is not dependent on k. Hence from Lemma 2.1(1), all roots of the characteristic equation (41) have negative real parts for all k (≥ 1) and all τ ≥ 0, so that both the Turing instability and Hopf bifurcation do not occur. However, when k = 0, according to Theorem 2.3, system (39) will undergo a Hopf bifurcation at τ = τ j near (u * , v * ).
From the proof of Theorem 3.1 (1), we know that the stability of E * is determined by the sign of the real eigenvalue of mode k = 0, which in fact corresponds to the FDE case. Thus we have only a Hopf bifurcation at E * for some large τ , and a homogeneous periodic solution under the given conditions. A simulation example is shown in Figure 2A. Further note that our simulation observations also suggest that when d 1 with ε < 1, we still have a homogeneous periodic solution though an asymmetric periodic solution appears initially as shown in Figure 2B.
In the following corollary, we consider the case of τ inf = τ 0 so that there does not exist a diffusion-driven Hopf bifurcation around E * . Corollary 1. (When τ inf = τ 0 , there is no diffusively driven Hopf bifurcation.) Suppose that u * = 2v * holds. Then, for ε sufficiently close to 1, we have τ inf = τ 0 so that there does not exist a Hopf bifurcation induced by diffusion around E * for all τ > 0.
Thus, we have τ 0 = τ inf . By the definition of τ k j in (44), τ k j is continuously dependent on the parameter ε. Thus, we can conclude that, τ inf = τ 0 still holds even for ε sufficiently close to 1.

Remark 2.
1. From the result of Theorem 3.1(1), there are no spatial patterns around E * for large d given (H 2 ) and ε < 1. This result also implies that if diffusion is sufficiently large (or if the spatial length is sufficiently small), Turing patterns never occur even for large delays. 2. If q k + r k < 0, that is R 2 fails, and (H 1 ) holds, then E * is always stable for the system (7)| τ =0 but unstable for any (d, ε) with the system (39)| τ =0 . Thus we have both possibility of a Turing instability and Hopf bifurcation for the system (39)| τ >0 within the parameter space q k + r k < 0. In fact, the numerical results of [37] show that the both cases actually occur. In the case that (39)| τ =0 has a Hopf bifurcation, we still have a Hopf bifurcation for E * in the system (39)| τ >0 in accordance with [4]. 3. If τ inf < τ 0 holds, then the instability of homogeneous steady state E * for τ ∈ (τ inf , τ 0 ) is induced by a joint influence of both diffusion and delay and we have a delay-diffusion driven instability. If either τ = 0 or d = 0, E * is always asymptotically stable. 4. In the case of τ = 0, E * is asymptotically stable for any d and ε in system (39) given the conditions (H 2 ) and either ε ≥ 1 or the condition (46). Thus the diffusion driven instability never occurs without delays under these conditions, and we cannot expect a Turing pattern, that is a non-constant stationary solution, induced by the destabilization of symmetric breaking by diffusion effects around E * . From Theorem 3.1, we also know that Turing patterns do not occur for all τ when ε < 1 given a sufficiently large diffusion coefficient. This implies that even though the two diffusion coefficients have different scales, pattern formation will fail on small spatial length scales. 5. Let us consider the case of 0 <ε < 1 and d < 1/(1 − ε)π 2 given the condition (H 2 ). In fact, this parameter space includes both regions of the Turing space and non-Turing space when τ = 0. When the diffusion rate is in the Turing space, we have observed the Turing instability failure with a large delay (see [37]). In contrast, we may have symmetric breaking triggered by the interplay of diffusion and delay, in the alternative case.

LI model.
In what follows, we show that the LI system, (4), has neither a Hopf bifurcation nor a Turing bifurcation for any τ ≥ 0. The rescaled LI system is given by

FENGQI YI, EAMONN A. GAFFNEY AND SUNGRIM SEIRIN LEE
The linearized system of (53) at E * is simply given by and thus the characteristic equation (41) is We have the following theorem which holds for all ε(≥ 0) and d(≥ 0).
Then E * is always asymptotically stable for all τ ≥ 0 so that there is never a Turing bifurcation nor a Hopf bifurcation for the PFDEs (53).
Here, the first term of (55) is from (H 1 ), the second term is from (H 3 ), and the final inequality, (56), is from (H 3 ). Thus the condition, (R 3 ), is satisfied. Note that all the above calculations hold for all ε(≥ 0) and d(≥ 0), or k ≥ 0. Hence, all roots of the characteristic equation, (41), have negative real parts for all τ ≥ 0 and k ≥ 0 from Lemma 2.1, so that E * is always asymptotically stable.
Remark 3.  [15], it has been shown that there is never a Hopf bifurcation at a point where the Turing instability occurs. Thus we conjecture that a Hopf bifurcation never occurs for all τ in LI PFDEs. 3. Note that the above results holds within the parameter space where the Turing instability does not occur in the absence of delay, with τ = 0 regardless of d and ε, because r k + p k > 0.
4. Concluding Remarks. By using bifurcation theory, normal form methods and center manifold theory, we performed a detailed stability and Hopf bifurcation analysis for the delayed reaction diffusion system associated with Schnakenberg kinetics, given the inclusion of gene expression delays for morphogen-induced-morphogen production. The results further allow an analytical basis for understanding how delayed biological self organization models behave and the parameter space for when bifurcations occur. In particular, for ligand internalization models where morphogens binding to a cell are internalised, the analysis demonstrates the absence of bifurcations, and thus no spatial organisation, emerges in the system dynamics. Hence, in generality, extensive ligand internalization appears to antagonize the induction of pattern formation in the context of gene-expression delayed Schnakenberg representations of the Turing instability.
In contrast, for reversible ligand binding models, the existence of spatially homogeneous and non-homogeneous periodic solutions was proved and temporally periodic solutions of both spatially homogeneous and heterogeneous solutions were confirmed by numerical simulations. In particular, a Turing instability can occur but for a sufficiently large delay, a Hopf instability can also be excited, leading to oscillations. Thus, spatial patterning without oscillations requires parameter fine tuning, in that the kinetic scales have to be tuned relative the delay size.
Furthermore, this analytical study provides analytical-level insight into the resulting Turing space of parameters; for instance it especially emphasises that the emergence of a Hopf bifurcation and thus oscillations with increasing delay is general for Schnakenberg kinetics with reversible ligand binding. Similarly, it demonstrates that extensive ligand internalisation antagonizes biological patterning in extensive generality for Schnakenberg models of the Turing instability. In particular, such observations are not restricted to a few parameter sets, a necessary limitation of simulation studies.
However, of course, the observation that the formation of heterogeneous steady patterning in the Schnakenberg model of the Turing instability is often inhibited by gene expression time delays does not necessarily entail that one should not expect to see a Turing instability in biological systems. Indeed, and in contrast, the empirical evidence is indicative of such an instability [5,18,19,20,22,23,24,34,32,39], even if molecular level details remain to be confirmed. Firstly, it may be that different kinetics are more robust to the inclusion of gene expression time-delays, motivating our future analytical studies of gene-expression delays for other models used to represent the Turing instability, such as the Gierer-Meinhardt system. More generally, the observations of this paper do not apply for the prospect of a generalized Turing instability resulting from more than two interacting morphogens plus a host of other prospective biological complexities. Indeed it is of curious interest that Turing did not immediately restrict his prospective mechanism to just two interacting chemical species in his original paper [42], whilst the observed tendency for two species models to fail to exhibit the Turing instability given gene expression time delays indicates that Turing's generality may yet prove important. Hence, the prospects of how gene expression time delays may impact on, or be compensated by, larger prospective biological self-organisation systems is still very much an uncharted area for theoretical study, both in terms of simulation and analysis, even though such systems are experimentally reported, at least to a limited extent, as illustrated with feather and hair patterning [10,26].