Instabilities in threshold-diffusion equations with delay

The introduction of delays into ordinary or partial differential equation models is well known to facilitate the production of rich dynamics, ranging from periodic solutions through to spatio-temporal chaos. In this paper, we consider a class of scalar partial differential equations with a delayed threshold nonlinearity which admits exact solutions for equilibria, periodic orbits and travelling waves. Importantly, we show how the spectra of periodic and travelling wave solutions can be determined in terms of the zeros of a complex analytic function. Using this as a computational tool to determine stability, we show that delays can have very different effects on threshold systems with negative as opposed to positive feedback. Direct numerical simulations are used to confirm our bifurcation analysis, and to probe some of the rich behaviour possible for mixed feedback.


Introduction
Delayed differential equations (DDEs) arise naturally as models of dynamical systems where memory effects are important. Indeed in models of population biology, ecology and epidemiology they are almost ubiquitous [1]. Delays also arise in many physiological systems as part of a feedback loop -one classic example being the pupil light reflex [2], though many others are wonderfully described in [3,4]. In contrast to models without delay the analysis of DDEs is notoriously hard and has generated considerable activity in the mathematics community (see for example [5]). This is directly attributable to the fact that the solution space for DDEs is infinite dimensional despite only a finite number of dynamical variables appearing in a model. Consider for example a simple scalar DDE that often arises in physiological modelling is a nonlinear function of v at time t − τ d . Even here the characteristic equation, 1/τ + λ − f (v) e −λτ d = 0, determining the stability of a fixed point v can have an infinite number of solutions making the spectral analysis challenging [6]. Probing the properties of fully nonlinear oscillations that are known to arise when f (v) is taken to be a "humped" function, such as the Hill function f (v) = θ n /(θ n +v n ), is harder still. However, progress is possible for this case in the limit n → ∞, so that f (v) becomes piecewise constant: f (v) = 1, v ∈ [−θ, θ] and is zero elsewhere. Information about the time where a solution crosses the threshold at θ is now enough to self-consistently determine a periodic oscillation in closed form [7]. Since threshold models are relatively common in the applied biological sciences, arising for example in models of calcium release [8], neural tissue [9] and gene networks [10], it is worthwhile developing a more comprehensive treatment of delayed threshold models. Since both positive and negative feedback are seen as important enhancers of the properties of biological systems [11] generic examples from both these classes, as well as a mixture of the two, should be considered. Moreover, in biological systems where transport is important it is common to use diffusion as a model for this process. With this in mind we concentrate in this paper on threshold-diffusion equations with delay.
In section 2 we introduce our model of choice, a delayed scalar PDE, and discuss some known results for smooth nonlinearities. We then introduce three distinct threshold nonlinearities describing negative, positive and mixed feedback models. As a precursor to the analysis of travelling waves we first neglect space and consider the generation of periodic oscillations in delayed threshold models in section 3. As well as briefly reviewing known results for the construction of such orbits we present the first treatment of their stability.
The Floquet exponents are given in terms of the zeros of a complex analytic function.
Importantly an examination of the spectrum shows that the model with negative feedback generates stable oscillations with period larger than twice the delay, whilst the positive feedback model generates unstable orbits with period less than twice the delay. In sections 4 and 5 we study travelling wave solutions of the full delayed PDE model, calculating wave speed and stability as a function of the delay. Once again we determine stability in terms of the zeros of a complex analytic function, which in the context of travelling waves we refer to as an Evans function. Computation of the Evans function shows that there is a stable front and an unstable standing bump in the model with positive feedback and that periodic waves can undergo instabilities in the model with negative feedback.
Direct numerical simulations are used to confirm our predictions, and also establish that such wave phenomenon are robust in the sense that they persist for smooth caricatures of the feedback nonlinearities. The case of mixed feedback is discussed in section 6. Here we show that, neglecting space, there can be more complex oscillations than occur in the case of purely negative feedback. We also demonstrate a spatio-temporal chaotic solution.
Finally we discuss natural extensions of the work in this paper.

