{'fraction_non_alphanumeric': 0.05528032942667089, 'fraction_numerical': 0.03325942350332594, 'mean_word_length': 3.7229968284363593, 'pattern_counts': {'":': 0, '<': 70, '<?xml version=': 0, '>': 138, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Periodic Solutions of a Singularly Perturbed Delay Differential Equation With Two State-Dependent Delays

Periodic orbits and associated bifurcations of singularly perturbed state-dependent delay differential equations (DDEs) are studied when the profiles of the periodic orbits contain jump discontinuities in the singular limit. A definition of singular solution is introduced which is based on a continuous parametrisation of the possibly discontinuous limiting solution. This reduces the construction of the limiting profiles to an algebraic problem. A model two state-dependent delay differential equation is studied in detail and periodic singular solutions are constructed with one and two local maxima per period. A complete characterisation of the conditions on the parameters for these singular solutions to exist facilitates an investigation of bifurcation structures in the singular case revealing folds and possible cusp bifurcations. Sophisticated boundary value techniques are used to numerically compute the bifurcation diagram of the state-dependent DDE when the perturbation parameter is close to zero. This confirms that the solutions and bifurcations constructed in the singular case persist when the perturbation parameter is nonzero, and hence demonstrates that the solutions constructed using our singular solution definition are useful and relevant to the singularly perturbed problem. Fold and cusp bifurcations are found very close to the parameter values predicted by the singular solution theory, and we also find period-doubling bifurcations as well as periodic orbits with more than two local maxima per period, and explain the alignment between the folds on different bifurcation branches.


Introduction
We consider singularly perturbed periodic solutions of the scalar state-dependent delay-differential equation (DDE) εu(t) = −u(t) − K 1 u(t − a 1 − cu(t)) − K 2 u(t − a 2 − cu(t)), (1.1) which has two linearly state-dependent delays, and no other nonlinearity other than the state-dependency of the delays. We consider ε 0, c > 0, a i > 0, K i > 0, and without loss of generality we order the terms so that a 2 > a 1 > 0. Equation (1.1) is an example of a singularly perturbed scalar DDE with N state-dependent delays of the form εu(t) = f (t, u(t), u(α 1 (t, u(t))), . . . , u(α N (t, u(t)))), u(t) ∈ R. (1.2) We will define a concept of singular solution for (1.2) based on a continuous parametrisation. This essentially entails defining a singular limit for the equation (1.2), resulting in an equation whose solutions can in principle be found algebraically. In the case of (1.1) we construct several such classes of singular periodic solution, and investigate the bifurcations that arise. Delay differential equations arise in many application fields including engineering, economics, life sciences and physics [12,24,36]. There is a well established theory for functional differential equations as infinite-dimensional dynamical systems on function spaces [9,16], which encompasses delay differential equations with constant or prescribed delay. However, many problems that arise in applications have delays which depend on the state of the system (see for example [13,20,23,37]). Such state-dependent DDEs fall outside of the scope of the previously developed theory and have been the subject of much study in recent years. See [17] for a relatively recent review of the general theory of state-dependent DDEs.
The study of singularly perturbed DDEs already stretches over several decades. As early as 1985 Magalhaẽs [25] recognised that singularly perturbed discrete DDEs have different asymptotics to singularly perturbed distributed DDEs. For equations with a single constant delay, in the singular limit the DDE reduces to a map (see (1.7) below) which describes the asymptotic behaviour when the limiting profiles are functions [21,34].
One of the principal difficulties in studying (1.2) in the singular limit is that while the solution u(t) is a graph for any ε > 0, this need not be so in the limit as ε → 0, when derivatives can become unbounded, and the resulting limiting solution can have jump discontinuities. Techniques for studying singularly perturbed DDEs with a single constant discrete delay can be found in [26,6]. In [26] slowly oscillating periodic solutions (SOPS) are proved to converge to a square wave in the singular limit, using layer equations to describe the solution in the transition layer. In [6] for monotone nonlinearities a homotopy method is used to show that the layer equations have a unique homoclinic orbit. Mallet-Paret and Nussbaum, in a series of papers [27,29,30] extend the study of SOPS to DDEs with a single state-dependent delay. In [27] SOPS are established to exist for all ε sufficiently small. These solutions are shown to have non-vanishing amplitude in the singular limit in [29], and under mild assumptions the discontinuity set of the limiting profile is shown to consist of isolated points. In [30] Max-plus operators are introduced to study the shapes of the limiting profiles. The DDE is considered as a particular example in [30]. This corresponds to (1.1) with K 2 = 0 and c = 1. It is established in [30] that the limiting profile is the "sawtooth" shown in Figure 1.1(ii) below. In [31] the SOPS of (1.3) are studied in detail and the shape of the solution near the local maxima and minima is determined for 0 < ε 1 as well as the width of the transition layer, and the "super-stability" of the solution. Other work on singularly perturbed state-dependent DDEs includes [14] where state-dependent DDEs arising from the regularisation of neutral state-dependent DDEs are studied, and also [33] where metastability of solutions of a singularly perturbed state-dependent DDE is studied in the case where the state-dependency vanishes in the limit as ε → 0.
The studies mentioned above all considered singularly perturbed DDEs in the case of one delay only, and either considered single solutions or a sequence of solutions as ε → 0. In the current work we will study the bifurcation diagram for (1.3) when 0 ε 1, regarding K 1 as a bifurcation parameter. Beyond the works mentioned previously, the only other work of which we know that tackles singularly perturbed bifurcations is [22], where the solutions of (1.3) with a 1 = c = 1 are studied close to the singular Hopf bifurcation. On the other hand, singularly perturbed ODEs frequently arise through mixed mode oscillations on multiple time-scales and their bifurcation analysis is well understood (see [8] for a review). Co-dimension two bifurcations have also been studied in singularly perturbed ODEs [3,5].
The development of bifurcation theory for state-dependent DDEs has been difficult because the centre manifolds have not been shown to have the necessary smoothness [17], and a rigorous Hopf bifurcation theorem for state-dependent DDEs was first proven only in the last decade [10] (see also [18,35,15]). The numerical analysis of state-dependent DDEs is more advanced with numerical techniques for solving both initial value problems [1,2] and for computing bifurcation diagrams [11]. Although it was developed before the Hopf bifurcation theorem was proved, DDEBiftool [11] is a very useful tool for computing Hopf bifurcations and continuation of solution branches in state-dependent DDEs, and it has been used to study the bifurcations that arise in (1.1) when ε = 1 [19,4]. John Mallet-Paret has presented numerical simulations of (1.1) in seminars, but the only other published work of which we are aware that encompasses (1.1) is [28]. There the existence of SOPs was proved for (1.2) with suitable nonlinearities when α i (t, u(t)) = t − τ i (u(t)) with τ i (0) = 0 for all i.
In [19] a largely numerical investigation of (1.1) with ε = O(1) revealed fold bifurcations on the branches of periodic orbits, resulting in parameter regions with bistability of periodic orbits. While the stable periodic orbits usually had one local maxima per period, the unstable periodic orbits in the these windows of bistability often had more than one local maxima per period. In the current work we will investigate these fold bifurcations and the profiles of the periodic orbits in the singular limit ε → 0.
To study (1.2) in the singular limit ε = 0 when the limiting profile may have jump discontinuities, we propose nested continuous parameterisations of the limiting singular solution. We will not restrict our attention to slowly oscillating periodic orbits, but will consider both long and short period orbits. We will study the case of two delay state-dependent DDE (1.1) in detail, and construct branches of singular periodic orbits with fold and cusp bifurcations. We will then use the predictions of this theory to guide a numerical study which will reveal branches of periodic orbits for 0 < ε 1 with profiles close to the singular limiting profiles and fold and cusp bifurcations close to the predicted parameter values. We will also find period-doubling bifurcations in the singularly perturbed problem. For our outer parametrisation we consider the solution profile as a parametric curve, Γ (µ) = (t(µ), u(µ)). This is a familiar concept from physics, where trajectories in space-time are parameterised, and has been used in the study of the DDEs arising from Wheeler-Feynman Electrodynamics [7]. However, in the current work we use the parametric curve Γ (µ) to enable us to consider continuous objects even in the singular limit. For any ε > 0 an injective parametrisation of the solution must have t(µ) strictly monotonic, but limiting profiles as ε → 0 may have t(µ) merely monotonic. This leads us to the parametric definition of an admissible singular solution profile in Definition 1.1. In Definition 1.2 we will introduce the inner parametrisation that allows us to define singular solutions of (1.2). Definition 1.1 Let Γ : I → R 2 be a continuous injective parametric curve defined on a nonempty interval I ⊆ R. For µ ∈ I let Γ (µ) = (t(µ), u(µ)). Then if t : I → R is monotonically increasing we say that Γ (I) is an admissible singular solution profile for (1.2).
Although t(µ) is not required to be a strictly monotonically increasing function to be an admissible singular solution profile, it is important to note that on any subinterval I i on which t(µ) is constant, the injectivity requirement ensures that u(µ) is strictly monotonic. Thus we partition the interval I as I = I + ∪ I − ∪ I * where 1. I * a disjoint union of open intervals and t(µ) is strictly monotonically increasing on each interval, 2. I ± are each disjoint unions of closed intervals with t(µ) constant on each such interval, and u(µ) strictly monotonically decreasing (respectively increasing) on each interval of I − (resp. I + ).
The partition of I generates a corresponding partition of Γ (I) as Γ (I) = Γ + ∪ Γ − ∪ Γ * . For (1.1) we will find that I + = / 0, and so I * and I − will both be unions of disjoint intervals which we may write as for a sequence of strictly increasing real numbers b i . See Fig. 1.1 for an example.
If J = R and there exists T > 0 and η T > 0 such that t(µ i (η + η T )) = t(µ i (η)) + T, u(µ i (η + η T )) = u(µ i (η)), ∀i = 1, . . . , N, ∀η ∈ R, (1.8) then we say that the singular solution is periodic. The period is the smallest T > 0 for which (1.8) holds. The main aim of this paper is to initiate a study of periodic solutions of the singularly perturbed two-delay DDE (1.1). We will construct singular periodic solutions (as per Definition 1.2), and will find both unimodal sawtooth solutions that correspond to the profile seen in Fig. 1.1 and bimodal solutions which have two "teeth" per period. The labels unimodal and bimodal are used throughout this work to indicate the number of local maxima of the solution per period. Although superficially the unimodal solutions look similar to those found in the one delay case, the interaction between the two statedependent delays adds both complications to the derivations and richness to the dynamics observed. We will demonstrate numerically using DDE-Biftool [11], a sophisticated numerical bifurcation package for delay differential equations, that the singular solutions and associated bifurcation structures that we find persist for ε > 0.
In Section 2 as an example we first consider (1.3) with one delay, for which Mallet-Paret and Nussbaum [30,31] have already established the so-called sawtooth limiting profile, as illustrated in Fig. 1.1(ii). We construct the corresponding singular solution following Definition 1.2. We then consider the two-delay problem (1.1) and in Theorem 2.2 establish conditions on the parameters for this to have a sawtooth solution. In (2.16) and (2.17) we introduce two admissible singular solution profiles which have two local maxima per period. Theorems 2.3 and 2.4 present singular solutions for these profiles and establish the constraints on the parameters for them to exist. Since these solutions have two local maxima per period, we refer to them as type I and type II bimodal (periodic) solutions.
In Section 3 we will treat K 1 as a bifurcation parameter and in Theorems 3.1, 3.3 and 3.4 identify intervals of the parameter K 1 for which unimodal, type I bimodal and type II bimodal solutions exist. We will also find singular fold bifurcations in Theorem 3.3 where solutions transition between unimodal and type I bimodal solutions. Theorem 3.4 as well as identifying a singular fold bifurcation between the unimodal and type II bimodal solutions also identifies a curve of parameter values at which a co-dimension two singular cusp bifurcation occurs. At this bifurcation the co-dimension one fold bifurcation unfolds and there is a transition between unimodal and type II bimodal solutions without a fold in the bifurcation branch.
The definition of singular solution introduced above, and the resulting solutions found are only useful if they tell us something about the dynamics of (1.1) when 0 < ε 1. In the case of one delay (1.3), Mallet-Paret and Nussbaum [30] proved the existence for ε > 0 of a singular solution which is a perturbation of the sawtooth profile. It is not readily apparent how to extend that proof to the two delay DDE (1.1). So in Section 4 we perform an extended numerical investigation of (1.1) with 1 ε > 0 close to the singular limit. We use DDEBiftool [11] to construct bifurcation diagrams and find periodic solutions in various parameter regimes guided by our results from Sections 2 and 3. In particular we show numerically that there are periodic solutions of (1.1) for 0 < ε 1 which are perturbations of the unimodal and type I and II bimodal solutions that we constructed in Section 2. Moreover, we find fold bifurcations close to the values predicted by our singular solutions and the co-dimension two bifurcation is also found for ε > 0. We also numerically investigate multimodal solutions, which are more complex than the singular solutions that we constructed algebraically. The existence of these seems to be generic on the unstable legs of the bifurcation branches between folds. On the principal branch stable periodic solutions with two or more local maxima per period are found through period doubling bifurcations, cusp bifurcations and in parameter regimes where stable unimodal singular solutions cannot occur. We also consider the alignment of the fold bifurcations on different solution branches and explain this using our results from Section 3. We finish in Section 5 with brief conclusions.

Singular Solutions
Before constructing singular solutions for (1.1), as an example we consider the singular solutions of the one delay DDE (1.3) which we write as We will construct periodic singular solutions following Definition 1.2 for (2.1) when K > 1 (required for instability of the trivial solution), with the following profile.
Definition 2.1 (Sawtooth Profile) For any n ∈ N 0 and period T > 0 the sawtooth profile is an admissible periodic singular solution profile on I = R defined by for each i ∈ Z.
Fig. 1.1 shows a part of this profile when a 1 = c = 1. Notice that I * is the union of the intervals (2i, 2i + 1) and on each such interval u increases from (−a 1 + nT )/c to (−a 1 + (n + 1)T )/c while t increases by T . I − is the union of the intervals [2i + 1, 2i + 2] and on each such interval u decreases from (−a 1 + (n + 1)T )/c to (−a 1 + nT )/c while t is fixed. Mallet-Paret and Nussbaum have considered this Γ (but not our parametrisation of it) extensively, and named it the "sawtooth profile" for the shape of Γ in Fig. 1.1(ii) [30,31].
The motivation for Definition 2.1 comes from numerical simulations, where we observe whenu(t) is finite that α(t, u(t)) is (almost) constant. The sawtooth profile can then be constructed for (1.3) by assuming that α(t, u(t)) is constant with α(t, u(t)) ∈ t(I − ) whenu(t) is finite (that is u(t) ∈ Γ * ). If the phase of the periodic solution is chosen so that t(I − ) = { jT : j ∈ Z} then for t ∈ (0, T ) we have −nT = α(t, u(t)) = t − a 1 − cu(t) for some n ∈ N 0 . Rearranging this leads to the formula for u in (2.2).
Then µ i (η) is continuous on the real line. It is a simple but tedious algebraic exercise to check that (1.5) holds for all η ∈ R. Notice in particular that for η ∈ [3 j, as required. Before checking the conditions on F(η), notice that µ 0 (η) ∈ I * for η ∈ (3 j, and hence F(η) < 0 for all η ∈ (3 j + 1, 3 j + 2], since K > 1. Finally on the interval [3 j + 2, 3 j + 3], we have u(µ 1 (η)) is a linear function of η, while u(µ 0 (η)) is constant, and hence F(η) is a linear function of η. By continuity and the previous calculations F(3 j + 2) = T (1 − K)/c < 0 and F(3 j + 3) = 0 hence F(η) < 0 for all η ∈ J − as required. Thus for each n ∈ N 0 and each K > 1 we have constructed a periodic singular solution of (1.3) defined by (2.4)-(2.5). The parametrisation leading to one of these solutions and the corresponding periodic singular solution is illustrated in Fig. 2.1. Fig. 2.2 shows the periods of the resulting branches as K varies. Using max-plus equations, in [30] this Γ is proved to be the limiting profile of the slowly oscillating periodic solutions (corresponding to n = 0) of (1.3) as ε → 0. In [31] higher order asymptotics reveal the shape of the periodic solution for 0 < ε 1. It is noted that that the asymptotic forms of the periodic solution are very different near to the local maximum and minimum of the solution, with the maximum corresponding to a regular point of the dynamics scaled by ε, while the minimum can be interpreted in the spirit of Fenichel as a turning point near a normally hyperbolic invariant manifold for an ordinary differential equation with a time scaling of ε 2 [31]. The singular solution (2.2)-(2.5) also reveals a difference between the dynamics near to the maximum and minimum of the periodic solution. The solution u(t(µ 0 (η))) has its maximum when η = 3 j + 1 (for any integer j), which is at the boundary between two of the linear segments in the solution parametrisation (2.5), corresponding to the boundary between J * and J − . In contrast u(t(µ 0 (η))) takes its minimum value on the entire interval η ∈ [3 j + 2, 3 j + 3] (but u(t(µ 1 (η))) is not constant on this interval). Note that while at first sight it may have appeared more natural in Definition 1.2 to define J * to be the set of η ∈ J such that µ 0 (η) ∈ I * (or equivalently such that u(µ 0 (η)) ∈ Γ * ), such a definition would be problematical in the example above because µ 0 (η) is constant on ∂ I * on the interval η ∈ [3 j + 2, 3 j + 3]. We will also find nontrivial intervals on which µ 0 (η) is constant on ∂ I * for singular solutions of (1.1). Now consider the two delay DDE (1.1). We assume several conditions on the positive parameters. Without loss of generality we assume that a 2 > a 1 (if not we can either swap the order of the terms, or reduce to an equation with one delay). Then letting α i (t, u(t)) = t − a i − cu(t) we see that α 2 (t, u(t)) < α 1 (t, u(t)) with α 1 (t, u(t)) − α 2 (t, u(t)) = a 2 − a 1 > 0, constant. So although the arguments α i (t, u(t)) are both linearly state-dependent, the difference between them is constant. The more general case where α i (t, u(t)) = t − a i − c i u(t) with c 1 = c 2 so that the difference between the delays is nonconstant would also be interesting, but in the current work we concentrate on understanding the simpler case, which already leads to very complicated dynamics.
It is useful to define the ratio A = a 2 /a 1 > 1 which will play an important role later. If K 1 + K 2 < 1 the trivial solution is asymptotically stable and there are no stable periodic solutions, so we assume that K 1 + K 2 > 1. Finally with a 2 > a 1 we assume that It is shown in [19] that this condition ensures that the DDE initial value problem is well-posed for (1.1), and in particular that the delay α 1 (t, u(t))) < t and so does not become advanced. It is also shown in [19] that when ε > 0 the function α i (t, u(t))) is a strictly monotonic increasing function of t for t a 2 + a 1 (K 1 + K 2 ). Hence α i (t, u(t))) must be a strictly monotonic increasing function of t on any periodic solution. Thus we will construct singular periodic solutions for which all the µ i (η) are monotonic increasing functions of η for all i, although Definition 1.2 only requires that µ 0 (η) be monotonic in general. We first construct singular periodic solutions for (1.1) which have the same sawtooth profile (2.2),(2.3) as the sawtooth solutions of the one delay DDE (1.3). Since these solutions have one local maxima per period we refer to them as unimodal. We will then construct two types singular periodic solution with two local maxima per step; type I and type II bimodal solutions. Each of the solutions that we construct of each type will be characterised by a pair (n, m) of non-negative integers which will have the same meaning in each case. The first number n is the integer number of periods in the past that the first delay falls, and the second number m is the integer number of periods between the two delay times α 1 (t, u(t))) and α 2 (t, u(t))). So for a singular solution of period T we always have t − α 1 (t, u(t)) ∈ [nT, (n + 1)T ], a 2 − a 1 = α 1 (t, u(t)) − α 2 (t, u(t)) ∈ (mT, (m + 1)T ). (2.7) Or using the parametrisation and t(µ 1 (η)) − t(µ 2 (η)) = α 1 (t(µ 0 (η)), u(µ 0 (η))) − α 2 (t(µ 0 (η)), u(µ 0 (η))) With n and m defined by (2.8) and (2.9) to construct unimodal singular solutions of (1.1) it is useful to define θ ∈ (0, 1) by t(µ 1 (η)) − t(µ 2 (η)) = a 2 − a 1 = (m + θ )T, θ ∈ (0, 1), so θ is the fractional part of a period between the two delays, which is assumed to be non-zero. (Although n and m will always have the same meaning, θ will be defined slightly differently for each type of bimodal solution.) As in the one delay case we will construct a solution with t(µ 1 (η)) = −nT while t(µ 0 (η)) ∈ (0, T ). The following theorem establishes conditions for such a solution to exist.
We see immediately from Theorem 2.2 that θ is bounded away from zero. We will see in Section 3 that only certain pairs of values of m, n ∈ N 0 will satisfy the bounds (2.12) in Theorem 2.2. In Theorem 3.1 we will determine which pairs (n, m) are possible and for which parameter ranges the conditions of Theorem 2.2 are satisfied to begin to construct a bifurcation diagram of solution branches. For now, we note that using (2.10) and (2.11) we can write where A = a 2 /a 1 . Using this, the condition θ > K 2 /(K 1 + K 2 − 1) can be rewritten as where When the parameters are such that the bounds on θ in (2.12) are violated other types of singular solution arise. We will construct two such classes of solutions which we refer to as type I and type II bimodal solutions, since each has two local maxima per period.

