DYNAMICS OF A DELAY DIFFERENTIAL EQUATION WITH MULTIPLE STATE-DEPENDENT DELAYS

We study the dynamics of a linear scalar delay differential equation εu̇(t) = −γu(t) − N ∑ i=1 κiu(t− ai − ciu(t)), which has trivial dynamics with fixed delays (ci = 0). We show that if the delays are allowed to be linearly state-dependent (ci 6= 0) then very complex dynamics can arise, when there are two or more delays. We present a numerical study of the bifurcation structures that arise in the dynamics, in the nonsingularly perturbed case, ε = 1. We concentrate on the case N = 2 and c1 = c2 = c and show the existence of bistability of periodic orbits, stable invariant tori, isola of periodic orbits arising as locked orbits on the torus, and period doubling bifurcations.


1.
Introduction. Delays are ubiquitous in biology, arising from maturation, transcription, incubation and nerve impulse transmission time, to name but a few situations. Classical examples of delay equations in mathematical biology include the Mackey-Glass equation [24], Nicholson's blowflies equation [15], and the delayed logistic equation, also known as Wright's equation, after a change of variables [19,38,23,36]. Such equations have inspired decades of mathematical research into delay differential equation (DDEs), and there is now a well established mathematical framework for problems with fixed or prescribed delays as infinite-dimensional dynamical systems on function spaces (see [17,6], or [36] for an easier treatment). However, some biological delays such as maturation and transcription delays would seem to be more naturally modelled as state-dependent delays, see for example [33] for evidence of state-dependency in the maturation time of neutrophil precursors. On a larger scale, the maturation age of juvenile seals and whales has been observed to depend on the abundance of krill [12]. Mathematical models with state-dependent delays appear in many contexts, for example in milling [21], control theory [37], economics [25] and population dynamics [2].
There has been much work in recent years to extend the general theory of DDEs to allow for state-dependent delays (see [18] for a review). However, many papers concentrate on problems in a singular limit, or near equilibrium or with a single state-dependent delay. Much less theory has been established for equations with multiple state-dependent delays. Mallet-Paret et al. [28] showed the existence of slowly oscillating periodic solutions of (1) when a i = a for all i, Györi and Hartung [16] considered the stability of equilibrium solutions, and Eichmann [8] proved a local Hopf bifurcation theorem for multiple state-dependent DDEs. Rigorous theorems for other bifurcations of periodic orbits in multiple state-dependent DDEs have yet to be proven, although Sieber [35] suggests an approach for doing so, as well as providing an alternative proof of the Hopf bifurcation theorem.
It is difficult to envisage physiological models with multiple state-dependent delays gaining much traction while the theory of these equations is so incomplete, but at the same time it is difficult to develop theory without models and examples to work from. Before breaking this impasse it is desirable to understand the dynamics that result from simple models with multiple state-dependent delays. Accordingly, we study the dynamics of the model scalar multiple state-dependent DDE where the coefficients ε, γ, N , c i , κ i and a i are all strictly positive. After developing some theory for the general case we present numerical computations of the bifurcations and invariant objects in the case N = 2 and c 1 = c 2 = c. We show that a wide range of dynamical behaviour is exhibited, including stable and bistable periodic solutions, period-doubled solutions, stable tori with quasi-periodic solutions and with phase-locked periodic orbits, together with the associated bifurcation structures. We choose the problem (1) because the state-dependent delays are essential for driving all of the dynamics seen, since there is no nonlinearity in this system apart from the state-dependency of the delays. Indeed setting c i = 0 in equation (1) results in a linear constant delay DDE with no interesting nontrivial dynamics. For c i > 0 the delays are merely linearly state-dependent, but having two or more such delays is sufficient to create the cornucopia of dynamics that we observe. Thus our main result is to report that the presence of multiple state-dependent delays is sufficient to generate very complex dynamics, even in the absence of any other nonlinearity in the model. Consequently, neglecting the state-dependency of delays in mathematical models has the potential to dramatically alter the dynamics.
Equations of the form (1) and (2) have been studied in a series of papers by Mallet-Paret and Nussbaum [26,27,28,29,30,31] mainly concentrating on establishing the existence, stability and shape of slowly oscillating periodic solutions in the singularly perturbed case 1 ≫ ε > 0 with a single state-dependent delay. In particular, the linearly state-dependent delay equation (1) with N = 1 is studied in [31]. The existence of a stable slowly oscillating periodic solution of (2) is established in [28] for an arbitrary number of delays under the condition that α i (t, 0) = t − a for all i (so all the delays are equal when u = 0). Equation (1) with a i = a for all i is considered as an example. Numerical computations of some initial value problems of the related equation (4) with N = 2 have been presented by John Mallet-Paret in seminars, but have not been published. Those computations inspired our more systematic study of the bifurcation structures of these equations. In Section 2 we establish existence and uniqueness of solutions of (1) as an initial value problem, under the bound γ > N i=2 κ i . In particular we show that none of the α i (t, u(t)) become advanced, and establish a bound on solutions. We also show that under the condition c i = c for all i, the α i (t, u(t)) are monotonically increasing functions of t, possibly after some initial transient.
In Section 3 we consider the linear stability of the equilibrium solution of (1) and show that the only possible local bifurcations from the equilibrium solution are Hopf bifurcations. We also establish a necessary condition for instability, and identify the asymptotic distribution of the Hopf bifurcation points. We then find bounds on the amplitude of any oscillatory solutions of (1) using a Gronwall argument.
In Section 4 we present a numerical study of the dynamics and bifurcations of (1) with κ 1 as a bifurcation parameter. In the case of one delay, the equilibrium solution loses stability in a supercritical Hopf bifurcation at κ 1 = κ * 1 > γ, and a branch of stable slowly oscillating periodic solutions (SOPS) is created which persists for all κ 1 > κ * 1 . There are infinitely many other Hopf bifurcations, but each results in a short period low amplitude unstable periodic orbit. There are no other bifurcations. In the rest of the paper we consider the dynamics when there are two delays in the non-singularly perturbed case, ε = 1. The bifurcation diagram in the two delay case is similar to that of the one delay case when the coefficient κ 2 of the second delay is small, and even if it has moderate values when the parameters are chosen so that the period T 0 at the first Hopf bifurcation satisfies T 0 > a 2 − a 1 .
In Section 4.2 we show that the branch of stable periodic solutions undergoes two saddle-node bifurcations resulting in bistability of periodic solutions for a range of values of κ 1 when κ 2 = 1 and a 2 is sufficiently large so that T 0 < a 2 − a 1 . This bistability region is seen in later sections too for all larger values of κ 2 , and so occurs for a region in the (κ 1 , κ 2 ) plane. Bistability is important in biology, where it is the dynamic origin of some biological switches. Two of the most studied and understood natural examples are the lactose operon and the phage-λ switch, both of which can be modelled with (constant delay) DDEs [34,39]. Bistability of periodic orbits has been observed in DDE models of dynamical diseases [11]. It is also possible to bioengineer bistability [10] by constructing systems with positive feedback loops, or double negative feedback loops. The model (1) gives another mechanism for bistability through the interaction of two linearly state-dependent delays in an otherwise linear system. Bistability is also important in the context of Wright's conjecture where, if found, it would disprove the conjecture [23].
In Section 4.3 we find regions of parameter space where the principal branch of periodic solutions loses stability through a torus bifurcation, and there is no stable periodic orbit or equilibrium solution. We are able to compute the resulting stable torus and follow the branch of torus solutions in parameter space between torus bifurcations on the principal branch of periodic solutions where the torus is born and dies. We verify the torus structure by computing a Poincaré section. In contrast to the ODE case, because the phase space of the DDE is an infinite-dimensional function space, the Poincaré section is itself an infinite-dimensional function space, and we plot a projection of the Poincaré section into R 2 to reveal the torus structure. Green et al. [13] studied tori and associated Poincaré sections and their projections (as well as bistability of periodic orbits) in the context of a fixed delay DDE arising in a model of a semiconductor laser with feedback.
In Section 4.4 we present other bifurcations and dynamics seen in this system, including examples of stable periodic orbits not on the principal bifurcation branch. In particular, there is a period-doubling bifurcation which results in a stable branch of period-doubled orbits. This gives an example of a stable periodic orbit on a secondary bifurcation branch in this system. We also find isolas of periodic orbits on the bifurcation diagram through phase locking on the torus. This gives another example of a stable periodic orbit not on a primary bifurcation branch, but this time isolated from the primary, secondary and subsequent branches. Given that we find tori, we also investigate and find the existence of double Hopf bifurcations. We conclude by showing an example where a small parameter perturbation changes the connectivity of the branches between the Hopf bifurcations. This change indicates that there are other bifurcations yet to be identified in the dynamics.
2. Existence of Solutions. Although we use the expression, "state-dependent delay differential equation" is dangerously misleading. In contrast to the case of equations with fixed or prescribed delays, if the offset arguments α i in equation (2) depend on the state u(t) of the system at time t then α i being a delay is a property of some solution trajectory under consideration, and not a property of the differential equation itself.
We consider (1) and without loss of generality order the arguments α i so that Since α i (t, u(t)) ≡ t − a i − c i u(t) < t provided u(t) > −a i /c i , we have that all the arguments α i (t, u(t)) are retarded when u(t) > −a 1 /c 1 . Although equations with advanced and retarded arguments are interesting in a number of settings, including Wheeler-Feynman electrodynamics [4,5] and travelling waves for lattice differential equations [1], in the current work we will restrict attention to the case where all the α i (t, u(t)) are retarded in equation (1). In Figure 1 we present numerically computed solutions for with N = 2 and −2 = −a 1 /c 1 > −a 2 /c 2 = −5. The definition of α i (t, u(t)) in (4) ensures that α i (t, u(t)) t so the computation does not fail because an argument becomes advanced. Numerically computed solutions for different initial conditions are shown in Figure 1 and appear to be converging to the same periodic solution (but with different phases). If these solutions satisfied u(t) −a 1 /c 1 they would also be solutions of (1). However, the solutions found have u(t) < −a 1 /c 1 on parts of the orbit, and so are not valid as solutions of (1). To guarantee that solutions to (1) do not terminate because an argument becomes advanced, we need to restrict the range of parameter values under consideration.
Since (1) and (4) are equivalent for u(t) ∈ (L 0 , M 0 ) the following result is immediate. In the case of one delay, N = 1, Corollary 2.2 bounds solutions in the interval (−a 1 /c 1 , κ 1 a 1 /γc 1 ) with no condition on the parameters, a result previously found in [31]. In the case of two delays the condition (5) becomes γ > κ 2 . For general N > 1 we can fix the parameters κ 2 ,. . . ,κ N and γ such that (5) holds and Corollary 2.2 guarantees the existence of bounded solutions for any κ 1 > 0. This suggests that κ 1 is a suitable bifurcation parameter for studying the dynamics of (1).
The proof of Corollary 2.2 does not require that the delays t − α i (t, u(t)) are bounded away from zero (or equivalently using (3) that u(t) is bounded away from −a 1 /c 1 ), but this can be easily established. To do so, suppose u(t) < 0 andu(t) 0 and hence u(α i (t 2 , u(t 2 ))) < 0 for all i. Then (1) implies εu(t 2 ) > 0, a contradiction, and hence t 2 − t 1 < max i {a i }. Thus each interval on which u andu are both negative is bounded with u(t) having a local minimum at the end of the interval, which by Theorem 2.1 is bounded away from −a 1 /c 1 . In Section 3 we will establish explicit bounds on solutions u(t) using a Gronwall type argument.
By Corollary 2.2, u(t) is continuously differentiable for t > 0, but in general we expect that lim tր0u (t) = lim tց0u (t), that is the derivatives of the given history function that defines the IVP, and the solution itself do not agree at t = 0. In this case t = 0 is called a 0-level primary discontinuity point and following standard notation (see [3]) we label it ξ 0,1 , where the first subscript indicates the level and the second is an index. At the points t = ξ 1,j where α i (ξ 1,j , u(ξ 1,j )) = 0 for some i, the discontinuity inu(t) at t = 0 will cause a discontinuity inü(t) at t = ξ 1,j , and these points are called 1-level primary discontinuity points. Similarly at points t = ξ 2,k such that α i (ξ 2,k , u(ξ 2,k )) = ξ 1,j for some i and j there can be a discontinuity in u (3) (t). In this way the discontinuity inu(t) at t = 0 propagates to discontinuities in higher derivatives of u(t) at later times. However, the bounds u(t) ∈ (L 0 , M 0 ) for solutions of (1) define upper and lower bounds on the delays α i (t, u(t)) via the linearity of the state-dependency of the latter and so we have t − α i (t, u(t)) ∈ (a i + c i L 0 , a i + c i M 0 ), which implies that α i (t, u(t)) → +∞ as t → ∞. These bounds also allow us to bound the discontinuity points (also called breaking points), since t − a i − c i M 0 α i (t, u(t)), we have ξ n+1,k ξ n,j + τ when α i (ξ n+1,k , u(ξ n+1,k )) = ξ n,j , and hence the n-level primary discontinuities satisfy ξ n,k nτ for all k. Thus u(t) ∈ C n+1 for all t nτ . Discontinuities in the derivative of the history function give rise to so-called secondary discontinuity points, which propagate in a similar manner, but since the 0-level secondary discontinuity points satisfy ξ 0,i < 0, we still have that u(t) ∈ C n+1 for all t nτ . See [3] for a full discussion of breaking points in DDEs.
Regarding (1) as a dynamical system, the phase space consists of function segments containing the necessary solution history to integrate for all future time. For fixed delays, c i = 0, equation (3) implies that the largest delay is a N so at any The situation is more complicated for state-dependent problems since the amount of history that needs to be retained will vary with the delay. If the α i are monotonic (increasing) functions of t then at , and in this case the amount of history required to integrate the solution depends on the as yet undetermined solution itself. Fortunately, for (1) under (5) and (3) Hence on any invariant set (and in particular on any periodic orbit) or on an arbitrary orbit for t > T solutions are contained in the space where we note that the length of the time interval is defined by the value of the function at the right-hand end point.