The model
We consider a delayed scalar PDE for the variable v = v(x, t) with x ∈ R and t ∈ R + , which takes the form ∂v where τ d > 0 represents a fixed delay. Here f (v) is a given nonlinear function, τ > 0 is a decay constant and D > 0 a diffusion coefficient. We may interpret (1) as perhaps the simplest delayed reaction-diffusion equation imaginable. In particular, with f (v) = pv e −v we recover the diffusive Nicholson's blowflies equation [12]. Making use of comparison principles, monotone-iteration schemes and fixed-point theorems several results for monotone travelling waves of (1) are now known. For a recent discussion and an extension to cover nonmonotone travelling waves we refer the reader to the recent article by Trofimchuk and Trofimchuk [13]. However, the focus of the work in this paper will be on choices of f for which we may explicitly construct solutions of (1), as well as determine their stability. Motivated by the interest in various forms of threshold feedback in biological and physiological systems we consider three distinct choices for the nonlinear function f : Here, H is a Heaviside step function and h, h 1 and h 2 are constant thresholds. In all the cases I-III the presence of a threshold nonlinearity means that (1) can effectively be treated as a piece-wise linear (PWL) system. The study of PWL systems has allowed for important advances not only in the understanding of excitable systems [14] but also in the field of engineering [15,16].
To see how delay-induced dynamics can differ between the various types of feedback it is first worthwhile to treat the case of zero diffusion and assess dynamics in a point model.

Periodic solutions: zero diffusion limit
In the case of zero diffusion we recover the widely studied delayed ODE model Many results concerning the existence of periodic orbits are now known for the cases of smooth monotone positive and negative feedback [17]. Moreover, a Poincaré-Bendixson type theorem exists showing that chaotic behaviour is not possible in these cases [18].
However, for mixed feedback, as would occur for example in the Mackey-Glass model [3], very complex dynamics is possible [19]. Recent work of Röst and Wu [20] has established that this can include heteroclinic orbits from the trivial equilibrium to a periodic orbit oscillating around the positive equilibrium. However, we shall focus here on simple Tperiodic solutions of the form v(t) = q(t) with q(t) = q(t + T ). Importantly, with direct access to the shape of the periodic orbit, we shall show how to pursue a full Floquet analysis, without recourse to any approximations [21].

Negative feedback
For the case of negative feedback, f (v) = H(h − v), and with 0 < h < τ , there are no fixed points of (2), and it is natural to look for periodic solutions like those depicted in Fig. 1 (left). Assuming that q(t) crosses h from below, then after a delay τ d , q(t) evolves according to q(t) = A + e −t/τ , where time is measured from when q = A + . Next q(t) crosses h from above, and after a delay τ d evolves according to q(t) = A − e −t/τ + τ (1 − e −t/τ ), where time is now measured from when q = A − . Here A ± are the maximum and minimum values of q(t). Introducing an origin of time such that q(0) = A − and times t 1 and t 2 as in Fig. 1 (left), then q(t 1 ) = q(t 1 + t 2 + τ d ) = h and solving these equations for t 1 and t 2 and matching solutions at t = 0 and t = t 1 + τ d gives the period of oscillation T = 2τ d + t 1 + t 2

as [7]
where R = τ /h. The amplitude of oscillation is given by plot of the period and amplitude of the periodic orbit as a function of the delay is shown in Fig. 2 (left). We observe that both period and amplitude increase with the delay and that T > 2τ d . To establish the stability of the periodic solution we need to determine the Floquet exponents. Although for arbitrary smooth f this is a hard problem, for which progress has been made only in some special cases [22], we may make explicit progress here for the case of a threshold model.
Introducing the T -periodic function z(t) and the Floquet exponent λ we now look for solutions of the form u(t) = z(t) e λt u(0). Substitution into (4) shows that z(t) is the T -periodic solution oḟ Using the fact that the Green's function for (5) is temporal convolution. Since the periodic orbit q(t) has just two threshold crossings (at is T -periodic and takes the form Hence for + λ > 0 (i.e. to the right of the essential spectrum) the periodic solution can be calculated as where η P (t) is the periodic extension of η(t) restricted to the domain [0, T ). Substitution of t = t 1 and t = T − τ d into (7) Here Re λ < 0. To compute the zeros of E(λ) it is practical to first decompose λ as λ = ν + iω. The pair (ν, ω) may then be found by the simultaneous solution of E R (ν, ω) = 0 and Hence an examination of a plot of the zero contours of E I,R can be used to reveal the point spectrum -with Floquet exponents occurring where the two contours intersect. A plot obtained in this fashion is shown in Fig. 3 (left). From time-translation symmetry of the periodic orbit we see one exponent with λ = 0 as expected, and some others just slightly to the left of the imaginary axis. For the solution branch shown in Fig. 2 (left) no exponents are ever found in the right hand complex plane and so this branch of periodic orbits is stable.