Fig. 2.4
A generic Type II admissible bimodal periodic solution profile as defined in (2.17). Also shown in colour are the ten stages of the parametrisation of (t(µ i (η), u(µ i (η)) from the proof of Theorem 2.4 for i = 0, 1, 2 and j = 0.
We see from the (2.16) and (2.17) that both solutions have global minima with u = (−a 1 + nT )/c. If the phase of the periodic solution is chosen so that these minima occur when t = jT , for integer j, then for type I bimodal solutions the first local maximum which occurs when t = jT + T 1 is also the global maximum, while for type II bimodal solutions the second local maximum on the period is equal to the global maximum.

21)
then (1.1) has a Type I bimodal singular solution of period T > 0 with solution profile given by (2.16).
and T 1 = T − T 2 . Let the parameters be chosen so that T 2 > 0 and θ ∈ (0, 1) where Then (1.1) has a Type II bimodal singular solution of period T > 0 with solution profile given by (2.17).
Proof This proof is similar to the proof of Theorem 2.3, differing only in the details and conditions, due to differences in the solution profiles and parameterisations. First note that (2.22) is also valid for this solution, while (2.27) can be rewritten as where .
For type I bimodal solutions to exist we require K 1 > 1 in Theorem 2.3. This condition is used twice in an essential way in the proof of that theorem, to show that s 13 > 0 and also that F(10 j + 3) < 0, and so Type I bimodal solutions can only exist for K 1 > 1. In contrast, Theorem 2.4 does not require K 1 > 1, and we will see examples later of type II bimodal solutions which exist for K 1 < 1.
The type I and type II bimodal solutions were constructed so that when α 1 = −nT the second delay α 2 falls in the first (type I) or second (type II) leg of the periodic solution, and for both solutions when the second delay satisfies α 2 = −(n+m)T the first delay satisfies α 1 ∈ (−nT, −nT + T 1 ) so so falls in the first leg of the solution. We also investigated solutions where the first delay satisfies α 1 ∈ (−nT + T 1 , −(n − 1)T ) when α 2 = −(n + m)T and so α 1 falls in the second leg of the solution. However, we did not find examples of such solutions on the branches, so will not present them here.

Bifurcation Branches
Theorems 2.2, 2.3 and 2.4 specify conditions on the parameters for unimodal and Type I and Type II bimodal singular solutions to exist for (1.1). In this section we will show that we can use those theorems to construct bifurcation branches. We require K 2 < 1 to ensure that (1.1) is well posed, while K 1 can be arbitrary large. Thus, it is natural to consider K 1 as a bifurcation parameter.
The unimodal and type I and type II bimodal solutions will be characterized by a pair of integers (n, m) as in the last section, where n and m are related to the delays via (2.7)-(2.9). We will see that each value of n defines a different branch of solutions, with each branch mainly made up of segments of unimodal and type I and II bimodal singular solutions for certain values of m. An example is shown in Fig. 3.1. To explain this example we need to study the parameter conditions from the three aforementioned theorems more closely.
First consider the bounds (2.12) on θ from Theorem 2.2 for the existence of unimodal solutions. By (2.11), the bound θ < 1 is equivalent to a 2 − a 1 < (m + 1)T . Using (2.10) with A = a 2 /a 1 this becomes and hence θ < 1 is equivalent to We already showed that the bound θ > K 2 /(K 1 + K 2 − 1) can be written as G nm (K 1 ) < 0 where G nm (K 1 ) is defined by (2.15). Notice that both bounds only depend on n and m through the common term m − n(A − 1). Let us consider the possible values of m and K 1 > 1 that satisfy these inequalities for fixed values of the other parameters. First define and hence the bound θ < 1 is satisfied for all Finally there is no unimodal solution satisfying the conditions of Theorem 2.2 if m n(A − 1) − 1, since then it is impossible to satisfy (3.1). Now to establish an interval of K 1 parameters on which a unimodal solution exists, we need to consider both the bounds θ < 1 and θ > K 2 /( In the following theorem we establish that for m = m 0 (n) there is a unimodal solution satisfying the conditions of Theorem 2.2 for all K 1 sufficiently large, while for each integer m ∈ (m 0 (n), m * * (n)), there is a non-empty bounded interval of values of K 1 for which (1.1) has a unimodal solution.
defined by (2.15) has distinct roots denote them as M − nm < M + nm , and let M nm be the root of G nm (K) when it has a unique root.  Proof First consider the case when m = m 0 (n) ∈ (n(A − 1) − 1, n(A − 1)). If m < m * (n) (which is always the case if A 3 + K 2 ) then we have θ < 1 for K 1 > L nm > 1. Now consider the polynomial G nm (K 1 ). In this case the coefficient of the quadratic term is negative, and it is easy to verify that G nm (0) > 0 > G nm (L nm ) and hence M + nm < L nm and G nm (K 1 ) < 0 for all K 1 L nm . It follows that (2.12) is satisfied for all K 1 > L nm > 1. On the other hand if m m * (n) then θ < 1 is satisfied for all K 1 > 1 L nm while the coefficient of the quadratic term of G nm (K 1 ) is still negative, but now G nm (1) > 0. In this case G nm (K 1 ) = 0 has a unique positive root K 1 = M + nm > 1 and (2.12) is satisfied for all K 1 > M + nm > 1. Next consider the case when m = m 0 (n) = n(A − 1), which can only arise when A is rational or when n = 0. In this case the quadratic term in G nm (K 1 ) vanishes and the condition G nm (K 1 ) < 0 becomes which can only be satisfied for In that case (3.5) is satisfied for If we also set m = n(A − 1) in (3.3) we obtain Now there are three cases. If A ∈ (1 + K 2 1+K 2 , 3 + K 2 ) then M 00 > 1 and by (3.2) we have m > m * (n) and hence (2.12) is satisfied for all K 1 > M 00 > 1. If A = 3 + K 2 then L 00 = M 00 = 1 and (2.12) is satisfied for all K 1 > 1. Finally if A > 3 + K 2 we have L 00 > 1 > M 00 and (2.12) is satisfied for all K 1 > L 00 > 1. This completes the proof of (i).
To prove (ii), first note that if A 3+K 2 then m * (n) n(A−1) < m 0 (n)+1 and so there is no integer m ∈ (m 0 (n), m * (n)) and nothing to prove. If A > 3 + K 2 then m * (n) > n(A − 1) and for m ∈ (n(A − 1), m * (n)) the bound θ < 1 is satisfied for all K 1 > L nm > 1. Moreover we find that G nm (L nm ) < 0, while the coefficient of the quadratic term is positive so (2.12) is satisfied for all K 1 ∈ (L nm , M + nm ), where M + nm is the largest root of G nm (K 1 ) = 0. In cases (iii) and (iv) we have m > n(A − 1) and m m * (n). The bound θ < 1 is satisfied for all K 1 > 1, while implies that .
In Theorem 3.1(i) we have shown that for m = m 0 (n), the smallest value of m for which a unimodal solution exists, the resulting solution exists for all K 1 sufficiently large. This holds for each integer n 0 and hence, as illustrated in Fig. 3.2(i), we have found the far end of infinitely many solution branches. We note from (2.10) that the period T increases linearly with K 1 on the first (n = 0) branch, but that for n > 0 we have lim K 1 →∞ T = a 1 /n.
The remainder of this work is devoted to the extension and study of these bifurcation branches as well as their persistence for ε > 0. Most of the rest of each solution branch will be composed of legs of other unimodal solutions (with m > n(A − 1)) and of bimodal solutions. Theorem 3.1(ii)-(iv) identifies the parts of the solution branch which are composed of unimodal solutions. This is illustrated in Fig. 3.2(ii).
From (2.10) we see that the unimodal solutions with the largest period occur on the branch n = 0. Let us consider this branch further. From Theorem 3.1(i) provided A = a 2 /a 1 > 1+ K 2 1+K 2 there is a leg of unimodal solutions for n = m = m 0 (0) = 0 for all K 1 > max{L 00 , M 00 }. In that case Theorem 3.1 also ensures there will be legs of unimodal periodic solutions for each integer m between 0 and m * * (0). Hence we require m * * (0) > 1 for there be a second leg of unimodal solutions with n = 0 and m = 1. Fig. 3.3 shows the dependence of m * * (0) on A and K 2 , from which we see that we require the ratio A = a 2 /a 1 > 3 for there to be a second, m = 1, leg of unimodal solutions for K 2 sufficiently small, while for A > 5 there is an m = 1 leg of unimodal solutions for all K 2 ∈ (0, 1). Arbitrary large values of m * * (0) are possible but require A 1. We will explore the case A 1 in Section 4. For other branches of solutions with n > 0, note that m 0 (n) ∈ (m 0 (0)+n(A−1)−1, m 0 (0)+n(A−1)] and from (3.4) we have m * * (n) = m * * (0) + n(A − 1), and so for fixed A and K 2 essentially the same number of legs of unimodal solutions appear for each value of n, but the corresponding values of m are shifted by n (A − 1).
To show for a given value of n that the legs for different values of m form part of a connected branch of solutions we need to join up the branches, which we will do using bimodal periodic solutions of type I and II, and multi-modal solutions.  First we note that if there is a continuous branch for fixed n then for A sufficiently large it must have fold bifurcations. To see this note that if m > m 0 (n) then the coefficient of the quadratic term in (2.15) is positive, while by (3.3) .
But since, M + nm is the largest zero of G nm this implies that L nm−1 < M + nm . By (3.3) we also have L nm < L nm−1 . Thus in the case of legs for m = m 0 (n), . . . , M with M < m * (n) + 1+K 2 2+K 2 we have that L nm−1 ∈ (L nm , M + nm ) for m = m 0 (n) + 1, . . . , M. This implies that the K 1 values of legs for adjacent values of m overlap, just as illustrated in Fig. 3.2(ii). Hence if the legs form part of a continuous branch, the branch must have folds, as seen in Fig. 3.1.  As Fig. 3.1 suggests, these will not be smooth fold bifurcations in the classical sense, but we will see in Section 4 that they are the ε → 0 singular limit of smooth fold bifurcations of periodic orbits, and so we will refer to them as fold bifurcations anyway. We will see that these singular fold bifurcations typically occur at K 1 = M + nm and K 1 = L nm . To show this we need the following lemma which determines the sign of the denominator of (2.18), (2.19) and (2.25).
Proof First note that from (3.3) is a nonincreasing function of K 1 , then it follows that s(K 1 ) < 0 for all K 1 L nm−1 . But this is trivially true if m n, which is always the case when n = 0, since m ∈ N 0 . For n > 0, provided A 2 − 1/n, we have m m 0 (n) + 1 > n(A − 1) n − 1 and hence m n. Finally m > m * (n) + 1+K 2 The following theorem establishes the existence of a fold bifurcation of periodic singular solutions at K 1 = M + nm . As noted before Lemma 3.2, this will not be a smooth bifurcation, but rather a leg of unimodal solutions and a leg of type I bimodal solutions will both exist for K 1 ∈ (M + nm − δ , M + nm ) and these solutions will coincide in the limit as K 1 → M + nm . By coincide, we mean that the limiting profiles and periods of both solutions will be identical. Proof Theorem 3.1 gives the existence of a leg of unimodal solutions for using (2.11),(2.12) and the definition of M + nm . Thus the bounds (2.21) are trivially satisfied when K 1 = M + nm . The bound θ < 1 − 1/K 1 also holds provided K 2 < (K 1 − 1) 2 , in particular whenever M + nm > 2, but M + nm > L nm−1 and m < m * (n) + 1 − Thus all the conditions for the existence of a type I unimodal solution from Theorem 2.3 are satisfied when K 1 = M + nm , and by continuity on an interval containing this point, except possibly for the condition T 2 > 0. But T 2 = a 1 G nm (K 1 )/s(K 1 ) by (2.19). Now, noting that G nm (K 1 ) < 0 for K 1 ∈ (M − nm , M + nm ), and G nm (K 1 ) > 0 for K 1 > M + nm , provided s(M + nm ) = 0, the conditions of Theorem 2.3 must be satisfied on some interval (M + nm − δ , M + nm ) or (M + nm , M + nm + δ ) by continuity of s(K 1 ). But, by Lemma 3.2 we have s(K 1 ) < 0 for all K 1 > L nm−1 , and since L nm−1 < M + nm it follows that for δ > 0 sufficiently small that Thus there is a unimodal solution for K 1 ∈ (L nm , M + nm ) and a bimodal solution on an interval (M + nm − δ , M + nm ) which coincide at a fold bifurcation at For values of m outside the range for which Theorem 3.3 is valid, it can still be possible to obtain type I bimodal and unimodal solutions which coincide at K 1 = M + nm without a fold bifurcation. An example of this will be seen later in Fig. 4.16. We will not determine here the size of δ > 0 such that Theorem 3.3 applies. However, we note that since the theorem guarantees the existence of the unimodal and type I bimodal solutions on some interval, it is a straightforward task to check the conditions of Theorems 2.2 and 2.3 to determine the interval on which each solution exists, and this is what we will do in later examples.
Since the proof of Theorem 3.3 is purely algebraic, it is interesting to consider the bifurcation from a dynamical viewpoint. For the leg of unimodal solutions θ approaches its lower bound in (2.12) as K 1 approaches M + nm . Indeed, since M ± nm are the zeros of G defined by (2.15), it follows that θ → K 2 /(K 1 + K 2 − 1) as K 1 → M ± nm for all of the unimodal solutions identified in Theorem 3.1. At K 1 = M ± nm we have F(5 j + 3) = 0 in the proof of Theorem 2.2. The condition F(5 j + 3) < 0 in that proof ensures that F remains negative while u(µ 2 (η)), the value of u at the second delay, decreases from its maximum value (−a 1 + (n + 1)T )/c to its minimum value (−a 1 + nT )/c. If F(η) = 0 for some η ∈ (5 j + 2, 5 j + 3) then the solution would reenter J * and we would expect another interval on which F(η) = 0. This is exactly what happens in the bifurcation to the type I bimodal solution in Theorem 3.3. For the type I bimodal solution which exists for K 1 < M + nm , from the proof of Theorem 2.3 we see that for η ∈ (10 j + 4, 10 j + 5) the solution at the second delay, u(µ 2 ) decreases from its maximum value (−a 1 +(n +1)T )/c to (−a 1 +nT + (1 −s 14 )(1 −θ )T 1 )/c with F(η) < 0, but F(η) → 0 as η → 10 j +5. For η ∈ (10 j + 5, 10 j + 6) we have η ∈ J * , F(η) = 0 and u(µ 2 ) further decreases to its minimum value (−a 1 + nT )/c. However as K 1 approaches M + nm we have G nm (K 1 ) → 0, and hence T 2 → 0 and s 14 → 1, so (−a 1 + nT + (1 − s 14 )(1 − θ )T 1 )/c is no longer larger than (−a 1 + nT )/c and F does not become zero before u(µ 2 ) reaches its global minimum. Hence the second interval of J * for η ∈ (10 j + 5, 10 j + 6) collapses, and as s 11 and s 12 both tend to 1, we find that five of the intervals of the parametrisation of the type II bimodal solution become trivial, and the remaining parts correspond to the unimodal solution. Thus at the bifurcation between the unimodal solution and the type I bimodal solution θ hits its lower bound for the unimodal solution, and T 2 → 0 for the bimodal solution. It will be interesting to investigate below what other bifurcations arise as other conditions in the theorems of Section 2 are violated. Now consider the case of K 1 = L nm−1 at the left-hand end of the interval of unimodal solutions for K 1 ∈ (L nm−1 , M + nm−1 ). We show that at this point there is a fold bifurcation and the solution transforms from a unimodal solution to a type II bimodal solution. By the definition of L nm for K 1 > L nm−1 the unimodal solution satisfies a 2 − a 1 = (m − 1 + θ )T with θ ∈ (0, 1) but as K 1 → L nm−1 we have θ → 1. But if θ were equal to 1, the difference a 2 − a 1 between the two delayed times would be exactly m periods. Perhaps not surprisingly, as the following theorem shows, this can result in a (type II) bimodal solution with the value of m increased by 1. Proof Theorem 3.1 gives the existence of a leg of unimodal solutions with m = p−1 for K 1 ∈ (L np−1 , +∞) when p = m 0 (n)+ 1, for K 1 ∈ (L np−1 , M + np−1 ) when p ∈ (m 0 (n) + 1, m * (n) + 1], and for K 1 ∈ (M − np−1 , M + np−1 ) when p ∈ [m * (n) + 1, m * * (n) + 1). To prove (i) and (ii), next we show that if there exists a leg of type II bimodal solutions with m = p for K 1 ∈ (L np−1 − δ , L np−1 ) or K 1 ∈ (L np−1 , L np−1 + δ ) then the unimodal and type II bimodal solutions must coincide in the limit as K 1 → L np−1 . To see this, compare the profile of the type II bimodal solution in (2.17) with the profile of the unimodal solution in equations (2.2),(2.3). The two solutions will coincide in the limit as K 1 → L np−1 if both the unimodal and type II bimodal solution have the same limiting period T and for the type II bimodal solution θ → 0 as K 1 → L np−1 . But it is simple to check that the value of T given by (2.18) for the type II bimodal solution with m = p agrees with that given by (2.10) for the unimodal solution with m = p − 1. The rest of this proof concerns the existence and properties of the type II bimodal solution with m = p, so we can use m and p interchangeably. To show that θ → 0 as K 1 → L np−1 for the type II bimodal solution with m = p, note that by (2.18) Now, s(L np−1 ) = 0 by Lemma 3.2, since p = m * (n) + 1+K 2 2+K 2 . Hence as required, provided T 2 > 0. To derive expressions for T 1 and T 2 when K 1 = L nm−1 , using (2.25) and (3.11) Thus T 2 = K 2 T and T 1 = (1 − K 2 )T when K 1 = L nm−1 . Since K 2 ∈ (0, 1) this implies that T 1 > 0 and T 2 > 0 in the limit as K 1 → L np−1 , which establishes that the unimodal solution and type II bimodal solution have the same limiting profiles as To prove (i) and (ii), it remains to verify the conditions of Theorem 2.4 to confirm the existence of the type II bimodal solution. First note that since T 1 > 0 and T 2 > 0 when K 1 = L nm−1 , and s(L nm−1 ) = 0 there exists δ > 0 such that T 1 and T 2 defined by (2.18),(2.25) vary continuously and are strictly positive for K 1 ∈ (L nm−1 − δ , L nm−1 + δ ). Now consider the condition θ > 0. From above Under (i) we have s(L nm−1 ) < 0 and hence s(K 1 ) < 0 for K 1 ∈ (L nm−1 , L nm−1 + δ ) for δ sufficiently small. Also by (3.11) for Similarly under (ii) we have s(L nm−1 ) > 0 and hence s(K 1 ) > 0 for K 1 ∈ (L nm−1 − δ , L nm−1 ) for δ sufficiently small, and by (3.11) we have . Moreover since θ → 0 as K 1 → L nm−1 in both cases, for δ > 0 sufficiently small we also have θ < 1.
To prove (iii) it remains only to establish the existence of the type II bimodal solution in that case, but this is similar to above, where we note that p < (n + 1 2 )(A − 1) implies L np−1 > 1 − K 2 and choosing δ sufficiently small so that L np−1 − δ > 1 − K 2 ensures that K 1 + K 2 > 1 for all K 1 ∈ (L np−1 − δ , L np−1 ). This implies that the second term on the right-hand side of (2.29) is strictly positive, while the first expression tends to (L np−1 + K 2 − 1)/(L np−1 K 2 ) > 0 in the limit as K 1 → L np−1 . Again, since lim K 1 →L np−1 θ = 0, for δ > 0 sufficiently small equation (2.29) is satisfied for K 1 ∈ (L np−1 − δ , L np−1 ). Theorem 3.4(i) establishes the existence of a fold bifurcation when K 1 = L np−1 for p ∈ m 0 (n), m * (n) + 1+K 2 2+K 2 . Interestingly, for p ∈ (m * (n) + 1+K 2 2+K 2 , m * (n) + 1] the fold disappears, but the two legs of periodic solutions continue to exist and coincide at K 1 = L np−1 , but now the type II bimodal solution exists for K 1 < L np−1 while the unimodal solution exists for K 1 > L np−1 . Essentially, the fold bifurcation unfolds suggesting a cusp-like bifurcation of periodic orbits, which we will investigate in Section 4. Theorem 3.4(iii) also indicates interesting behaviour. When m = m * (n) + 1, or equivalently, we have L nm−1 = L nm * = 1. Noting that m * (n) + 1 < (n + 1/2)(A − 1), for m ∈ (m * (n) + 1, (n + 1/2)(A − 1)), or equivalently for we have 1 > L nm−1 > 1 − K 2 and Theorem 3.4(iii) ensures the existence of type II bimodal solutions for In contrast the construction of the unimodal and type I bimodal solutions in Theorems 2.2 and 2.3 requires K 1 > 1 for those solutions to exist. Dynamically, we see in the proof of Theorem 3.4 that for both the unimodal and type II bimodal solution we have T → (a 2 − a 1 )/m as K 1 → L nm−1 . For the unimodal solution (a 2 − a 1 ) = (m − 1 + θ )T < mT and θ → 1 as K 1 → L nm−1 , while for the type II bimodal solution (a 2 − a 1 ) = mT + θ T 2 > mT and θ → 0 as K 1 → L nm−1 . Whether or not there is a (non-smooth) fold bifurcation at K 1 = L nm−1 depends on whether m is greater or smaller than m * (n) + (1 + K 2 )/(2 + K 2 ). That the value of m increases close to K 1 = L nm−1 was already observed numerically for ε > 0 in [19].
Theorem 2.2 identified upper and lower bounds on θ for a unimodal solution to exist. In Theorems 3.3 and 2.4 we have shown bifurcations to type I or type II bimodal solutions when one of these bounds is violated. In Section 4 we will investigate the solutions that can arise when the parameters bounds identified in Theorems 2.3 and 2.4 for type I and type II bimodal solutions are violated.

Singularly Perturbed Solution Branches
We are interested in solutions of (1.2) when 0 < ε 1. However, so far we have only constructed ε = 0 singular solutions, in the sense of Definition 1.2. It would be desirable to prove that (1.2) has solutions close to the constructed singular solutions for all ε sufficiently small. Mallet-Paret and Nussbaum [31] proved that the sawtooth is indeed the limiting profile as ε → 0 for the state-dependent DDE (1.3) which has one delay. However, for the two delay problem (1.1), Theorems 3.3 and 3.4 lead us to expect fold bifurcations of periodic orbits. Indeed such bifurcations and resulting intervals of co-existing stable periodic solutions were already observed for ε = O(1) in [19]. 'Superstability' is central to the results of [31], and without further insight it is difficult to see how to modify the techniques of [31] to rigorously prove the persistence of the singular solutions for ε > 0 given that it is possible for (1.1) to have co-existing stable periodic orbits.
Given the analytical difficulties, in the current work we will instead pursue a numerical approach, and compute the periodic solutions and bifurcation structures for (1.1) numerically for 1 ε > 0. From the numerical solutions we will infer that over wide parameter ranges the singular solutions identified in the theorems of Sections 2 and 3 are indeed the limits of the solutions of (1.1) as ε → 0. Moreover we will find that that (1.1) has fold bifurcations of periodic orbits at K 1 values which converge to K 1 = L nm and K 1 = M nm + in the limit as ε → 0.
We will numerically compute the bifurcation branches using DDEBiftool [11]. This computer package consists of a suite of matlab [32] routines for computing solution branches and bifurcations of DDEs using path following and branch switching techniques. Periodic orbits are found as the solution of a boundary value problem (BVP), using collocation techniques. The numerical analysis details are well described in [11] and elsewhere, so we will not repeat them here. We emphasise however, that periodic orbits are found by solving BVPs, and not by solving initial value problems. This allows us to find unstable orbits just as easily as stable ones. DDEBiftool can also determine the stability of periodic orbits by computing their Floquet multipliers which also allows us to detect bifurcations. In [19] we already used DDEBiftool to investigate the dynamics of (1.1) in the non-singular case ε = 1.
Although DDEBiftool can compute unstable periodic solutions just as easily as stable ones, we will mainly concentrate our attention on the principal branch of periodic solutions, which is the only branch on which we found large amplitude stable periodic solutions for ε > 0. By the principal branch, we mean the branch of periodic orbits which has the largest period among all the Hopf bifurcations, both at the bifurcation and in the limit as K 1 → ∞. This usually also corresponds to the Hopf bifurcation with the smallest value of K 1 , but due to the vagaries of the behaviour of the characteristic values in DDEs for ε very close to zero, it is sometimes possible for a shorter period orbit to bifurcate first. If that happens then the periodic orbit on the principal branch is initially unstable but we found numerically that it becomes stable in a torus bifurcation while its amplitude is still very small. In the current work we will not study small amplitude solutions or invariant tori (see [4] for a study of the invariant tori of (1.2), and [22] for a study of solutions of (1.3) close to the singular Hopf bifurcation). The principal branch will always correspond to the choice n = 0 for the singular solutions and hence m 0 = 0.
We begin our exploration of the periodic solutions of (1.1) in Fig. 4.1 with an example which for ε = 0.1 and ε = 0.05 shows the amplitude and period of the periodic orbits on the principal branch. Here, and throughout this work, the amplitude of a periodic orbit of period T is defined simply as the difference between the maximum and minimum values of u(t) over the period. We will take c = 1 and K 2 = 0.5 in all our numerical computations, and a 1 = 1, so A = a 2 /a 1 = a 2 .
Also shown in Fig. 4 solutions exist their period is very close to that of the numerically computed ε > 0 solutions. The agreement in the amplitude is not quite so good, particularly closer to the Hopf bifurcation, with the numerically computed orbits having slightly smaller amplitudes (which is to be expected since the singular solutions are of sawtooth shape, while for ε > 0 the periodic orbits are smooth, and some amplitude is lost in the "smoothing" of the sawteeth for ε > 0). We see in Fig. 4.1 that the ε > 0 branches pass continuously through the gap in the singular solution branch, and so we next investigate periodic solution profiles as ε → 0 paying particular attention to the legs between the fold bifurcations where those gaps occur. Fig. 4.2 shows the amplitude and stability of periodic solutions on the two branches of periodic solutions bifurcating from the trivial solution at the first two Hopf bifurcations with ε = 0.1 and A = 6. The stability of the orbits is determined from the Floquet multipliers computed in DDEBiftool. There is an interval of bistability of periodic solutions between the two folds on the first branch (first noted in [19]). Periodic solutions on the leg of the first branch between the two folds are always unstable, and are bimodal for at least part of the leg. As the first inset shows, the amplitude of the periodic solutions does not vary monotonically with K 1 on this leg. Similar behaviour is seen on the second branch, except this branch is always unstable (there are two Floquet multipliers greater than one on the leg between the folds and one elsewhere). Fig. 4.3 shows the profiles and periods of the numerically found periodic orbits on the principal branch for the same parameter values as in Fig. 4.2. There are two legs of stable unimodal periodic solutions, and one leg of unstable solutions. We remark that the leg of stable unimodal solutions existing for all K 1 sufficiently large, corresponds to the unimodal  singular solutions existing for all K 1 sufficiently large identified in Theorem 3.1(i). Indeed, for the parameters considered here those unimodal singular solutions exist for all K 1 > L 00 where L 00 = A − 2 − K 2 = 3.5 by (3.3), with period T = a 2 − a 1 = 5 when K 1 = L 00 (using (2.10)), and period T > a 2 − a 1 = 5 for K 1 > L 00 (using (2.11) with m = 0 and θ < 1). The other leg of unimodal solutions starting at the Hopf bifurcation, corresponds to the unimodal singular solution identified in Theorem 3.1(iii) with m = m * = 1 and T < a 2 − a 1 < 2T (ie a 2 − a 1 = mT + θ with m = 1 and θ ∈ (0, 1)).
On the leg of unstable solutions in Fig. 4.3 we see two types of bimodal periodic solutions for ε = 0.1, where either the first or second local maximum after the solution minimum is higher (see insets with K 1 = 4.8001 and K 1 = 4.1993). Note the resemblance between the profiles of the ε = 0.1 solutions and the singular solutions shown in Figs. 2.3 and 2.4, which is not coincidental; the construction of the singular solutions in Section 2 was guided by preliminary numerical computations.
With ε = 0.1 there is a smooth transition between the two types of bimodal solutions along the unstable leg, whereas Fig. 3.1 suggests that there should be a gap between the intervals where these two types of solutions exist in the limit as ε → 0. To investigate the transition between these solutions we decrease ε progressively. Fig. 4.4 shows that the numerically computed solution for ε = 0.05 has periodic solutions with three local maxima per period (which we call trimodal solutions) for an interval of K 1 values which falls within the gap in the K 1 values between the type I and type II bimodal singular  solutions. Fig. 4.5 shows the profiles of two of these trimodal periodic solutions. We note that both profiles are similar to bimodal solutions, but that in both cases the first local maxima of the bimodal solution has split into two local maxima.
For parameter values close to where the type I bimodal solutions exist (including K 1 = 4.659) the first two local maxima of the solution resemble those of a type I bimodal solution (with the first local maxima after the global minima being the global maxima), while for parameter values close to the type II bimodal solutions (including K 1 = 4.6294) the first two local maxima of the solution resemble those of a type II bimodal solution (with the second local maxima after the global minima being the global maxima).
To investigate further we reduce ε to 0.01. This is smallest value of ε that we can achieve for these parameters with DDEBiftool. Indeed the bifurcation branch shown in Fig. 4.6 terminates shortly after passing the region of interest, when the numerical algorithm fails to find a solution to the specified numerical tolerance. DDEBiftool computes periodic solutions as smooth functions using collocation, but in the limit as ε → 0 the solutions are forming vertical segments and so are not smooth, which restricts how small we can take ε in numerical computations.
In Fig. 4.6 with ε = 0.01 Type I-like bimodal solutions are seen on the branch in the approximate range K 1 ∈ (4.6998, 4.9802). At K 1 ≈ 4.6998 there is a transition to a trimodal solution, and trimodal solutions exist in the interval K 1 ∈ (4.6735, 4.6998). The numerically found trimodal solution for K 1 = 4.6908 is illustrated in Fig. 4.7(i). Again we see (in the inset) that it is the first maximum of the solution which splits into two to form the trimodal solution. Around K 1 ≈ 4.673 there is a brief interval of quadrimodal solutions, where the first maximum of the trimodal solution splits into two as illustrated Fig. 4.7(iii). There is then another interval of trimodal solutions for K 1 ∈ (4.5746, 4.673), with the solution for K 1 = 4.6266 shown in Fig. 4.7(ii). Finally for K 1 < 4.5746 the solutions are bimodal (and type II-like). Comparing the trimodal solutions in Fig. 4.7 with those in Fig. 4.5 we see that the trimodality is much more clearly defined for the smaller value of ε with the profiles in Fig. 4.7 much more 'sawtooth-like' than the smoother profiles seen in Fig. 4.5. Moreover, we see that for both ε = 0.05 and ε = 0.01 the trimodal solution in the interval adjacent to the type I bimodal solutions has a larger first peak than second peak, just as the type I bimodal solutions do, and similarly for type II bimodal solutions and the trimodal solutions in the adjacent parameter interval the second peak is larger. This could lead us to define type I and type II trimodal solutions which could then be found algebraically in the ε = 0 limit using the singular solution techniques of Section 2, however each would involve about 15 intervals of parametrisation, which would be tedious beyond belief. Moreover, Fig. 4.7 suggests that the quadrimodal solutions also come in both types, and we suspect that as ε → 0 there is a cascade of solutions with arbitrary many maxima and minima, and some form of self-similarity to the evolution of the periodic solution profile in the limit as ε → 0. However all these multimodal solutions lie on the unstable leg of the bifurcation branch, and we will not pursue their construction here.
Although we omit the algebraic construction of the trimodal solutions, our algebraic construction of the bimodal solutions is sufficient to compare the bimodal to trimodal solution transition with the unimodal to bimodal transitions/bifurcations already studied algebraically in Section 3. The interval of values of the parameter K 1 for which type I bimodal solutions exist was found by checking the conditions of Theorem 2.3. In all the examples shown above, and indeed in the other examples of type I bimodal solutions that occur later in this paper we only find two different behaviours which arise at the ends of these intervals. One case is when K 1 = M + nm indicating a transition or bifurcation between a unimodal and type I bimodal solution as studied in Theorem 3.3 and the comments after the theorem. At the other end of the leg of type I bimodal solutions where the ε > 0 numerics indicates a transition to a trimodal solution, algebraically in all the examples shown here we find that the lower bound on θ T 1 in (2.21) fails. From the proof of Theorem 3.3 we see that equality in this bound corresponds to F(10 j +3) = 0. This situation is similar to the transition from a unimodal to type I bimodal solution at K 1 = M + nm as described just after Theorem 3.3. In that case we saw that the failure of the condition F(5 j + 3) < 0 in the unimodal singular solution led to the creation of a second subinterval of J * in the periodic orbit. In an analogous manner the failure of the condition F(10 j + 3) < 0 in the type I bimodal singular solution can lead to a solution where J * consists of three disjoint intervals per period and the resulting solution is trimodal.
The transition from type II bimodal to trimodal solutions also appears to be similar to the transition from unimodal to type II bimodal solutions. After Theorem 3.4 we noted that at the transition between unimodal and type II bimodal solutions at K 1 = L nm−1 we have θ → 1 for the unimodal solution and θ → 0 for the type II bimodal solution as K 1 → L nm−1 . Checking the conditions of Theorem 2.4 we find that in all the examples above we have θ → 0 as K 1 → L nm−1 and θ → 1 at the other end of each leg of type II bimodal solutions. We would expect the solution to transition to a type II trimodal solution with θ = 0 at this point.

Cusp-like Bifurcations
As noted after Theorem 3.4, that theorem establishes the existence of a fold bifurcation of singular solutions for m < m * (n) + 1+K 2 2+K 2 which disappears in a cusp-like bifurcation when m = m * (n) + 1+K 2 2+K 2 and K 1 = L nm−1 , or equivalently when For the principal branch n = 0, so the cusps occur when A = 2 + m(2 + K 2 ) and L 0m−1 = 1 + 1/m. Taking K 2 = 0.5, the cusp-like bifurcations for the singular solutions occur when A = 4.5, 7, 9.5, ... and L 0m−1 = 2, 3/2, 4/3, . . .. Here we will investigate the first two cusp bifurcations numerically for 0 < ε 1. Figs. 4.8-4.11 illustrate the change in the dynamics near to K 1 = L 00 as A passes through 4.5, showing amplitude and period plots of the principal branch of periodic solutions computed with ε = 0.02 for A = 4.48 and A = 4.54 and some numerically computed solution profiles shown as insets. For comparison the legs of unimodal and bimodal singular solutions are also plotted. By Theorem 3.4 we have a leg of unimodal singular solutions for K 1 > L 00 = A − 2 − K 2 with period T → A − 1 as K 1 → L 00 , and since c = 1, amplitude equal to the period. For A = 4.48 and K 2 = 0.5 this gives L 00 = 1.98 and period and amplitude of unimodal solutions tending to 3.48 as K 1 → L 00 . Moreover (3.2) implies m * (0) + 1+K 1 1+K 2 < 1 and Theorem 3.4(ii) gives the existence of a leg of type II bimodal solutions for K 1 < L 00 . We see from Figs. 4.8 and 4.9 that the solution branch for ε = 0.02 behaves similarly though the transition point is perturbed to K 1 ≈ 1.9221 which is slightly less than L 00 = 1.98.
For A = 4.54, m 0 (0) = 0 and m * (0) + 1+K 1 1+K 2 > 1 and Theorem 3.4(i) gives the existence of a leg of type II bimodal solutions for K 1 > L 00 , resulting in a fold bifurcation. Figs. 4.10 and 4.11 show that the ε = 0.02 solution branch also has a    One important aspect of this cusp-like bifurcation is that it has the potential to create stable bimodal and multimodal periodic solutions. In Figs. 4.1-4.7 all of the solutions with more than one local maxima per period were unstable occurring on the unstable leg of the bifurcation branch between the fold bifurcations. The bimodal and trimodal solutions occurring       a bimodal to a unimodal solution close to K 1 = L 01 without a fold bifurcation, while for A = 7.08 the same transition is associated with a fold bifurcation. As was the case with the first cusp-like bifurcation for ε = 0.02 this transition occurs for a value of K 1 slightly less than L 01 both when A = 6.98 and A = 7.08. Also, whereas the fold appears when A   Fig. 4.13 shows that the ε = 0.02 branch has solutions whose K 1 values and periods almost exactly agree with those of the unimodal m = 2 singular solutions while the inset for K 1 = 1.601 in Fig. 4.12 shows that while the ε = 0.02 solution has smaller amplitude and is bimodal, its profile is close to unimodal. The inset for K 1 = 1.2962 shows a trimodal solution for ε = 0.02 occurring before the pair of fold bifurcations. Such a solution is also seen in Fig. 4.11, and they seem to be ubiquitous, also arising even when the folds disappear (see inset for K 1 = 1.5997 in Fig. 4.9).

Other Solutions and Bifurcations
So far we have concentrated our attention on the fold and cusp-like bifurcations predicted by Theorems 3.3 and 3.4, but the folds occur between legs of unimodal solutions, and as noted in the discussion after Theorem 3.1, we require A > 3 to have more than one leg of unimodal solutions on the principal n = 0 branch of periodic solutions. We now consider the dynamics when A ∈ (1, 3). Noting that Theorem 3.1 guarantees the existence of unimodal solutions for all K 1 sufficiently large unless m = m 0 (n) = n(A − 1) and A 1 + K 2 /(1 + K 2 ) we first consider this exceptional case. On the principal branch this occurs when m = n = 0 and taking K 2 = 0.5 with A ∈ [1, 4/3]. Consequently in Fig. 4.16 we consider the dynamics with A = 7/6.
With A = 7/6 and ε = 0.05 DDEBiftool reveals a Hopf bifurcation at K 1 ≈ 0.5373 leading to a branch of stable solutions which exists for all larger values of K 1 . Close to the Hopf bifurcation these periodic solutions are unimodal and sinusoidal, but for all K 1 > 1.5167 these solutions have two local maxima per period, and closely resemble type II bimodal solutions (see K 1 = 1.9817 and K 1 = 6.8797 insets in Fig. 4.16). Indeed, verifying the conditions of Theorem 2.4 we find that there is a type II bimodal solution with n = m = 0 when A = 7/6 for all K 1 > 0.59816. This is also shown in Fig. 4.16 and has very good agreement with the amplitude of the numerically found ε = 0.05 solutions when K 1 > 1. Type II bimodal singular solutions and their ε > 0 counterparts are also found for all K 1 sufficiently large for other values of A ∈ (1, 1 + K 2 /(1 + K 2 )) when m = m 0 (n) = n(A − 1). Fig. 4.16 also shows the n = 1 branch for ε = 0.05 when A = 7/6. DDEBiftool finds the solutions on this branch to all be unstable. By Theorem 3.1(i) there is a unimodal singular solution with n = 1, m = 0 for all K 1 > M + 10 = 2.7625, while verifying the conditions of Theorem 2.3 reveals that there is a type I bimodal solution for K 1 ∈ (2.0481, M + 10 ). At K 1 = M + 10 the two solutions coincide, with T 2 → 0 as K 1 → M + 10 for the type I bimodal solution. Theorem 3.3 deals with unimodal and type I bimodal solutions coinciding at K 1 = M + nm in a fold-like bifurcation. That theorem does not apply here because we have m = 0 < n(A − 1) = 1/6 outside its range of validity, nevertheless we still have the transition between the two types of solution, but here it occurs without the fold-like bifurcation.    Next we consider A ∈ (4/3, 3) for which Theorem 3.1 guarantees the existence of a unimodal singular solution with n = m = 0 on the principal branch for all K 1 sufficiently large, but for which with n = 0 there is no value of m which satisfies the conditions of Theorem 3.3 or 3.4. Taking A = 1.5, Theorem 3.1(i) gives the existence of unimodal singular solutions with n = m = 0 for K 1 > M 00 = 5. Similarly to the n = 1 branch of the previous example, verifying the conditions of Theorem 2.3 we find a type I bimodal singular solution with n = m = 0 for K 1 ∈ (3.3508, M 00 ). For ε = 0.05, using DDEBiftool we numerically compute the principal solution branch from the Hopf bifurcation (at K 1 ≈ 0.7363), finding stable bimodal periodic solutions for K 1 ∈ (3.3414, 5.0543), and stable unimodal solutions for K 1 > 5.0543 as shown in Fig. 4.17. The parameter ranges and amplitudes of the solutions with ε = 0.05 are seen to be very close to the ε = 0 singular solutions. Solutions are also stable on the initial part of the branch and unimodal for K 1 ∈ (0.7363, 0.8586) and bimodal for K 1 ∈ (0.8586, 2.3347). However, the solutions on principal branch are unstable in the range K 1 ∈ (D 1 , D 2 ) = (2.3347, 3.3414) with DDEBiftool detecting period doubling bifurcations (characterized by a real Floquet multiplier passing through −1) at both ends of this interval. On the principal branch trimodal solutions are found for K 1 ∈ (2.9016, 3.3112) while the solutions are bimodal in the rest of the interval (D 1 , D 2 ). In Fig. 4.18 we show the resulting branch of stable period-doubled solutions for K 1 ∈ (D 1 , D 2 ). The branch is computed from K 1 = D 1 and appears to terminate at K 1 = D 2 , though numerical computation of the branch is very difficult near to K 1 = D 2 . Insets show profiles of the resulting stable periodic solutions which all have period close to 7 and mainly have four local maxima per period, except for K 1 ∈ (2.5247, 2.6909) where the first peak splits into two (reminiscent of the earlier transitions from bimodal to trimodal solutions) resulting in periodic solutions with five local maxima per period.
We do not have a characterization from the singular solutions of when to expect period doubling bifurcations with ε > 0. To determine the parameter ranges where period doubled orbits can appear we could parameterise the period doubled singular periodic orbits. This would be similar to a perturbation of two copies of the parameterisations illustrated in Figs. 2.3 and 2.4 and would involve twenty parametrisation intervals and more algebraic manipulation than we care to contemplate. We note that at the end of the interval of validity of the type I bimodal solution at K 1 = 3.3508 the left inequality in (2.21) is tight, and fails for smaller values of K 1 . This is the same inequality that failed at the transition between type I bimodal and trimodal solutions between the fold bifurcations at K 1 = L nm−1 and K 1 = M + nm which we studied earlier. Given the proximity of the end of the interval of type I bimodal singular solution at K 1 = 3.3508 to the period doubling bifurcation with ε = 0.05 at K 1 = 3.3414 and the transition from bimodal to trimodal solutions at K 1 = 3.3112 we suspect that in the limit as ε → 0 the singular solutions undergo a period doubling bifurcation at the same parameter value where the periodic solution transitions from type I bimodal to trimodal. The behaviour seen in this example contrasts with earlier examples where no period doubling bifurcations were detected between the fold bifurcations.
We have already seen examples corresponding to Theorem 3.4(i) and (ii) with unimodal and type II bimodal solutions which coincide at K 1 = L nm−1 with or without a fold bifurcation. Fig. 4.19 illustrates Theorem 3.4(iii), showing a type II bimodal singular solution which exists for K 1 < L nm−1 and a unimodal solution for When n = 0 have L nm * = 1 at A = 1 + m(2 + K 2 ) so A = 3 + K 2 , 5 + 2K 2 , 7 + 3K 2 , . . . and so separated unimodal and type II bimodal solutions will occur for A slightly smaller than these values. In Fig. 4.19 we consider K 2 = 0.5 and A = 5.75 < 6 = 5 + 2K 2 .
Theorems 2.2 and 2.3 require K 1 > 1 for the construction of unimodal and type I bimodal solutions, but Theorem 2.4 only requires K 1 > 1 − K 2 for the construction type II bimodal solutions, and in Fig. 4.19 we have an example of type II bimodal solutions which exist for K 1 < 1. Fig. 4.19 also shows a numerically computed branch of periodic orbits with ε = 0.005 which passes very close to the legs of bimodal type I and unimodal singular solutions. While the type II bimodal solution exists for K 1 ∈ (0.6621, 0.875) the ε = 0.005 branch has bimodal solutions for K 1 ∈ (0.6444, 0.7822) with the inset solution profile for K 1 = 0.7522 showing that these resemble type II singular solutions. The ε = 0.005 branch also has unimodal solutions for K 1 > 1.0126 which approximate the unimodal singular solutions existing for K 1 > 1.0348. It was found that the period of the solutions on the ε = 0.005 increases monotonically from T ≈ 2.2478 at the Hopf bifurcation and crosses T = 2.375 at K 1 ≈ 0.84265. At this value of K 1 the period T satisfies 2T = 4.75 = a 2 − a 1 , that is the difference between the delays is exactly two periods. For the singular solutions 2T = a 2 − a 1 when K 1 = L 01 = 0.875 at the end of the interval of bimodal type II solutions. Fig. 4.19 suggests that there is not a bifurcation near to K 1 = L nm−1 < 1 when the conditions of Theorem 3.4(iii) are satisfied, but neither do the solutions transition directly from type II bimodal to unimodal solutions, as occurs in Theorem 3.4(i) and (ii). We do not have an explanation for the dynamics for K 1 ∈ (L 01 , M − 01 ) in Fig. 4.19, but note with ε = 0.005 bimodal solutions are seen for K 1 ∈ (0.9606, 1.0066) and multimodal solutions for K 1 ∈ (1.0066, 1.0126). The multimodal solution for K 1 = 1.0115 is illustrated in Fig. 4.19. At K 1 = 1.0126 the numerical solution transitions directly from multimodal to unimodal at the kink in the bifurcation branch.
Another example where bimodal type II singular solutions could exist for K 1 < 1 was already seen in Fig. 4.4. Fig. 4.4 illustrated the boundary between Theorem 3.4(ii) and (iii) with m = m * (n) and L nm = M − nm = 1. We have concentrated our attention on unimodal and bimodal solutions in this work, but noting that trimodal and even quadrimodal solutions arise between legs of type I and type II bimodal solutions. Fig. 4.20 shows examples of multimodal solutions with up to seven local minima per period (see the K 1 = 1.372 inset). The parameters in Fig. 4.20 are the same as those considered in Figs. 4.8-4.9, when we were studying the cusp-like bifurcation occurring at K 1 = L 00 = 2 with A = 4.5. Fig. 4.20 shows that even for A < 4.5 when there are no fold bifurcations near to K 1 = L 00 , there are still six fold bifurcations earlier on the principal branch, and there are multimodal solution profiles near to each of these folds. The multimodal solution profiles shown in the figure for K 1 ∈ (1.3, 1.44) all have well-defined 'sawteeth' with the periodic solution profile having gradient close to 1/c = 1 before each local maxima and large negative gradient afterwards. It seems plausible that the fold bifurcations associated with the transitions between unimodal and bimodal solution profiles that we studied earlier are just the simplest example of a sequence of such bifurcations that occur at points where the number of local maxima in the periodic solution profile changes. In principle Definition 1.1 and the techniques of this paper could be used to locate these bifurcations in the ε → 0 limit. Fig. 4.21 illustrates the convergence of the solution branches as ε → 0. We observe that the branch converges over its entire length as ε → 0, including the Hopf points and the the smooth small amplitude solutions close to the Hopf points that we have not studied in this work. Insets show solution profiles on the unstable middle leg of solutions for K 1 = 4.4, K 2 = 0.5, converging to a type II bimodal solution as ε → 0. Even for ε as large as 0.2 the bimodal sawtooth structure of the solution profile is very clearly seen. For larger ε the solution profiles are smoother, particularly near the local maxima, but the fold structure on the solution branch persists even when ε = 1.
So far we have mainly considered the principal branch of periodic solutions corresponding to singular solutions with n = 0, since this is the branch on which stable solutions can be observed. However, it was demonstrated in [19] that there are infinitely many Hopf bifurcations for ε > 0 and we finish this work by considering the alignment of the bifurcations on the different branches. This is illustrated in Figs. 4.22 and 4.23 for ε = 0.05. With A = 6 we saw earlier that in the limit as ε → 0 the fold bifurcations occur on the principal n = 0 branch at K 1 = L 00 = 3.5 and K 1 = M + 01 = 5.  To explain this alignment notice that in the singular limit ε → 0 by Theorems 3.3 and 3.4(i) for suitable integer value(s) of m there are fold bifurcations on the n-th branch at K 1 = L nm and K 1 = M + nm−1 . But L nm defined by (3.3) and M + nm defined as the larger zero of G nm (K 1 ) (see (2.15)) both depend on n and m only through the common term m − n(A − 1). Hence if A = p/q is rational then defining n k = n 0 + kq, m k = m 0 + k(p − q), k ∈ N we see that m k − n k (A − 1) = m 0 − n 0 (A − 1), and hence L n k m k = L n 0 m 0 and M + n k m k = M + n 0 m 0 for each integer k and for each n 0 = 0, 1, 2, . . . , q − 1. Hence these singular fold bifurcations align on every q-th branch when A = p/q is rational. Thus, when A is integer these bifurcations align on all the branches (eg A = 6, see Fig. 4.22(i)), when A = p/2 the bifurcations align on every second branch (eg A = 5.5, see Fig. 4.23), and when A is irrational there is no alignment between the bifurcations (see Fig. 4.22(ii)).
Moving to the ε > 0 case, we see from the figures with ε = 0.05 that the fold bifurcations which should align exactly in the limit as ε = 0, actually appear to occur within shrinking subintervals of [L n 0 m 0 , M + n 0 m 0 −1 ] as n k is increased, and for sufficiently large n k the fold bifurcations disappear entirely. Although for each fixed n k the folds occur for all ε sufficiently small and converge to K 1 = L n 0 m 0 and K 2 = M + n 0 m 0 −1 as ε → 0 the convergence is clearly not uniform, with smaller values of ε required to create the fold bifurcations for larger values of n k . This is not surprising, since the larger the value of n k the smaller the period and amplitude of the ε = 0 singular solutions defined by Theorems 2.2, 2.3 and 2.4, but when solving with ε > 0 the smaller amplitude solutions appear smoother and more sinusoidal than the larger amplitude solutions and the fold bifurcations do not occur unless ε is reduced sufficiently to resolve the sawteeth in the solution.

Conclusions
Through Definitions 1.1 and 1.2 we have introduced a new definition of singular solution via a double parametrisation which enables us to define a continuous parametrisation even when the limiting profile of the singular solution is not a continuous function. This reduces the problem of constructing singular solutions to a purely algebraic problem. For the DDE (1.1) with two state-dependent delays we constructed three different solution profiles in Section 2 and in Theorems 2.2, 2.3 and 2.4 identified parameter constraints for these unimodal and type I and type II bimodal singular solutions to exist. In Section 3 we investigated the parameter constraints for the singular solutions constructed in the previous section, and treating K 1 as a bifurcation parameter in Theorem 3.1 identified intervals of K 1 for which unimodal singular solutions exist. Theorem 3.3 identifies intervals on which type I bimodal solutions exist, and also a singular fold bifurcation where the solution profile also transitions between unimodal and type I bimodal. Theorem 3.4 identifies intervals on which type II bimodal solutions exist, and a point where the solution profile transitions between unimodal and type II bimodal, with or without a singular fold bifurcation, and we hence identify a singular co-dimension two bifurcation.
The results in Sections 2 and 3 all follow from our definition of singular solution following purely algebraic arguments. Although we do not prove analytically that the singularly perturbed DDE (1.1) has corresponding periodic solutions for 0 < ε 1, in Section 4 we demonstrate numerically using DDEBiftool that the singular periodic solutions that we found do persist for ε > 0. Not only do individual periodic orbits persist, but our theory is robust enough to predict the positions of fold bifurcations of periodic orbits on the bifurcation branches and even the location of a co-dimension two cusp bifurcation associated with these folds. Theorem 3.1(i) and Theorem 3.4 also lead us to numerically find stable periodic orbits with two local maxima per period when ε > 0 (see Figs 4.8,4.12 and 4.16). A period doubling bifurcation also gives rise to stable periodic orbits with up to 5 local maxima per period. In contrast the one delay DDE (1.3) has only been seen to have periodic orbits with one local maxima per period, and no secondary bifurcations on the branches of periodic orbits [19]. We are also able to use our singular solution theory to predict the alignment of the fold bifurcations on different solution branches.
In conclusion, the state-dependent DDE (1.1) has very rich and interesting dynamics in the ε → 0 singular limit, and the concept of singular solution that we introduce in Definitions 1.1 and 1.2 is a useful tool in the study of those dynamics. While we have not proved rigourously that the singular solutions that we construct persist for ε > 0, we have shown numerically that they do, and identified where bifurcations occur. A useful first step in proving convergence as ε → 0 is to identify what the singularly perturbed solutions should converge to. With this work that first step is resolved.

Acknowledgments
Tony Humphries would like to thank John Mallet-Paret and Roger Nussbaum for introducing him to this problem and patiently explaining their results in the one delay case. He also thanks NSERC (Canada) for funding through the Discovery Grants program. Renato Calleja is grateful to NSERC and the Centre de recherches mathématiques, Montréal for funding and to FQRNT, Québec for a PBEEE award. Daniel Bernucci and Michael Snarski thank NSERC for Undergraduate Student Research Awards. Namdar Homayounfar is grateful to the Institut des Sciences Mathématiques, Montréal for an Undergraduate Summer Scholarship.