Linear Stability and Hopf
Bifurcations. Linearization of state-dependent DDEs was for a long time done heuristically by freezing the value of the statedependent delay at its equilibrium value. Justification of this procedure is more recent [18]. Linearizing (1) about the trivial equlibrium solution u(t) = 0 we obtain with solutions of the form u(t) = j β j e λj t where the λ j are solutions of the characteristic equation Since, with positive coefficients, g(λ) > 0 for all λ 0 we conclude that any real characteristic values are negative, and so there are no local bifurcations from u(t) = 0 where a real eigenvalue crosses zero. (Since g ′′ (λ) > 0 for all real λ, we also note that there are at most two real negative characteristic values.) Considering complex characteristic values, λ = x + iy, take real and imaginary parts of (12), square and add to find that Thus x 0 the right-hand side of (13) is less than or equal to ( N i=1 κ i ) 2 while left-hand side is strictly greater than γ 2 , thus the characteristic values lie on a curve completely contained in the left half-plane. In [16] it is shown that the equilibrium solution of (11) is exponentially stable if and only if the equilibrium solution of (1) is exponentially stable. Hence is a necessary condition for both the trivial solutions of (1) and (11) to be unstable. This condition is not sufficient for instability, since even when the curve crosses into the right half-plane, the characteristic values may still lie in the left halfplane. To determine when the equlibrium solution is unstable, we need to find when complex conjugate pairs of characteristic values cross the imaginary axis. These crossing result in Hopf bifurcation in the full state-dependent problem (1), as follows from the recently elaborated theory of Hopf bifurcations for state-dependent DDEs [8,20,35].
In the case of a single delay, N = 1, the Hopf bifurcations of (1) are well known to lie on the paramaterized curves (15) where for fixed γ > 0, as κ 1 is increased there is a infinite sequence of Hopf bifurcations with one complex conjugate pair of characteristic values crossing to the right half-plane as each curve Γ n is crossed in turn for n = 0, 1, 2 . . .. The curves Γ n do not intersect for γ > 0, so there are no double Hopf bifurcations.
Hence, by the Intermediate Value Theorem there exists (at least) one root of ϕ(y) = 0 for a 1 y ∈ (2n + 1 2 )π, (2n + 3 2 )π . Call this root y (n) and let κ (n) 1 be the corresponding value of κ 1 defined by (17). Thus we have found the existence of infinitely many Hopf bifurcation points. Moreover since (and in particular if N = 1, or if N = 2 and κ 2 is sufficiently small), and hence in this case there is a unique root y (n) in each interval. Moreover, it can be shown that for n sufficiently large there exists ξ ∈ (1/2, 1) depending on ε, γ, κ i , a i , N but independent of n, such that ϕ ′ (y) > 0 for a 1 y ∈ (2n + 1/2)π, (2n + ξ)π and ϕ(y) > 0 for a 1 y ∈ (2n + ξ)π, (2n + 3/2)π , and thus for general parameter values there is a unique y (n) in each interval for n sufficiently large. Now consider the behaviour of κ (n) as n → ∞. If n > a1 Hence ϕ(y (n) ) = 0 implies tan(a 1 y (n) ) < 0 using (18), and so Moreover, εy (n) → ∞ as n → ∞ implies that to satisfy ϕ(y (n) ) = 0 we require tan(a 1 y (n) ) → −∞, so a 1 y (n) → (2n + 1 2 )π. Thus sin(a 1 y (n) ) → 1 and κ . Hence for large n Hopf bifurcations obey Equation (21) implies that the period of the solution created in the Hopf bifurcation satisfies T (n) ∼ a 1 /(n + 1/4) for large n, so these solutions become more and more highly oscillatory. We will be more interested in solutions with large periods. When the trivial solution is unstable, Corollary 2.2 implies that any oscillatory solutions satisfy L 0 lim inf t→∞ u(t) < 0 < lim sup t→∞ u(t) M 0 . A more sophisticated solution bound can be found using a Gronwall argument. Suppose u(t 0 ) = v at a local minimum, and that But and substituting from (22) the minimum must satisfy But h(v, m) is continuous in v and using (5) and (6) h This bounds the solution u(t) away from −a 1 /c 1 and hence bounds the delays where This provides a bound on the amplitude of any periodic orbits, which we use below.
4. Oscillatory Dynamics. We now present numerical computations of the bifurcation structures, invariant sets and their stability for equation (1). These computations are performed with DDE-Biftool [9], which is able to detect Hopf bifurcations, switch to and follow branches of periodic orbits (arising from Hopf bifurcations or otherwise) and compute Floquet multipliers and hence stability. We emphasise that periodic orbits are found by solving a boundary value problem, and not by evolving the differential equation, so this package is able to follow branches of unstable periodic solutions as well as stable ones. Condition (5) leads us to use κ 1 as a bifurcation parameter, and we will consider the bifurcations that occur as κ 1 is varied, first with N = 1, and then with N = 2 and κ 2 initially small, but increasing it in subsequent examples. If for some parameters all the orbits found by DDE-Biftool are unstable, then the stable dynamics can be revealed by solving an initial value problem using the MATLAB [32] state-dependent DDE solver ddesd. In this way, we are able to find stable invariant tori. We refer to the branches of periodic orbits that bifurcate from the trivial solution in Hopf bifurcations as the primary branches, and to any branches of periodic orbits that bifurcate from the primary branches as secondary branches. The primary branch bifurcating from the first Hopf point (with smallest value of κ 1 ) is referred to as the principal branch. This first Hopf bifurcation is supercritical and the equilibrium solution loses stability to the periodic orbit in the bifurcation. Consequently, at least near the beginning of the principal branch, periodic orbits are stable. The results of Section 2 imply that u(t) ∈ C ∞ on any periodic orbit, but it is not known if these orbits are analytic in general.
We briefly consider the case of one delay, when (1) becomes In the rest of the paper we concentrate on the case of two delays under the restriction that c 1 = c 2 = c. We do not consider singularly perturbed problems in that case, and so without loss of generality we can divide equation (1) by ε and rescale γ and κ i accordingly. Thus we restrict attention to the case where ε = 1, and (1) becomeṡ where γ, κ i , a i and c are all strictly positive. By Theorem 2.3 the α i (t, u(t)) = t − a i − cu(t) are all monotonic functions of t on any invariant set, which will simplify the representation of the solutions of (26). Another reason for restricting attention to the case c = c 1 = c 2 is that this already leads to very rich dynamics, which we would like to understand before considering the general case. Parameters are also always chosen to be positive and satisfy (3), (5) and (14).
4.1. One Delay and Two Delays. In Figure 2 we present bifurcation curves for periodic solutions of (25). The periodic orbits on the branch bifurcating from the first Hopf bifurcation at κ ≈ 5.1922 are always stable (as indicated from the computed, but not shown, Floquet multipliers). The other Hopf bifurcations lead to periodic orbits that are always unstable (all these bifurcations occur when an already unstable equilibrium solution has an additional pair of complex conjugate characteristic values cross the imaginary axis from the left to right half-plane). The bifurcation points are given by (15). As expected from (21), the periods of the orbits decrease with each subsequent Hopf bifurcation; and these unstable orbits become increasingly highly oscillatory. Corollary 2.2 implies that the amplitude A = max t∈[0,T ] u(t) − min t∈[0,T ] u(t) of a periodic orbit of period T of (1) is bounded above by A < a1 c1 1 + 1 γ N i=1 κ i which is a linear function of κ 1 . This linear bound is indicated on Figure 2(i) as a straight line. This bound is not sharp (especially for κ small when the equilibrium solution is stable) but the amplitude of the stable branch of solutions does approach this bound as κ becomes large. The more sophisticated Gronwall bound on the amplitude M ∞ − L ∞ from (24) is indicated as a curve in Figure 2(i), and is seen to give a much better bound for the amplitude of the branch of stable periodic orbits.
In [29,31] slowly oscillating periodic solutions (SOPS) of (25) are studied in the singular limit ε → 0. A SOPS of (2) with N = 1 is a periodic solution with u(t) = 0 at t = T i for i ∈ N and α 1 (T i+1 , 0) > T i , so the delay T i+1 − α 1 (T i+1 , 0) at the (i + 1) st zero is smaller than the time interval T i+1 − T i between the i th and (i + 1) st zeros. For (25) we thus require T i+1 − T i > a for a SOPS. The SOPS of (25) is known [29,31] to have exactly one local maximum and minimum per period T with limiting profile, in the limit ε → 0, with a vertical transition from aκ/cγ to −a/c at t = T . Figure 3 shows the profiles of the stable periodic orbits of (25) on the principal branch bifurcating from κ ≈ 5.1922 with ε = 1. As κ → ∞ with ε = 1 the shape of the periodic orbit of (25) is seen to converge to the same saw-tooth shape (27) found by Mallet-Paret and Nussbaum [31] for the SOPS in the singular limit ε → 0 with γ and κ fixed. Rescaling (25) as we see that taking κ large is equivalent to taking both ε and γ small with κ fixed. Although the results of Mallet-Paret and Nussbaum apply for ε → 0 with γ fixed, nevertheless, numerically we still observe a SOPS of the form (27) with T large when we take κ large with ε = 1 and γ > 0 fixed. (Note that in Figure 3 the periods of the orbits are rescaled to 1; see Figure 2 for the actual periods.) By (20) the first Hopf bifurcation of (25) satisfies ay (0) < π (because for N = 1 we trivially have 0 > a1 2επ N i=2 κ i − 1 4 ). Hence, the period T 0 of the periodic orbit of (25) created at the first Hopf bifurcation satisfies  and so this periodic solution is already a SOPS at the Hopf bifurcation. It also converges to a singularly perturbed SOPS as κ becomes large, and numerically the periodic orbit is seen to remain a SOPS along the whole of the principal branch.
As κ becomes large, the unstable periodic orbits on the other branches are seen to converge to similar looking profiles defined on the n th branch by Notice that all these periodic solutions have the same limiting slope 1/c in the slowly evolving part of the orbit, but for higher Hopf bifurcation numbers n the periods become ever smaller, and none of these orbits for n 1 are SOPS. These orbits are constructed by looking for a solution where α 1 (t, u(t)) falls in one of the preceding vertical transition layers while t ∈ (0, T ), so t − a − cu(t) = −nT for t ∈ (0, T ). A rigorous derivation for the case n = 0 can be found in [29,31]. Not surprisingly, similar bifurcation diagrams are observed for (26) if we choose the coefficient of the second delay κ 2 to be small. Letting c 1 = c 2 = c and a 2 > a 1 (to satisfy (3)) we find that the same picture also persists qualitatively for κ 2 ≫ 0 if a 2 − a 1 is small. We consider the case a 2 − a 1 ≫ 0 below.