Positive feedback
For the case of positive feedback, f (v) = H(v − h), and with 0 < h < τ there are two stable fixed points of (2) at v = 0 and v = τ . Assuming that a periodic orbit can coexist with these fixed points we look for a solution like that depicted in Fig. 1 (right). Here we have introduced the four unknowns A ± and T ± , which denote the largest (A + ) and smallest (A − ) values of the trajectory and the times spent above (T + ) and below (T − ) the threshold h. The trajectory increases from A − for a duration T + and decreases from A + for a duration T − . The values for these four unknowns are found by enforcing periodicity of the solution and requiring it to cross threshold twice as in Fig. 1 (right). The details of this calculation are presented in [23], where it is found that the period of oscillation where R = τ /h. The amplitude of oscillation is given by A plot of the period and amplitude as a function of τ d is shown in Fig. 2 (right). In contrast to the model with negative feedback we observe that τ d < T < 2τ d .
In a similar fashion to the calculation for negative feedback in section 3.1 we may calculate the Floquet exponents in terms of the zeros of a complex analytic function E(λ), which this time takes the form Here analysed in a similar way to the negative and positive feedback models described above.
For example orbits with v > h 1 that repeatedly cross h 2 will have properties analogous to those of the model with negative feedback, and so will be stable if τ > h 2 . We will not pursue this further here. Rather we move on to the question as to whether periodic orbits The existence of travelling front solutions in reaction-diffusion systems with delay is now reasonably well understood [24][25][26]. Recent work by Samaey and Sandstede has also shown how to determine the stability of waves [27]. These authors emphasise that the analytical determination of spectra is such a hard problem in general that resorting to numerical computation is sensible. However, for the special case of threshold nonlinearity we will show here that it is relatively easy to pursue questions relating to wave speed and stability.
In this section we shall first study travelling front solutions before moving on to standing bumps.

Travelling front
In a co-moving frame a stationary solution that connects the fixed point at v = 0 to the one at v = τ describes a travelling front (heteroclinic connection). Introducing ξ = x + ct, where v = v(ξ, t). A travelling wave solution q(ξ) is obtained upon letting ∂v/∂t → 0, so that q satisfies Now consider a monotone front solution where q(ξ) ≥ h for ξ ≥ 0 and q(ξ) < h for ξ < 0.
In this case where Continuity of q and q at ξ = cτ d gives the coefficients A and B as The speed of the front is then determined by demanding that q(0) = h, giving an implicit equation for c: subject to q(0) ≤ q(cτ d ). For τ d = 0 equation (16) may be re-arranged to give the speed explicitly as A plot of the front speed as a function of the delay τ d is shown in Fig. 4, showing that as the delay increase the speed decreases. Note that standing fronts with c = 0 occur when τ = 2h for any value of τ d . Our results are in contrast to those of, for example, the delayed Fisher equation [24], for which there is a range of speeds at which a front may travel. In our system the fixed points of the piecewise linear system (12), at q = 0 and q = τ , are both saddles, and hence we expect heteroclinic connections between them to occur only at isolated values of c.
To determine the stability of the front we consider perturbations of the form v(ξ, t) = q(ξ) + u(ξ, t). To first order we have from (11) that Perturbations of the form u(ξ, t) = u(ξ) e λt lead to the eigenvalue problem The Green's function, η, of the linear differential operator Q (Qη = δ) may be calculated (using Fourier transforms) as where Here we have assumed that λ + > 0 (i.e. we are to the right of the essential spectrum).
We may now solve (19) Demanding a non-trivial solution at ξ = 0 gives the condition E(λ) = 0, where Using the fact that k ± (0) = m ± we see that E(0) = 0, as expected for a system with translation invariance. The function E(λ) is real-valued if λ is real. Furthermore, the complex number λ is an eigenvalue of (19) if and only if E(λ) = 0. The function E(λ) thus provides us with an analytic tool for locating the point spectrum (isolated eigenvalues) of the linearised problem obtained after considering perturbations around a travelling wave solution. As such we shall refer to it as an Evans function (even though we do not attempt here to prove that the algebraic multiplicity of an eigenvalue is equal to the order of the zero of the Evans function). The Evans function has now been used to study the stability of travelling waves in a number of PDE models (see [28] for a recent discussion), though it is fair to say that for dissipative systems there are very few which possess explicit Evans functions [29]. For τ d = 0 then there is a single solution of E(λ) = 0 with λ = 0.
Indeed, for this case there are no solutions of E(λ) = 0 in the right half of the complex plane. Since E (0) > 0 the zero-eigenvalue is simple, and so the travelling front solution is stable. For τ d = 0, there are always non-zero solutions of E(λ) = 0 and a dynamic instability is in principle possible. It may be possible to show analytically that in practice these bifurcations do not occur, but here we just state that a numerical calculation of the spectrum shows that the front is stable even in the presence of delays. (Note here that all of our stability results relate to local stability; proving that a particular solution is globally attracting is more difficult [30].)

A standing bump
Here we consider time-independent standing waves v(x, t) = q(x) that satisfy For a symmetric bump solution that connects the fixed point at v = 0 to itself (homoclinic connection), crossing through the threshold h only twice, we may write with m ± = ± /D and q(−x) = q(x). Matching the solution and its first derivative at x = ∆, enforcing the threshold condition q(∆) = h, and fixing an origin such that For positive ∆ we require τ > 2h. The stability of the bump is determined along similar lines to section 4.1 using c = 0 and (Note that |q (−∆)| = |q (∆)|.) Perturbations to the bump in the form v(x, t) = q(x) + u(x, t) can then be shown to satisfy the equation where η(x) = e −k|x| /(2Dk), with k = ( + λ)/D. Substitution of x = ±∆ into (28) yields a pair of linear equations for the unknown amplitudes (u(∆), u(−∆)). Demanding a non-trivial solution of this system gives the spectrum as the zeros of the Evans function: It is simple to check that E(0) = 0 as expected from spatial translation symmetry. A further calculation also establishes that E (0) < 0 and that for λ ∈ R then lim λ→∞ E(λ) = 1. Hence, the Evans function must have at least one positive zero on the real axis, showing that the bump is always unstable.

Negative feedback: periodic travelling waves
The model with negative feedback and D = 0 was shown to support stable periodic oscillations in section 3.1. For non-zero diffusion it is thus natural to look for periodic travelling waves and see whether they are also stable. We shall show in this section that although a family of such waves can be supported, not all of them are stable. Our discussion treats separately two different types of domain: i) the unbounded real line; and ii) a finite system with periodic boundary conditions. In both these cases it is convenient to introduce the linear differential operator with m ± given by (14), so that the travelling wave ODE for q(ξ) may be written P q = H(h − q(ξ − cτ d )).

Infinite domain
For the case of an infinite domain the Green's function of P is given by where η ± (ξ) = exp(m ± ξ)/[D(m + − m − )]. We may now solve for q(ξ) using a convolution to give For a periodic wave oscillating around h we have that H(h − q(ξ)) = H(φ∆ − ξ)H(ξ) for some φ < 1 and q(ξ) periodic in ∆. We may then calculate q(ξ and Enforcing the two threshold conditions gives us two simultaneous equations for the two unknowns (φ, c), which can be solved numerically to obtain the dispersion relationships c = c(∆) and φ = φ(∆). Interestingly there can be many disconnected branches of solution to these dispersion curves. An illustrative plot of some of the dispersion relationships is shown in Fig. 5, together with results from direct numerical simulation. However, not all solutions of (35) correspond to valid periodic travelling waves. Each alternate family (shown with dashed lines in Fig. 5) actually corresponds to a wave for which q > h on (0, φ∆), whereas in the derivation of (33) and (34), we assumed that q < h on (0, φ∆). In Fig. 6 (left) we show a plot of a periodic wave obtained using the above prescription. Although it is possible to develop a stability condition along the lines for a front, this does not lead to any closed form expressions for the spectrum (as we get an infinite set of relations between perturbations at the threshold crossing points). However for a finite domain it is possible to make explicit progress as we show next.