4.2.
Bistability of Periodic Solutions. Consider (26) with a 2 > a 1 > 0 chosen so that a 2 −a 1 > T 0 , the period of the stable orbit at the first Hopf bifurcation. If κ 2 is small, the behaviour is similar to that seen for (25) in Figure 2. However, for κ 1 sufficiently large two saddle node bifurcations are created on the principal branch of periodic solutions and a region of bistability of periodic solutions is created. This is illustrated in Figure 4(i) for κ 2 = 1 with a 2 − a 1 = 4.7 and T 0 ≈ 2.86. We see that for each κ 1 with 10.643 κ 1 11.0907 there are three periodic orbits on the principal branch, two stable and one unstable. These orbits are shown as functions of time over one period in Figure 4(ii) for κ 1 = 10.95 when the stable orbits (as indicated by the Floquet multipliers) have periods T 1 ≈ 4.1954 and T 3 ≈ 4.7262, and the unstable orbit has period T 2 ≈ 4.501, with T 1 < T 2 < T 3 . Figure 4(ii) also shows two non-periodic orbits computed using ddesd by taking very small perturbations of the unstable orbit, then plotting the slices of the resulting solution trajectory between each downward crossing of u = 0. We see that these orbits converge to the stable periodic orbits, which independently confirms the stability of those orbits. The two nonperiodic orbits lie in the unstable manifold of the unstable periodic orbit, but this manifold is difficult to represent graphically since periodic orbits of (1) with c i = c for all i and their unstable manifolds lie in the infinite-dimensional function space defined by (10). In Figure 5 we present two projections of the periodic orbits into R 3 that reveal some of the structure of the unstable manifold of the unstable periodic orbit which lies between the two stable orbits. In Figure 5(i) we plot solutions in R 3 using the values of the solution at the current and delayed times; (u(t), u(α 1 (t, u(t)), u(α 2 (t, u(t))). For problems with a single fixed or prescribed delay, it has become common to plot u(α(t)) against u(t) since the work of Mackey and Glass [24] showed this projection to yield insight into the orbit's shape. Plotting (u(t), u(α 1 (t, u(t)), u(α 2 (t, u(t))) is a natural generalisation of this, where we note, by Theorem 2.3, that the α i (t, u(t))  Figure 4 as (i) (u(t), u(α 1 (t, u(t)), u(α 2 (t, u(t))) and (ii) (u(t),u(t),ü(t)). The stable periodic orbits (red) appear to define the boundaries of a surface which the unstable periodic orbit (black) and its unstable manifold lie on. The two non-periodic orbits (green) both lie in the unstable manifold of the unstable periodic orbit, and the stable manifold of one of the stable periodic orbits.
are strictly monotonically increasing functions of t; this projection would be problematical if they were not. Since an analytic function is completely described by its value and derivatives, another natural projection to finite dimensions is to plot the coefficients of the Taylor polynomial of relevant degree, and in Figure 5(ii) we plot (u(t),u(t),ü(t)) in R 3 . Both projections clearly show the unstable orbit and its unstable manifold filling a surface between the two stable periodic orbits.  Figure 4. Also shown is the line T = a 2 − a 1 that separates solutions for which the difference between the delays is more than or less than one period.
For a possible explanation of the bistability consider Figure 6, which shows the periods of the orbits. Notice that T < a 2 − a 1 = 4.7 on the first stable segment of the principal branch, while T > a 2 − a 1 = 4.7 on the second stable segment, and that by (8) we have α 1 (t, u(t)) − α 2 (t, u(t)) = a 2 − a 1 . Although there is no general definition of a SOPS for a problem with multiple delays, and the construction (27) does not apply, we see somewhat similar behaviour with the periodic orbit developing 'sawtooth' shape for κ 1 large, and the solution evolving with α 1 (t, u(t)) in the sharp transition layer. Now if a 2 − a 1 is close to T , there are two possible cases. If a 2 − a 1 < T then α 1 (t, u(t)) − α 2 (t, u(t)) < T so when α 1 (t, u(t)) is in the transition layer u(α 2 (t, u(t))) < 0. If a 2 − a 1 > T then α 1 (t, u(t)) − α 2 (t, u(t)) > T so when α 1 (t, u(t)) is in the transition layer u(α 2 (t, u(t))) > 0. It should not be surprising that changing the sign of u(α 2 (t, u(t))) over a large part of the period would have a significant effect on the solution, and this interaction between the two delay terms is implicated in the bistability. A more systematic construction of the singularly perturbed solutions in this case is a current research topic.

Stable Tori.
As well as bistability of periodic orbits, we also find regions of parameter space where the trivial solution is unstable, but there is no stable periodic orbit. Maintaining the parameter values of Section 4.2, but increasing κ 2 > 2, we find first one then two intervals of κ 1 values for which there is no stable periodic orbit. For κ 2 = 2.3 the principal branch of periodic orbits loses stability in two intervals, approximately κ 1 ∈ (5.377, 6.6793) and κ 1 ∈ (7.0806, 7.6939) with a complex Floquet multiplier of modulus 1 at the end point of each interval (λ ≈ 0.014411 ± 0.99989i at κ 1 = 5.377 and λ ≈ 0.61505 ± 0.78849i at κ 1 = 6.6793). In general, this indicates a torus bifurcation otherwise known as a Neimark-Sacker bifurcation [22]. Although there is no rigorous theory for torus bifurcations in statedependent DDEs and no general algorithm for finding invariant tori in the (infinitedimensional) phase space of state-dependent DDEs, by Corollary 2.2 solutions of (1) remain bounded so there must be some stable invariant object in the flow after the bifurcation. The simplest scenario would be a supercritical bifurcation to a stable torus, and below we show numerically that this indeed occurs.
In parameter regions where there is no stable periodic orbit on any of the primary branches of periodic solutions created in the Hopf bifurcations, we identify the stable dynamics in the flow using a small perturbation of the unstable periodic orbit on the principal branch found by DDE-Biftool as the initial history function for an initial value problem integration of (1) using ddesd. Discarding the transient dynamics reveals the stable invariant object. The results of this computation are shown in Figure 7(i) with parameter values κ 1 = 5.95 and κ 2 = 2.3. This reveals a very structured object, which in the (u(t), u(α 1 (t, u(t)), u(α 2 (t, u(t))) projection looks very like a classical torus.
In Figure 7(ii) we plot a representation of a Poincaré section of the dynamics to verify the existence of the torus seen in Figure 7(i). In Section 2 we showed that the largest time interval on which u(t) andu(t) have the same sign is bounded above by max i {a i }. This implies that all orbits are either eventually monotonic (and converging to u = 0) or cross u = 0 infinitely many times, but eventually monotonic orbits do not occur because the equilibrium solution u = 0 does not have any real negative characteristic values with the chosen parameter values. Thus, any non-trivial orbit crosses u = 0 infinitely often, and this is a natural surface to consider for a Poincaré section. This Poincaré section is infinite-dimensional since the phase space of the dynamical system is defined by (10), but plotting function segments is not very revealing (especially as on this Poincaré section they all have value u = 0 at the left-hand end) and so we project the Poincaré section into R 2 . In Figure 7(ii) at each time that u(t) = 0 withu(t) < 0 we plot the point u(α 1 (t, u(t))) u(α 2 (t, u(t))) (i) (ii) Figure 7. (i) (u(t), u(α 1 (t, u(t)), u(α 2 (t, u(t))) projection of a single trajectory (in blue) of (26) filling a stable torus, computed using ddesd, for κ 1 = 5.95, κ 2 = 2.3, c 1 = c 2 = 1, γ = 4.75, a 1 = 1.3 and a 2 = 6. The unstable periodic orbit on the principal branch found by DDE-Biftool for the same parameters is also shown (red).
(ii) Projection of the Poincaré section for the torus in (i). The (blue) dots are points of the trajectory on the torus obtained by plotting (u(α 1 (t, u(t))), u(α 2 (t, u(t)))) at each time that u(t) = 0 withu(t) < 0. The (red) star represents the unstable periodic orbit.
(u(α 1 (t, u(t))), u(α 2 (t, u(t)))) (equivalently (u(t − a 1 ), u(t − a 2 )) since u(t) = 0) in R 2 , revealing a smooth closed curve. In simple textbook examples in R 3 , this would be a simple connected curve in R 2 with the periodic orbit in its interior, whereas the curve in Figure 7(ii) has a point of self intersection. This point of self intersection arises from the projection to R 2 , which can be chosen in many different ways, and the projection to R 2 using the delayed values, as it turns out, does not reveal a familiar projection of the torus. Nevertheless, since the point of self intersection is due to projection, Figure 7(ii) clearly reveals that there is a stable smooth invariant torus with quasi-periodic dynamics or periodic dynamics with very high period. Figure 8 shows the bifurcation diagram as κ 1 is varied for κ 2 = 2.3 and the other parameters as previously stated, with the branches of periodic orbits and their stability computed using DDE-Biftool, while the tori are computed using ddesd as described above for each set of parameters values. The amplitude of the torus for the bifurcation diagram is computed as A = max t∈[0,T ] u(t) − min t∈[0,T ] u(t) for an orbit that makes 400 revolutions around the torus. There are four torus bifurcations on the principal branch of periodic orbits with κ 1 < 9 (detected by computing the Floquet multipliers of the periodic orbit). These torus bifurcations border two intervals where the periodic orbit on the principal branch is unstable. In these two intervals, κ 1 ∈ (5.377, 6.6793) and κ 1 ∈ (7.0806, 7.6939), the amplitude of the branch of torus solutions is shown, and we see that these tori both bifurcate from the stable periodic orbit at each end of these intervals, and do not persist outside these parameter intervals, so each of these torus bifurcations is indeed supercritical. As for the case κ 2 = 1, we see a region of bistability of periodic orbits when κ 2 = 2.3 between the saddle node bifurcations of periodic orbits in the approximate interval κ 1 ∈ (9.1800, 10.0483). However, as the inset zoom in Figure 8 shows, the bifurcation diagram is more complicated than in the previous example, with additional torus bifurcations and saddle-node bifurcations of periodic orbits indicated. DDE-Biftool also detects two pairs of Floquet multipliers crossing the unit circle at two points on the second branch of periodic solutions, indicating additional torus bifurcations. However, these torus bifurcations are from unstable periodic orbits, and we expect unstable tori to be created in these bifurcations. Since there is yet no method for computing unstable tori in state-dependent DDEs, we are unable to investigate these tori.

Other Bifurcations and Stable Solutions on Secondary Branches.
Maintaining the values of the other parameters from Section 4.2, but increasing κ 2 to κ 2 = 3 reveals additional bifurcations as κ 1 is varied, as shown by the bifurcation diagram in Figure 9. There is still a region of bistability of periodic orbits on the principal branch for κ 1 ∈ (7.82, 8.2585), but now there is a much larger part of the principal branch where the periodic solution on this branch is unstable, and for some parameter values we also encounter stable periodic solutions that are not on one of the primary branches of periodic solutions created in the Hopf bifurcations. For κ 1 ∈ (9.086, 9.363) the principal branch of periodic orbits loses stability in a pair of period-doubling bifurcations to a stable branch of period-doubled solutions, as shown in Figure 10. The stable period-doubled branch only exists for a small range of parameter values, and the periodic orbits on the branch are seen in Figure 10(ii) to all be close to two copies of the periodic orbit on the principal branch, but as shown by Figure 10(i), with a slightly increased amplitude. These stable period-doubled orbits are our first example of stable periodic orbits on a secondary bifurcation branch for (26). This shows that the dynamics in the multiple-delay case is richer and fundamentally different from the single delay case (25) where stable periodic orbits are only found on the principal primary branch of periodic solutions emanating from the first Hopf bifurcation. The existence of period-doubled solutions, opens the question as to whether (26) could exhibit a period-doubling cascade in some parameter range. However, for the parameter values considered here, computation of the Floquet multipliers reveals that the period-doubled branch is stable with no secondary bifurcations between the two period-doubling bifurcations.
The first Hopf bifurcation at κ 1 ≈ 3.2061 is found to be supercritical, as in all the previous cases, but the stable periodic orbit thus created soon loses stability in a torus bifurcation at κ 1 ≈ 3.656 (with characteristic values λ ≈ −0.44504 ± 0.89548i). Unlike the previous example shown in Figure 8 where the stable torus soon contracted back on to the periodic orbit on the principal branch, the stable torus created at κ 1 ≈ 3.656 is seen in Figure 9 to exist and be stable for a large parameter interval up to κ 1 ≈ 7.6818. There are three parameter intervals along the branch of stable tori for which we find phase-locked periodic orbits on the torus that persist for a significant parameter interval, namely for κ 1 ∈ (5.6727, 5.7937), κ 1 ∈ (6.8252, 6.9763) and κ 1 ∈ (7.5796, 7.6818). The behaviour in the second of these three intervals is illustrated in Figure 11. While the torus itself remains stable for these parameter values, at κ 1 ≈ 6.8252 a stable and unstable periodic orbit are created on the torus through a saddle node bifurcation of periodic orbits, and then destroyed in a second saddle node bifurcation at κ 1 ≈ 6.9763. For intervening parameter values there is exactly one stable and one unstable periodic orbit on the torus, each of which form a (1:4)-torus knot with the two knots interleaved around the torus. A closeup of the bifurcation diagram for these periodic orbits is shown in Figure 11(i). There are no secondary bifurcations, and so these periodic orbits form an isola (isolated closed curve) on the bifurcation diagram of periodic orbits, and cannot be reached by following primary and secondary and subsequent bifurcations of periodic orbits starting from the equilibrium solution. Thus the upper stable branch of the isola in Figure 11(i) gives another example of a stable periodic solution that is not on the principal branch; in contrast to the period-doubled orbit found earlier, the stable periodic orbit on the torus is an example of a stable periodic orbit that cannot be found by simply following all the bifurcations of periodic solutions starting from the trivial solution. The existence of an isola of periodic solutions for (26), is interesting in the context of Wright's conjecture [23] where the possibility of the existence of an isola of SOPS is one of the stumbling blocks in establishing the conjecture. The profiles of a number of the periodic orbits on the isola are overlayed in Figure 11(ii). These were computed starting at κ 1 = 6.9 on the stable branch of the isola and performing one circuit of the isola. The periodic orbits vary continuously as the isola is traversed, but the phase shifts; the first and last computed stable periodic orbits in the figure have the same profile, but are phase shifted.
These periodic orbits and the isola were found by performing an initial value problem integration of (1) using ddesd to reveal the dynamics on the torus. For parameter values for which there is no phase-locking the resulting orbit fills the torus and its 'amplitude' as displayed in Figure 9 is computed as A = max t∈[0,T ] u(t) − min t∈[0,T ] u(t). For parameters values for which there is phase-locking the ddesd computed orbit, no longer fills the torus, but instead converges to the stable periodic orbit on the torus. Hence it is not possible to numerically compute a torus 'amplitude' for these parameter values, which is why the torus curve in Figure 9 has gaps. However, having found a stable periodic orbit on the torus by initial value integration, the same solution can be imported into DDE-Biftool and then followed to obtain the resulting branch of periodic orbits; its Floquet multipliers can also be computed to confirm its stability. In this way we find the two saddle-node bifurcations of periodic orbits and the counterpart unstable periodic orbit on the torus (which could not otherwise have been found, being unstable) and confirm that there are no other bifurcations from the isola of periodic orbits. We mention that in the inset in Figure 9, near κ 1 = 5.6727 and κ 1 = 5.7937 the amplitude of the torus is not continuous with that of the phase locked periodic orbits created in the saddle node bifurcations on the torus. This is because the amplitude of the torus is measured from an orbit that fills the torus, whereas the phase locked periodic orbits lie on the torus but do not fill it, and so in general will have a smaller amplitude. Figure 12 illustrates the different dynamics on the torus inside and outside the phase-locked parameter regions with the distinct stable and unstable periodic orbits visible for κ 1 = 6.9 while for κ 1 = 7 a single orbit fills the torus. Similar dynamics are observed for the isola of periodic orbits in the first interval κ 1 ∈ (6.8252, 6.9763), except that the periodic orbit closes after three revolutions of the torus, not four. However the behaviour is somewhat different in the final interval κ 1 ∈ (7.5796, 7.6818) where the stable torus coexists with a stable periodic orbit on the principal branch of periodic solutions. The structure of the principal branch of periodic solutions is quite complicated in this region of parameter space with several saddle node of periodic orbit bifurcations and several torus bifurcations, as shown in Figure 9. As for the stable torus, numerics indicates that there may be a secondary bifurcation and break up of the torus. Torus break-up in this and other state-dependent DDEs is an interesting topic for future research. Tori are often associated with double Hopf bifurcations [14,22], so in Figure 13 we plot the loci of the Hopf bifurcations as κ 1 and κ 2 are varied. This reveals four double Hopf bifurcations which all occur with κ 2 ∈ (γ/2, γ) and three of which lie on a branch of Hopf bifurcations that does not intersect the κ 1 axis. This branch only exists for κ 2 > 2.627, and manifested itself already in Figure 9 where for κ 2 = 3 we saw a bounded branch of periodic orbits joining Hopf bifurcations from the trivial solution at κ 1 = 3.973 and κ 1 = 7.46. These three double Hopf bifurcations occur with κ 2 > 2.6436, and yet we already observed stable tori for κ 2 = 2.3, which is not only smaller than the values of κ 2 for which the double Hopf bifurcations occur, but smaller than the values of κ 2 for which the branch they occur on even exists. Thus if the tori are created at these double Hopf bifurcations, they persist for smaller values of κ 2 than the double Hopf bifurcations. If the dynamics are studied by gradually increasing κ 2 , as we have done, tori can be found before discovering the codimension two bifurcation that gives rise to them. Unfortunately, DDE-Biftool does not have a routine for following curves of torus bifurcations to confirm this, but the usual unfolding of the double Hopf bifurcation [22], although not yet analyzed in the case of state-dependent delay, and general bifurcation theory leads us to expect this leaching of the effects of bifurcations to nearby regions of parameter space.
Notice that equation (21) gives an asymptotic formula for the values of κ 1 at the Hopf bifurcations when κ 1 is large, This shows that when κ 2 = 0 and κ 1 ≫ 0, the bifurcations are spaced at intervals of 2επ/a 1 , while κ (n) 1 and κ 2 are linearly related for κ 2 = 0. Thus, for κ 1 sufficiently large there are no double Hopf bifurcations for κ 2 < επ/a 1 , which is κ 2 < 2.4166 for the parameters in our example. The double Hopf bifurcation at (κ 1 , κ 2 ) ≈ (47.64, 2.6063) respects these bounds. However, possible double Hopf bifurcations and invariant tori for κ 1 large are not very interesting, as the corresponding frequencies from the linear theory reveal them all to be highly oscillatory, and numerics do not reveal any highly oscillatory dynamics which is stable. Finally, Figure 14 shows the bifurcation diagram for the periods of the periodic orbits as κ 1 is varied when (i) κ 2 = 3, and (ii) κ 2 = 3.1. The locations of the Hopf bifurcations and the periods and stability of the branches near these bifurcations are very similar in both cases, but the global bifurcation structure appears very different, as the branches connect to different Hopf points. Most likely, there is either a global bifurcation that changes the branch connectivity as κ 2 varies between 3 and 3.1, and/or there are other missing branches of periodic orbits in our bifurcation diagram (adding such branches could make the two diagrams topologically equivalent). The bifurcation structures when κ 2 is large is a topic for future research.

5.
Conclusions. In this work, we found a condition on the parameters (5) which ensures existence of solutions to the state-dependent DDE (1) and in particular that α i (t, u(t)) < t for all t > 0. We then used DDE-Bifftool and ddesd to numerically explore the non-singularly perturbed two delay problem (26). In that problem we found a multitude of bifurcations including saddle-node bifurcations, perioddoubling bifurcations, torus bifurcations, an isola of periodic solutions, double Hopf bifurcations and a possible global bifurcation which changes the branch connectivity. Such bifurcation structures are familiar from ODEs and fixed delay DDEs, but a rigorous theory for such bifurcations has yet to be developed for state-dependent DDEs. The fact that we are able to numerically compute the bifurcations suggests that it should be possible to extend the DDE bifurcation theory to state-dependent DDEs, and Sieber [35] suggests a method for proving such results.
The saddle-node bifurcations of periodic orbits give rise to parameter regions with bistability of periodic orbits on the principal branch. Period-doubling bifurcations result in stable periodic solutions on a secondary branch of periodic orbits. The torus bifurcations lead to parameter regions where the equilibrium solution and the periodic orbits on the primary, secondary and subsequent branches were all unstable, but the dynamics are contained a stable torus. Locating intervals of phase locking on a stable torus allowed us to locate a branch of stable periodic solutions on an isola of periodic solutions in the bifurcation diagram.
All of these dynamics were observed in (26), which is not singularly perturbed, has two linearly state-dependent delays, with the difference between the delays fixed, and no other nonlinearity. Since the only nonlinearity in (26) arises through the state-dependency of the delays we conclude that for equations with multiple delays, state-dependency of the delays on its own is sufficient to generate very complex dynamics. Consequently, neglecting the state-dependency of delays in mathematical models has the potential to suppress interesting dynamics that are inherent to the system being modelled. The linear state-dependency that we consider is, of course, just the first two terms of a Taylor expansion of a general state-dependent delay, and we would expect similar dynamics to be found in related multiple statedependent DDE models. It would be interesting therefore to study state-dependent versions of DDEs arising in applications where the state-dependency of the system has previously been excluded from the model.
Our study of this equation suggests a number of other possible avenues for future research. It would be interesting to study the dynamics of the tori encountered in Sections 4.3 and 4.4 including the structure of the Arnold tongues that the phase locked orbits lie in and the possible break-up of the torus. The torus dynamics are of particular interest, as they have not, to our knowledge, been studied for statedependent DDEs. Torus dynamics and break-up is studied in [13] for a fixed delay DDE arising from a laser model. Numerical analysis issues in the computation of invariant tori for state-dependent DDEs are also of interest. While torus break-up is one potential route to chaos in (26), other avenues, including the possibility of a period-doubling cascade are also worth investigating. It would also be interesting to extend the singularly perturbed analysis of Mallet-Paret and Nussbaum [31] to the case of (1) with two delays. The results of Section 4.2 on the bistability of periodic solutions, and the dependence of the behaviour on the difference a 2 − a 1 between the delays suggests that the results would be much more complicated than in the case of (25). Even in the nonsingularly perturbed limit, studying (26) when (5) does not hold, or considering (1) with two delays and ε = 1 but c 1 = c 2 could lead to new and interesting dynamics. The dynamics for more than two delays has also yet to be explored. 6. Acknowledgments. Tony Humphries is grateful to John Mallet-Paret and Roger Nussbaum for introducing him to this problem, and patiently explaining their results for the in the case N = 1. He is also grateful to NSERC (Canada) for funding through the Discovery Award program. We are grateful to Dave Barton for his assistance in computing the branch of period-doubled solutions seen in Figures