Finite domain
Here we consider a finite domain of length ∆ such that the wave is periodic with q(ξ) = q(ξ + ∆). Moreover we shall focus on the case for which q(ξ) crosses threshold only twice as shown in Fig. 6 (right), though more general solutions can be constructed in a similar fashion. In this case q(ξ) is parametrised by the six unknowns (A 1 , A 2 , A 3 , A 4 , c, φ) as where φ is the fraction of the period for which q < h. Matching the solution and its first derivative at ξ = (1 − φ)∆ and ξ = ∆ can be used to determine the four amplitudes.
Enforcing the threshold conditions at ξ = (1 − φ)∆ − cτ d and ξ = ∆ − cτ d gives two further implicit equations for the pair (c, φ) that we may solve numerically to determine the wave properties. On doing this we recover the families of periodic travelling waves on an infinite domain found in section 5.1. For the calculation of stability we follow along similar lines to section 4.1 and make explicit progress using the result that the periodic Green's function of Q (given by the right hand side of (19)) can be obtained in closed form as with η P (ξ) = η P (ξ+∆) and k ± defined as in (21). The linearised equation for the evolution Using the result that we have where Here the derivative of the wave, q , is easily calculated from (36).
We in Fig. 3, and since they are isolated, zeros can be followed as parameters are varied.
In Fig. 7  Hopf bifurcation. Figure 8 shows this instability in a direct numerical simulation. We chose a domain of size 10 and an initial condition that led to a travelling wave with only two threshold crossings, as analysed above, so that initially ∆ = 10. With D = 1, this solution is known to be stable. In order to see the instability predicted by Fig. 7 we need to decrease ∆, i.e. decrease the domain size while keeping a periodic travelling wave with only two threshold crossings as the solution. However, from (1) it is clear that rescaling x : x → ax, is equivalent to rescaling D: D → D/a 2 . In the simulation shown in Fig. 8 we switched D from D = 1 to D = (10/7) 2 = 100/49 at t = 120, which is equivalent to switching ∆ from 10 to 7. As predicted by Fig. 8, the travelling wave becomes unstable and the resulting solution is a spatially homogeneous oscillation of the form described in In the final section we discuss solutions that emerge for mixed feedback and see that in some sense they can be regarded as hybrids of behaviour found for either positive or negative feedback alone. We conclude with some discussion about potential further work. known that (1) can support a travelling wave that is a hybrid of a monotone front and a periodic. This non-monotone front has an exponentially decaying profile in one direction that connects to a periodic oscillation in the other [13]. It is easy to replicate such solutions where F (v, a) = f (v) − a or F (v, a) = f (v − a(t − τ d )) for example, and τ a and θ are parameters. In models like this, a can be thought of as a negative feedback term. Such well as fronts) and having piece-wise constant nonlinear terms, can be analysed in the same way as described here. See [31,32] for a related analysis of systems with nonlocal interactions. Interestingly, delay induced periodic oscillations are also possible in nonsmooth models for a relay controller of an externally forced system with delayed feedback [33]. It would be interesting to assess the stability of such solutions with the techniques we have developed here.
All of the solutions we have discussed are structurally stable and are thus expected to persist for f smooth but sufficiently step-like. (For example, in the case of negative feedback one could take f (v) = 1 + e −β(h−v) −1 ; our results correspond to taking β → ∞.) We have numerically verified this for all solutions discussed (results not shown). For a smooth function f , the software package DDE-BIFTOOL [34] could be used to follow the front and standing bump solutions as heteroclinic and homoclinic connecting orbits, respectively [35]. The package can also follow periodic orbits, so could be used to study both the periodic orbits that occur in the zero diffusion limit and the periodic orbits corresponding to period travelling waves that we found. Another interesting possibility would involve the use of the function f (v) = H(v − h)e −r/(v−h) 2 (where r is a positive constant) in the positive feedback case, as used in [36]. This function is exactly zero for v < h, but is infinitely differentiable everywhere. Part of the travelling front described in Sec. 4.1 could be constructed in the same way: while the remainder of the front could be found by numerically shooting; solving the on cτ d < ξ (with a history given by (45)), and searching for the value of c for which the solution of (46) tends to the upper fixed point q * , which is a solution of q * = τ exp [−r/(q * − h) 2 ]. All of the above are topics of ongoing work and will be reported upon elsewhere.