Regions of Stability for a Linear Differential Equation with Two Rationally Dependent Delays

Stability analysis is performed for a linear differential equation with two delays. Geometric arguments show that when the two delays are rationally dependent, then the region of stability increases. When the ratio has the form 1/n, this study finds the asymptotic shape and size of the stability region. For example, a delay ration of 1/3 asymptotically produces a stability region 44.3% larger than any nearby delay ratios, showing extreme sensitivity in the delays. The study provides a systematic and geometric approach to finding the eigenvalues on the boundary of stability for this delay differential equation. A nonlinear model with two delays illustrates how our methods can be applied.


(Communicated by Peter Bates)
Abstract. Stability analysis is performed for a linear differential equation with two delays. Geometric arguments show that when the two delays are rationally dependent, then the region of stability increases. When the ratio has the form 1/n, this study finds the asymptotic shape and size of the stability region. For example, a delay ratio of 1/3 asymptotically produces a stability region about 44.3% larger than any nearby delay ratios, showing extreme sensitivity in the delays. The study provides a systematic and geometric approach to finding the eigenvalues on the boundary of stability for this delay differential equation. A nonlinear model with two delays illustrates how our methods can be applied.
Violet line gives the real root crossing, while the blue line gives the imaginary root crossing.
The bifurcation analysis of the one-delay version of (1) (B = 0) began with the work of Hayes [24]. The complete stability region in the AC-plane has been characterized by several authors [7,9,10,16,19] with the stability boundary easily parameterized by the delay R (which can be scaled out). Fig. 1 shows this region of stability. The boundary of the stability region for the two-delay equation (1) has been studied by many researchers [1,17,21,22,28,37,33]. Several authors study the special case where A = 0 [20,27,28,37,38,39,40,41,43]. Hale and Huang [21] performed a stability analysis of the two-delay problem, y(t) + a y(t) + b y(t − r 1 ) + c y(t − r 2 ) = 0, (2) where they fixed the parameters, a, b, and c, then constructed the boundary of stability in the r 1 r 2 delay space. Braddock and van Driessche [12] completely determined the stability of (2) when b = c, and partially extended the results outside that special case. Most of these analyses have studied the 2D stability structure of either (1) or (2) with one parameter equal to zero or fixing some of the parameters. Often the 2D analyses result in observing disconnected stability regions for (2). Elsken [17] has proved that the stability region of (2) is connected in the abc-parameter space with fixed r 1 and r 2 . Recently, Bortz [11] developed an asymptotic expansion using Lambert W functions to efficiently compute roots of the characteristic equation for (2) with some restrictions. Mahaffy et al. [32,33] studied (1) for specific values of R, examining the 2D crosssections for fixed A and developing 3D bifurcation surfaces in the ABC-parameter space. The work below shows why delays of the form R = 1 n have enlarged regions of stability. More generally, we illustrate that when R is rational, then the boundary of the stability region is created by a finite set of families of curves. We conjecture that this finite set limits the ability of the boundary of the stability region to converge to a minimum region of stability. Thus, the stability region for rational R is larger.
2. Motivating example. In 1987, Bélair and Mackey [2] developed a two delay model for platelet production. The time delays resulted from a delay of maturation and another delay representing the finite life-span of platelets. The resulting numerical simulation for certain parameters produced fairly complex dynamics. Here we examine a slight modification of their model and demonstrate the extreme sensitivity of the model behavior near rationally dependent delays.  [2]. Simulations show the model sensitivity to the delays R = 1 3 and 1 2 with nearby delays showing complex oscillations.
The modified model that we consider is given by: dP dt = −γP (t) + β(P (t − R)) − f · β(P (t − 1)), where β(P ) = β0θ n P θ n +P n . This model has a standard linear decay term and nonlinear delayed production and destruction terms. The total lifespan is normalized to one, while the maturation time is R. This model differs from the platelet model by choosing an arbitrary fractional multiplier, f , instead of having a delay dependent fraction. For our simulations we fixed γ = 100, β 0 = 168.6, n = 4, θ = 10, and f = 0.35. This gives the equilibria P e = 0 and P e ≈ 5.565 with β (5.565) ≈ 100. Fig. 2 shows six simulations near R = 1 2 and R = 1 3 , where the model is asymptotically stable. However, fairly small perturbations of the delay away from these values result in unstable oscillatory solutions, as is readily seen in the figure. The oscillating solutions are visibly complex. Our modeling parameters were fairly close to the original work after scaling, but our desire to fix f violates the original model's term for population conservation, which in turn allows negative solutions. This paper will explain the high sensitivity in the delay parameter and provide some information on the observed complex oscillations shown in Fig. 2. 3. Background.
3.1. Minimum region of stability and D-partitioning curves. There are a number of key definitions and theorems that are needed to build the background for our study. Our analysis centers around finding the stability of (1). The primary technique follows the D-partitioning method of El'sgol'ts and Norkin [16]. Stability analysis of a linear DDE begins with the characteristic equation, which is found in a manner similar to ordinary differential equations by seeking solutions of the form y(t) = c e λt . The characteristic equation for (1) is given by: This is an exponential polynomial, which has infinitely many solutions, as one would expect because a DDE is infinite dimensional. Stability occurs if all of the eigenvalues satisfy Re(λ) ≤ 0, for all λ.
By inserting λ = 0 into (3), we obtain the plane A + B + C = 0, Λ 0 , which according to the D-partitioning method divides the ABC-parameter space. Crossing Λ 0 results in a real eigenvalue λ passing between positive and negative. The D-partitioning method centers on the geometric image of the characteristic equation with λ = iω in the parameter space. Our analysis concentrates on finding regions of stability for (1) in either the BC-plane with A and R fixed or the ABCparameter space with R fixed. For sake of clarity we briefly review the Method of D-Partitions [16] for our problem with A and R fixed though the concept easily extends to higher dimensions. We have already noted that when λ = 0, the line B + C = −A divides the parameter space so that the region B + C < −A has a positive real root, which implies that (1) is unstable. When λ = iω is inserted into (3), this constraint generates a collection of curves passing through the BCparameter plane, as ω varies. The set of all these curves (hypersurfaces in higher dimensions), including B + C = −A, creates a D-partition of the parameter space based on the exponential polynomial (3), where this collection of curves partitions the BC-parameter space into regions, U k . Every region U k of the D-partition can be assigned a number k, which is the number of zeroes of (3) with positive real part, and the boundary of U k consists of part of the image of (3) with λ = iω, ω ≥ 0. In this decomposition of the parameter space, if there exists a region U 0 , then for all parameters BC in U 0 , (1) is asymptotically stable. It follows by continuity of (3) that crossing any partitioning curve from U 0 (except possibly a degenerate double crossing) results in either one positive root or two complex roots with positive real parts for (3), causing (1) to become unstable. Thus, the regions neighboring U 0 (except at degenerate crossings) are U 1 and U 2 by continuity of (3). This Method of D-Partitions provides a valuable geometric tool for locating regions of stability for delay differential equations [16].
To understand what is meant by rationally dependent delays resulting in larger regions of asymptotic stability, we need an important theorem about the minimum region of stability for (1). The proof of the MRS Theorem can be found in both Zaron [44] and Boese [8]. Note that one face of this MRS is formed by the plane A + B + C = 0, Λ 0 , where the zero root crossing occurs.
The other way that (1) can lose stability is by roots passing through the imaginary axis or λ = iω. This parametrization of λ is substituted into (3) to compute the D-partitioning curves (or hypersurfaces). Since the real and imaginary parts are zero, we obtain a parametric representation of the D-partitioning curves (or hypersurfaces) for B(ω) and C(ω). These are given by the expressions: where (j−1)π 1−R < ω < jπ 1−R , and j ∈ Z + . Clearly, there are singularities for B(ω) and C(ω) at ω = jπ 1−R . This leads to the following definition for D-partitioning surfaces. Definition 3.2. When a value of R in the interval (0, 1) is chosen, a D-partitioning Surface j, Λ j , is determined by Eqns. (4) and (5), and is defined parametrically for (j−1)π 1−R < ω < jπ 1−R and A ∈ R. The surface, Λ j , creates a separate parameterized surface representing solutions of the characteristic equation, (3), with λ = iω, which can be sketched in the ABC coefficient-parameter space of (1), for each positive integer, j. Because the MRS is centered on the A-axis, we often choose to fix A and view the cross-section of the D-partitioning surfaces. Thus, we have the related definition: Definition 3.3. D-partitioning Curve j, Γ j , is determined by Eqns. (4) and (5) and is defined parametrically for (j−1)π 1−R < ω < jπ 1−R with the values of R and A fixed.
The curve, Γ j , creates a parametric curve, which can be drawn in the BC-plane for each j, representing solutions of the characteristic equation, (3), with λ = iω. For A fixed, we let Γ 0 be the D-partitioning curve for λ = 0 in (3). For most values of A, the D-partitioning curves in the BC plane generated by (4) and (5) tend to infinity parallel to the lines B +C = 0 or B −C = 0. When A and R are fixed, one can show using the partitioning method of El'sgol'ts and Norkin [16] that a finite number of D-partitioning curves will intersect in the BC parameter space to form the boundary of the stability region. It is along this boundary where eigenvalues of (3) cross the imaginary axis in the complex plane. Fig. 3 gives three examples showing the first 100 D-partitioning curves for A = 1000 and delays of R = 1 3 , R = 0.45, and R = 1 2 . Assuming that 100 D-partitioning curves give a good representation of the stability region, this figure shows how different the regions of stability are for the different delays. From Theorem 3.1, the interior of the MRS is stable, and the D-partitioning method requires crossing some curve Γ i for (1) to lose stability. Fig. 3 shows an ordering of the D-partitioning curves for R = 1 3 and R = 1 2 , which appear to prevent any Γ i from approaching the MRS on two sides, enlarging the region of stability. Thus, the stability region for R = 1 2 is significantly greater than the others, while the stability region of R = 0.45 is very close to the MRS.
As seen in Fig. 3, the D-partitioning curves can intersect often along the boundary of the region of stability, which creates challenges in describing the evolution of the complete stability surface for (1) in the ABC-parameter space. We need to discuss how we construct the 3D D-partitioning surface using a few more defined quantities.

3.2.
Evolution of the stability surface. Fig. 3 shows that the boundary of the stability region for (1) is very complex, varying significantly as the parameters change. Our goal in this paper is to show how the stability surface changes with the parameters. Generally, R is fixed, and we follow the evolution of the stability region in the BC-parameter space as A increases. Through extensive numerical studies, we have only observed a limited number of ways that the stability region can change. From the observed properties of the D-partitioning curves, we conjecture that our definitions below provide the only means for the boundary of the stability region to change for a range of R values. This study is mostly limited to R ∈ 1 10 , 1 2 . (For example, we expect very different techniques for R → 0 or 1.) To understand the changes in the boundary of stability, we have needed to create some terminology to describe the geometric changes. Below we provide definitions for our terminology and illustrate how the stability boundary evolves in the ABC-parameter space.
In the theorem below, Mahaffy et al. [33] proved that if R > R 0 for R 0 ≈ 0.012, then the stability surface comes to a point with a smallest value, A 0 . For some range of A values with A > A 0 , the stability surface is exclusively composed of Λ 1 and Λ 0 [33]. As ω → 0 + , Λ 1 intersects Λ 0 . The surface Λ 1 bends back and intersects Λ 0 again, enclosing the stability region. As A increases, Λ 2 approaches Λ 1 , and at least for a range of R, Λ 2 self-intersects. In the 2D BCparameter plane, this creates a disconnected stability region, which later joins the main stability surface emanating from the starting point. The A value, where this self-intersecting D-partitioning surface, Λ 2 , joins, is the 1 st transition. Transitions are one of the most important occurrences that affect the shape of the stability surface.
Proposition 1 (BC Transition). At A * j , the D-partitioning curves Γ j and Γ j+1 coincide at the specific point (B * j , C * j ), where Proof. Let A * j satisfy (6) and take the limit as ω → jπ 1−R in (4) and (5). With l'Hôpital's rule the formulas easily follow.
Proposition 2 (Degeneracy Line). Let A = A * j . The characteristic equation (3) is satisfied by the eigenvalues λ = ± jπi 1−R for all values of B and C, such that which we define to be the degeneracy line, ∆ j .
Proof. Insert A = A * j and λ = ± jπi 1−R into (3). It follows that from both the real and imaginary parts. With the definitions of B * and C * from (7), the degeneracy line (8) easily follows.
If Λ j is on the boundary of the stability region for A slightly less than A * j , then ∆ j becomes part of the stability region's boundary, at Transition j. Subsequently, Λ j+1 enters the boundary of the stability region. These transitions create the greatest distortion to the stability surface and attach stability spurs. It is important to note that many transitions occur outside the stability surface and only affect the organization of the D-partitioning curves. Fig. 4 shows cross-sections in the BCplane as A goes through a transition. Definition 3.6 (Stability Spur). The stability spur j, Sp j (R), satisfies the following properties: i The D-partitioning surface, Λ j+1 , self-intersects and encloses a region of stability for (1). ii Sp j (R) connects with the main stability surface at A = A * j via the degeneracy line, ∆ j . iii Self-intersecting Λ j+1 begins at a cusp with A = A p j , and the cross-sectional area monotonically increases with A for A p j ≤ A ≤ A * j , creating Sp j (R). Remark 1. The stability spur, Sp j (R), forms a quasi-cone shaped stable region from the D-partitioning surface, Λ j+1 , which in the BC-cross-sectional plane for A p j < A < A * j creates a disjoint stability region. At a transition A = A * j , this stable region joins the main stability region, allowing Λ j+1 to become part of the main stability region. We note that stability spurs have not been proven to exist analytically. However, extensive numerical studies have established the features described above. The stability spurs, Sp j (R), joining Λ j+1 to the main stability surface, result in the greatest observed distensions to the stability region from Dpartitioning surfaces. We also note that numerical studies show that stability spurs play an increasing role in the stability region as R decreases.
Numerical studies have shown only a couple of other ways for D-partitioning surfaces to enter (or leave) the boundary of the main stability region as A increases. We define these means of altering the boundary as transferrals and tangencies, which relate to higher frequency eigenvalues becoming part of the boundary (or being lost) as A increases.

Definition 3.7 (Transferral and Reverse Transferral). The transferral value of
i,j , Γ j enters the boundary of the stability region. If Γ j is part of the boundary of the stability region for some A <Ã z j,i , a reverse transferral, A z j,i , occurs when Γ j leaves the boundary of the stability region at the intersection of Γ i with Γ 0 . Definition 3.8 (Tangency and Reverse Tangency). The tangency value of A = A t i,j is the value of A corresponding to two D-partitioning curves, Γ i and Γ j , becoming tangent on the boundary of the stability region. For at least some A > A t i,j , Λ j becomes part of the boundary of the stability region, and Γ j separates segments of the D-partitioning curve Γ i to which it was tangent. Often as A is increased, the D-partitioning curve Γ j becomes tangent again to Γ i in a BC cross-section of the stability surface at some A =Ã t j,i , a reverse tangency. Subsequently, as A  Figure 5. Transferral 1, A z 1,6 . The stability boundary for R = 0.25 located in the BC plane before, during and after the first transferral is given by the closed region enveloping the MRS (bold dashed line). That closed region is formed by Γ 0 (violet), Γ 1 (blue), Γ 2 (green), Γ 3 (black), Γ 5 (dark green), and Γ 6 (orange).
increases further, the D-partitioning surface Λ j leaves the boundary of the stability region. Fig. 5 shows an example of the transferral, A z 1,6 , where D-partitioning curve, Γ 6 , enters the stability surface for A > A z 1,6 when R = 1 4 . We can readily see this change in the stability region near where Γ 0 and Γ 1 intersect. Fig. 6 shows an example of a tangency, A t 3,9 , where D-partitioning curve, Γ 9 , enters the stability surface for A > A t 3,9 when R = 1 4 . In this case, Γ 9 becomes tangent to Γ 3 and for larger A values becomes part of the stability boundary. Numerical studies show that the majority of the changes to the stability surface come from tangencies (or reverse tangencies).  Figure 6. Tangency, A t 3,9 . The stability boundary in the BC plane for R = 0.25 is given before, during and after A t 3,9 ≈ 49. The color scheme for the curves is: Λ 0 (violet), Γ 1 (blue), Γ 2 (green), Γ 3 (black), Γ 6 (orange), and Γ 9 (gray). The boundary is the closed region composed of portions of various competing curves that enshroud the MRS (bold dashed line).

4.
Examples from numerical studies. We have developed a number of tools in MatLab to facilitate our stability studies of (1). The ability to rapidly generate D-partitioning curves has led to significant insight into how the stability region evolves in A and R as viewed in the BC-parameter plane. In this section we begin with some 3D stability surfaces for R = 1 5 and R = 1 4 to illustrate how the region of stability varies with A. We present a diagram, which was numerically generated, to illustrate the systematic ordering of transitions, transferrals, and tangencies for a range of R values. Finally, we end this section with detailed numerical studies for specific values of R as R → 1 4 . Fig. 7 provides 3D views of the stability region of (1) for R = 1 5 and R = 1 4 with A ≤ 21. For R = 1 5 , there are only three transitions with a 4 th transition occurring at A = +∞. Fig. 7 shows the starting point of the stability surface at −6, − 1 4 , 25 4 . Initially, Λ 1 (blue) and Λ 0 (violet) bound the stability region. Next the first stability spur (green) enters and attaches Λ 2 to the stability region at the first transition, A * 1 . Subsequently, Λ 3 (black) and Λ 4 (red) adjoin the boundary of stability. At A ≈ 21 a transferral occurs bringing Λ 7 onto the boundary, and around A ≈ 70, there is a tangency of Λ 11 interrupting Λ 4 . For this stability surface with no other transitions affecting the boundary, we only observe additional tangencies, which change the boundary of the stability surface for larger A. Fig. 7 shows the MRS (gray) interior to the stability surface, and visually it is apparent how much the transitions and stability spurs distort the stability surface away from the MRS. It is worth noting that the stability spurs are shrinking in size as A increases. Fig. 7 with R = 1 4 is initially very similar to the case with R = 1 5 . However, when R = 1 4 , there are only two transitions, so no additional stability spurs occur. The next two changes to this surface are a transferral, A z 1,6 with Λ 6 (orange) entering around A = 14, then a tangency, A t 3,9 with Λ 9 entering around A = 50. Numerically, we can easily track the changes as A increases. Readers can access the MatLab figures for additional viewing angles and programs for generating the D-partitioning curves 1 .
with the other two transitions, their stability spurs, and one transferral all occurring before A = 15. Afterwards, when R = 1 4 , all remaining changes to the stability surface occur through tangencies, reverse tangencies, or a reverse transferral. None of these events dramatically change the shape of the boundary of the stability region. What is significant is that the 3 rd transition does not occur until A * 3 = +∞, causing the distortion of this transition to persist, while nearby delays do not have this distortion for sufficiently large A.

Stability region near
The characteristic equation (3) is an analytic function, so there is continuity of the stability surfaces as the parameters vary. To study what happens to the stability surface for R = 1 4 , we explore in some detail the stability surface for R = 0.249. Not surprisingly, there are many similarities between these surfaces until the singularity occurs at the transition, A * 3 = 749.93 for R = 0.249. We provide details of the evolution of stability surface for R = 0.249, which suggests why the region of stability for R = 1 4 remains larger than the MRS as A increases. at R = 0.249, we see there is a transferral, A z 1,6 ≈ 13.3. At this stage, the BC cross-section of the stability surface is very similar to the images in Fig. 5. The boundary at the transferral is comprised of Γ 0 , Γ 1 , Γ 2 , and Γ 3 . Subsequently, Γ 6 enters the boundary near the intersection of Γ 0 and Γ 1 .
The next change in the stability surface for R = 0.249 is a tangency, which occurs at A t 3,9 ≈ 49.4. This tangency has little effect on the shape of the boundary of the stability region, but allows higher frequency eigenvalues to participate in destabilizing (1). This tangency is very similar to the one depicted in Fig. 6, occurring in the 1 st quadrant. As A increases, a series of tangencies happens, alternating between the 1 st and 4 th quadrants. There are a total of 11 tangencies that occur before A * 3 ≈ 749.93, with the last one being A t 33,39 ≈ 462.063. Following A t 33,39 , all of these tangencies undergo a reverse tangency (in the reverse sequential order) with higher frequency eigenvalues leaving the boundary of the region of stability. Table 1 summarizes all of these events. The onset of reverse tangencies, which occur prior to A * 3 , create significant expansion of the stability region, primarily in the 1 st and 4 th quadrants of the BC-plane. At the same time the stability surface becomes much larger than the MRS. We note that rapidly after A * 3 , there are many tangencies, which occur for increasing A and reduce the size of the boundary of stability for R = 0.249. This results in the region of stability shrinking back to being near the MRS for large A. Numerically at A = 10, 000, our program with a large number of D-partitioning curves visually shows little difference between the region of stability for R = 0.249 and the MRS. Since R = 0.249 is rational, we conjecture that the boundary of the stability region never asymptotically approaches the MRS, yet it is substantially closer than for R = 1 4 . For R = 0.249 at A * 3 = 749.93, the boundary of the region of stability is reduced to only five D-partitioning curves. There are lines from Λ 0 (violet) and the degeneracy line, ∆ 3 , with λ = ± 3iπ 0.751 at A * 3 . The boundaries in the 1 st and 4 th quadrants are formed from Γ 3 (black) and Γ 1 (blue), respectively. Finally, there is a very small segment of the boundary formed by Γ 2 (green). This stability region is shown in Fig. 9. We see distinct bulges away from the MRS in the 1 st and 4 th quadrants, increasing the size of the region of stability. This simple boundary easily allows the computation of the area of the region of stability. A numerical integration gives that the region of stability outside the MRS is approximately 26.85% of the area of the MRS, which is a substantial increase in the region of stability.

5.
Analysis. The objective of this section is to convince the reader that the asymptotic stability region for R = 1 n with A → +∞ reduces to a simple set of curves bounded away from the MRS. As R → 1 n , there is a critical transition with A * n−1 → +∞, which leaves the boundary of the stability region for (1) composed of only Λ 0 , Γ 1 , Γ n−1 , and the degeneracy line, ∆ n−1 . (There may be small segments of other D-partitioning curves for any R < 1 n at A * n−1 .) Mahaffy et al. [32,33] showed that as ω → 0 + , Γ 1 intersects Λ 0 along the line For now assume that as R → 1 n from below, the shape of the stability region at A * n−1 is very similar to Fig. 9 (or Fig. 10) with a very limited number of Dpartitioning curves on the boundary of the stability region. The length of any edge of the MRS is easily seen to be A √ 2. We examine the edge of the stability region composed of the D-partitioning curve, Γ 0 , where real roots cross. The upper left  The cases for R = 1 2n and 1 2n−1 give different shaped regions, but ultimately as A → +∞, the regions are bounded by only four curves with two lying on the boundary of the MRS and two bulging away from this region. We prove analytically the position of some key points, then rely on the limited number of families of curves and their distinct ordering to give the simple structure. The numerically observed orderly appearance and disappearance of the tangencies on the boundary of the region of stability provide our argument for the ultimate simple structure of the stability region and the continued bulge away from the MRS for these particular rational delays. From the monotonicity of the transition curves, we note that all limiting arguments require R approaching 1 n from below.
5.1. Stability region for R = 1 2n . In this section we provide more details on the increased size of the region of stability for delays of the form R = 1 2n as A → +∞. Earlier we presented evidence that as R → 1 4 , the third transition, A * 3 , tended to infinity, and at that transition the boundary of the stability region in the BC-plane reduced to just four curves. Λ 0 provides the lower left boundary along the MRS in the 3 rd quadrant. With R < 1 4 , ∆ 3 , which occurs at A * 3 with eigenvalues λ = 3π 1−R , creates a boundary parallel to the MRS in the 2 nd quadrant. This line approaches the MRS as R → 1 4 from below. Γ 1 forms a boundary in the 4 th quadrant, which significantly bulges away from the MRS. Its intersection at Λ 0 extends 1 3 times the length of the MRS along the line of this plane into the 4 th quadrant. Γ 3 creates an almost mirror image across the C-axis in the 1 st quadrant, bulging away from the MRS the same distance. Fig. 9 illustrates this expanded stability region very well. Below we prove some results about the four curves on the boundary of the stability region and give additional information on why we believe the stability region has its enlarged character.
The case R = 1 4 extends generically to the case R = 1 2n . Λ 0 is always one part of the stability boundary. Symmetric to this boundary across the B-axis is ∆ 2n−1 at A * 2n−1 , which occurs with A * 2n−1 → ∞ as R → 1 2n . Below we demonstrate that ∆ 2n−1 approaches the line C − B = A * 2n−1 in the BC-plane as A * 2n−1 → ∞, which is one side of the MRS. The other two sides are composed of Γ 1 in the 4 th quadrant and symmetric to this, Γ 2n−1 in the 1 st quadrant.
To help obtain the symmetrical shape discussed above (and exclude the Dpartitioning curves that pass near the point (B, C) = (−A * 2n−1 , 0)), there are alternate forms of the Eqns. (7). We multiply and divide the cosecant terms by cos jRπ 1−R , then use the definition of A * j to obtain Proof. For R < 1 2n and R → 1 2n , we consider the transition A * 2n−1 . The degeneracy line, ∆ 2n−1 , satisfies Since j = 2n − 1 and A * 2n−1 = − (2n−1)π 1−R cot (2n−1)Rπ 1−R , Eqns. (10) give , so ∆ 2n−1 tends towards one edge of the MRS. We note that the small distance between ∆ 2n−1 and the MRS may allow a very small segment of Γ 2 to be part of the stability region for all R < 1 2n .
It follows that as R → 1 2n from below, then A * 2n−1 → ∞ and Γ 2n−1 passes arbitrarily close to the point (A * 2n−1 , 0). The lemmas above prove some of the features of the stability region in Fig. 9 for R = 1 2n with A → +∞. In particular, we see the left boundaries aligning with the MRS in the 2 nd and 3 rd quadrants. We also proved that Γ 1 and Γ 2n−1 intersect near B = A * 2n−1 when C = 0. Finally, we proved that Γ 1 and Γ 2n−1 intersect Λ 0 and ∆ 2n−1 , respectively, in a symmetric manner at as R → 1 2n from below and A * 2n−1 → +∞. It remains to show that Γ 1 and Γ 2n−1 are the only D-partitioning curves on the boundary as R → 1 2n from below and A → +∞. In this work we will only present numerical evidence for this result, providing some ideas for how a rigorous proof might proceed. The numerical results will include how much larger the asymptotic region of stability is for rational delays of the form R = 1 2n .
As noted earlier, we developed a MatLab code for easily generating and observing D-partitioning curves at various values of A and R in the BC-plane. Mahaffy et al. [33] showed that the rational delays, R, due to the periodic nature of the sinusoidal functions, result in the D-partitioning curves ordering themselves into families of curves.
Remark 2. The equations of Definition 5.4 show that B i+2j (µ) follows the same trajectory as B i (ω) with a shift of 2nπ cos( kω n )/ sin( jω n ) for ω ∈ (j−1)π 1−R , jπ 1−R , while C i+2j (µ) follows the same trajectory as C i (ω) with a shift of 2nπ cos(ω)/ sin( jω n ) over the same values of ω. Thus, there is a quasi-periodicity among the D-partitioning curves when R is rational.
The organization of these families of curves and systematic transitions allow one to observe how the different D-partitioning curves enter and leave the boundary of the stability region. (See Fig. 8.) We provide more details following some analytic results for delays of the form R = 1 2n+1 .

5.2.
Stability region for R = 1 2n+1 . When R = 1 2n+1 and A → +∞, the stability region again bulges away from the MRS. The stability region again asymptotically reduces to just four D-partitioning curves. However, the regions of stability, which lie outside the MRS, now appear in the 2 nd and 4 th quadrants. Fig. 10 shows R = 0.199 and A * 4 = 799.9 with primarily four D-partitioning surfaces comprising the boundary of the region of stability. This figure is quite representative of any R → 1 2n+1 from below as A * 2n → +∞. Below we present several lemmas, which analytically show results similar to Lemmas 4.1-3 for R = 1 2n . The next section will complete the study with numerical results.
For any delay, R, Λ 0 is always one part of the stability boundary, appearing in the 3 rd quadrant of the BC-plane along the MRS with B + C = −A. When R → 1 2n+1 from below, there is a transition A * 2n → ∞, which results in a degeneracy line that approaches is parallel to Λ 0 , and lies on the opposite side of the MRS. As in the previous case (R = 1 2n ), Γ 1 in the 4 th quadrant creates another edge of the stability region. Finally, the stability region for R = 1 2n+1 , asymptotically finds Γ 2n mirroring Γ 1 in the 2 nd quadrant to complete the simple enlarged stability region.

STABILITY FOR A TWO DELAY DE 4973
The next lemma shows that as R → 1 2n+1 from below and A * 2n → +∞, the point (B * 2n , C * 2n ) tends to a value symmetric with the origin to the intersection of Γ 1 and Λ 0 .
Lemma 5.6. For R < 1 2n+1 and R → 1 2n+1 , the D-partitioning curve Γ 2n comes to the point Proof. This proof is very similar to the proof of Lemma 5.2, so we omit it here. Figure 10 shows that Γ 1 and Γ 2n cross the C-axis near the corners of the MRS. The next lemma proves this asymptotic limit, giving more information about the symmetric shape of the region.
Lemma 5.7. For R < 1 2n+1 and R → 1 2n+1 , the D-partitioning curves, Γ 1 and Γ 2n , pass arbitrarily close to the points (A * 2n , 0) and (−A * 2n , 0), respectively, in the BC-plane with A * 2n → +∞. Proof. The argument for Γ 1 passing arbitrarily close to (A * 2n , 0) is almost identical to the argument given in Lemma 5.3. Similarly, Γ 2n has (2n−1)π 1−R < ω < (2n)π 1−R , which has ω = 2nπ inside the interval. Since this is an even multiple of π, the cos(ω) → 1. Otherwise, the arguments parallel Lemma 5.3 and so Γ 2n passes arbitrarily close to (−A * 2n , 0). In summary of the analytic results, we have established convergence and symmetry of key points for the enlarged region of stability for R = 1 n , as depicted in Figs. 9 and 10. We showed the degeneracy line, ∆ n−1 approaches the MRS either parallel for n odd or perpendicular for n even to Γ 0 . The region of stability extends out to where Γ 1 intersects Γ 0 and a point symmetric to the origin for n odd and symmetric to the B-axis for n even. Some limited results have proved special cases that other related family curves lie outside these symmetric points [42]. Finally, we proved that Γ 1 intersects ∆ n−1 (n odd) or Γ n−1 (n even) near the corner of the MRS on the C-axis. Similarly, Γ n−1 intersects Γ 0 near one corner of the MRS when n is odd. There remain significant analytic results to prove the proposed shape of the enlarged regions of stability. In particular, our results, except for ∆ n−1 , occur only at points instead of complete curves. Also, the simplified shape with primarily four curves has only been shown at A * n−1 for R < 1 n and not extended to the limit, A → +∞, when R = 1 n . The next section provides more numerical details to support our claims of increased stability regions for delays of the form R = 1 n with A large.
6. Asymptotic stability region for R = 1 n . The previous section presents the simple stability regions for R near 1 n , and we gave some details for R = 0.249 on how the stability region evolves as A increases. In this section we present more details from numerical studies to convince the reader of the enlarged stability region for these specific rational delays and provide some measure of the size increase compared to the MRS.
Geometrically, for a fixed A the MRS is a square region with one side bounded by Γ 0 . Definition 5.4 shows that rational delays, R, create families of smooth curves with similar properties and following similar trajectories. The limited number of family members for a fixed R = 1 n creates a type of harmonic, which prevents the complete collection of D-partitioning curves from asymptotically approaching all sides of the MRS. Hence, the stability region for R = 1 n is enlarged with larger regions for smaller values of n. This section provides some details and numerical results for the structure and size of the asymptotic stability region.
6.1. R near 1 4 . To illustrate our analysis we concentrate on the case R = 1 4 and note that the arguments generalize to other delays of the form R = 1 n . For R = 1 n , Definition 5.4 gives 2(n − 1) distinct family members, which implies R = 1 4 has six distinct family members. As noted above the key D-partitioning curves on the boundary of the stability region asymptotically as A * 3 → +∞ with R → 1 4 from below are Γ 1 and Γ 3 with Λ 0 and ∆ 3 , creating the other two boundaries along the MRS. Γ 1 and Γ 3 are obviously the first D-partitioning curves of the first and third families, while ∆ 3 arises from the singular point ω = 3π 1−R between Γ 3 and Γ 4 . Fig. 11 shows 60 D-partitioning curves for R = 0.249 at A * 3 = 749.93. By continuity of the characteristic equation (3), we expect similarities between R = 1 4 and R = 0.249, particularly for A < A * 3 . Fig. 11 shows clearly the six family structure we predict for R = 1 4 . (Note that R = 0.249 is predicted to have 1502 families by Definition 5.4, which will ultimately result in a much closer approach of the stability region to the MRS as A → +∞.) The figure shows the distinct ordering of the family members within the six families and the characteristic pattern of each of the six families. The coloring pattern in Figure 11 follows Λ 0 in violet and then the six successive families in blue, green, black, red, gray, and orange. Fig. 9 shows the five curves on the boundary of the stability region. As noted before, Λ 0 is always on the boundary in the 3 rd quadrant. It connects to Γ 1 , the 1 st member of the first family. The first close-up in Fig. 11 shows the organization of the first family with all members lying outside the boundary of the region of stability with each successive member further away. This first close-up also shows the 6 th family, which parallels the first family along the boundary in the 4 th quadrant, then diverges opposite the first family below Λ 0 . Thus, in the 4 th quadrant the boundary of the stability region consists only of the curves from Λ 0 and Γ 1 .
In the 2 nd quadrant, the primary boundary is ∆ 3 . By Lemma 5.1, as R → 1 4 , ∆ 3 approaches the line C − B = A * 3 , which is visible in Fig. 11. Since R = 0.249 < 1 4 , there is a small gap between ∆ 3 and the boundary of the MRS, so we see a small segment of Γ 2 on the boundary of the stability region, which is visible in the final close-up of Fig. 11. The line ∆ 3 extends into the 1 st quadrant. We see that outside the tiny segment of Γ 2 (left of the MRS), the second and fifth families are outside the stability region running symmetrically about the C-axis through the 2 nd and 3 rd quadrants.
In the 1 st quadrant, we see Γ 3 on the boundary of the stability region. The third close-up of Fig. 11 shows the ordering of the third family with all members successively outside Γ 3 . This figure also shows all members of the fourth family outside (and intertwined with) the third family, paralleling each other in the first quadrant. The fourth family diverges opposite the third family parallel to ∆ 3 .
Numerically, our MatLab programs allow us to observe the stability region for any delay R slightly less than 1 4 at A * 3 , and the boundary of the stability region is almost identical to Fig. 11, except that A * 3 increases with R → 1 4 , expanding the scales of the B and C axes. To further understand the evolution of this stability surface as R → 1 4 from below, we detail how Fig. 8 can be used to explain the process. As noted earlier, the key elements causing the bulges in the stability region are the limited number of families of D-partitioning curves and the transitions, which distort the boundary. At R = 1 4 , all the transitions A * 3n → +∞, n = 1, 2, ..., which is easily verified by Definition 6. Furthermore, it can be shown using perturbation analysis that the transitions have a distinct ordering, .. for delays R < 1 4 . (Details of this proof are not included, but depend on n and R, as would be expected [42].) This sequence allows one to readily determine changes to the boundary of the stability region. Table 1 provides the complete evolution of the stability surface for R = 0.249 from its beginning at A 0 ≈ −5.016 until A * 3 ≈ 749.93. We call attention to the sequence of reverse tangencies and the reverse transferral for A < A * 3 . These events region is relative to the MRS. Table 2 gives the relative increase of this region of stability for several values of R → 1 4 , indicating the predicted asymptotic increase in size of the stability region for R = 1 4 at A * 3 = +∞. Since the transitions, A * 3n , all occur at +∞ when R = 1 4 , continuity of the characteristic equation suggests this enlarged stability region persists for R = 1 4 and should be 26.86% larger than the MRS.
Area Ratio 0.249 749.9 1.2687437 0.24999 75,000 1.2686377 0.2499 7500 1.2686388 0.249999 750,000 1.2686377 Table 2. Numerical computation of the region of stability at A * 3 near R = 1 4 . Area given as ratio of stability region to MRS.
For any R < 1 4 and A > A * 3 , the six family structure breaks down, leading to new tangencies and a new ordering of the larger number of families, which results in significantly smaller regions of stability. We noted that R = 0.249 has 1502 families of curves, so after A * 3 these many families each approach the MRS in slightly different ways. With this large number of families of curves, the D-partitioning curves can very closely approach the MRS, which is impossible for R = 1 4 with only six families of curves. (The reader should examine Fig. 3, where R = 0.45, which has 22 families of curves and very closely approaches the MRS.) We conjecture that for R = 0.249 as A → +∞, the stability region closely approaches the MRS. However, since 1502 families of curves is finite, there is not a perfect asymptotic approach to the MRS though visually the stability region for large A and R = 0.249 will appear the same as the MRS (like Fig. 3 for R = 0.45). Similarly, this will be the case for any rational R near R = 1 4 with A A * 3 . For irrational R there are an infinite number of families of curves, which would allow asymptotic approach to the MRS as A → +∞.
6.2. Increased area for R = 1 n . Figs. 9 and 10 show the simple D-partitioning curve structure on the typical boundary of stability region for R near 1 n at A * n−1 . The symmetrical shape varies depending on whether n is even or odd, but all these stability regions at A * n−1 reduce to having Λ 0 on one edge, ∆ n−1 on another, and the D-partitioning curves Γ 1 and Γ n−1 comprising the majority of the remaining two edges of the boundary of the region of stability. (The gap between ∆ n−1 and the MRS allows small segments of Γ 2 and possibly Γ n−2 to remain on the boundary of the stability region at A * n−1 , shrinking as R → 1 n .) Γ 1 and Γ n−1 bulge out from the MRS, leaving an increased region of stability, which is readily computed.
For R = 1 n , Def. 5.4 gives 2(n − 1) families of D-partitioning curves. These families organize much in the same way as shown in the previous section (R near 1 4 ) to help maintain the increased regions of stability for R = 1 n over the MRS. The orderly family structure allows one to study each R = 1 n much as we did in the previous section and observe similar sequences of transferrals, tangencies, reverse tangencies, and reverse transferrals, which ultimately lead to the simple structure of the stability region seen in Figs. 9 and 10 at A * n−1 for R < 1 n , sufficiently close. Using the continuity (pointwise) of the characteristic equation, we claim that the bulge in the region of stability persists for R = 1 n .
R Area Ratio Linear Extension R Area Ratio Linear Extension 1.1386 0.2000 Table 3. Increases in the Region of Stability for R = 1 n . Area ratio is the ratio of the stability area to the MRS. The linear extension is the ratio of how far the stability region extends past the MRS.
We showed the region of stability extends linearly along the MRS by a factor of 1 n−1 for R = 1 n . The simple D-partitioning curve structure allows easy numerical computation of this area bulging from the MRS. Table 3 gives the size of the increased region of stability for various R = 1 n . Asymptotically, the region of stability for R = 1 2 is triangular with only Λ 0 , Γ 1 , and ∆ 1 . This results in a region that is twice the size of the MRS, asymptotically. As the denominator increases, the asymptotic region of stability decreases relative to the MRS, yet it is still over 10% larger than the MRS when R = 1 7 . In a small neighborhood about R = 1 n , all values of R = 1 n are either irrational or have a reduced rational form with a large number of families of curves. We conjecture that given any R = 1 n in this neighborhood, then we can find an A sufficiently large that the corresponding region of stability is close to the MRS. For example, if we choose a delay, R , close to R = 1 3 , then for A sufficiently large the region of stability for R = 1 3 is approximately 44.3% larger than the region of stability for R , according to Table 3, since the region of stability for R is close to the MRS. 6.3. Stability spurs. The discussion above shows how transitions increase the area of the region of stability for R = 1 n . Just before a transition, A * j , D-partitioning curve Γ j extends out toward D-partitioning curve Γ j+1 , expanding the main region of stability. Numerically, we observe that just prior to A * j , Γ j+1 self-intersects creating an island of stability that is disconnected in the BC-plane for a fixed A. Definition 3.6 describes the stability spurs, which connect to the main stability surface at transitions, A * j . In this section we provide some details from our numerical simulations about the stability spurs. Numerically, we observe exactly n − 1 stability spurs for R ∈ 1 n+1 , 1 n . The largest stability spurs, Sp 1 (R), correspond to A * 1 for at least a limited range of R. We performed an extensive study of the stability spurs for R ∈ 1 5 , 1 4 . Over this range there appear to be exactly three stability spurs. Fig. 13 shows the variation in length of the three stability spurs, which are all monotonically decreasing in R. The length of the first spur for R ∈ 1 5 , 1 4 ranges from |Sp 1 (0.2)| = 0.4064 adjoining the main stability surface at A * 1 (0.2) = −3.927 to |Sp 1 (0.25)| = 0.2840 adjoining the main stability surface at A * 1 (0.25) = −2.418. Where this first stability spur adjoins the main stability surface, the cross-sectional area in the BC-plane of Sp 1 (0.2) is 17.65% of the total stable region, while the area of Sp 1 (0.25) at A * 1 (0.25) is 7.08% of the total stable cross-section. (See Fig. 4.) The sizes of the second and third stability spurs are significantly smaller. The length of Sp 2 (0.2) is 0.1136, while the length of Sp 2 (0.25) is 0.0305. The crosssectional area of Sp 2 (0.2) is 0.3003% of the total stable region at A * 2 (0.2) = 0, while the area of Sp 2 (0.25) is 0.01483% of the total stable region at A * 2 (0.25) = 4.837. the length of Sp 3 (0.2) is 0.0110, and this spur adjoins the main stability region at A * 3 (0.2) = 11.78. Clearly, no spur occurs at R = 1 4 , as A * 3 (0.25) = +∞.
The stability spurs are easy to observe in Fig. 7 for R = 1 5 and R = 1 4 . From a numerical perspective the length of a stability spur is difficult to accurately compute because of the cusp point, A p j , which is a singularity. Both the length and the selfintersection of Sp 3 (R) as R → 1 4 becomes impossible to compute, and effectively, Sp 3 (R) vanishes as R → 1 4 . However, as we have already seen, the main stability region bulges out in the 1 st quadrant because of A * 3 (R), but we have been unable to confirm that Γ 4 always self-intersects as A * 3 (R) → +∞ with R → 1 4 .
7. Example continued. In Section 2, we examined an example motivated by work of Bélair and Mackey [2], and Fig. 2 showed how rational delays stabilized the model. Here we use some of the information above to provide more details about the complex behavior observed in Fig. 2. With the parameters given in Section 2, the model is linearized about the nontrivial equilibrium, P e ≈ 5.565. If y(t) = P (t)−P e , then the approximate linearized model becomes: This is Eqn. (1) with (A, B, C) = (100, 35, −100). We use our MatLab program to generate plots of 40 D-partitioning curves in the BC-plane with A = 100 for various values of R. Fig. 14 shows the distinct four family feature for the delay R = 1 3 and the enlarged region of stability. The linearized model is in the region of stability, which agrees with the numerical simulation in Fig. 2, where the delay R = 1 3 gives a stable solution. When the delay is decreased to R = 0.318, there is a transition at A * 2 = 42.8. The four family structure rapidly unravels, and the D-partitioning curves begin approaching the MRS more closely. Fig. 15 shows that Eqn. (12) has its equilibrium outside the curves Γ 9 , Γ 13 , and Γ 17 . With the help of Maple (using information from Fig. 15), the eigenvalues of Eqn. The leading eigenvalue comes from the equilibrium point being furthest from Γ 13 , and its frequency suggests a period of 2π/58.36 ≈ 0.108, which is close to the period of oscillation seen in Fig. 2 for R = 0.318. A similar analysis can be performed for R = 0.34. Since this delay is larger than R = 1 3 , there is not the four family structure. Fig. 16 shows four D-partitioning curves, Γ 8 , Γ 12 , Γ 16 , and Γ 19 , between the region of stability and where the linearized model is located. As before, it is easy to use this information to determine the eigenvalues: The equilibrium point is furthest from Γ 12 , and the period from λ 1 is 2π/53.51 ≈ 0.117, which is similar to the period of oscillation seen in Fig. 2 for R = 0.34. The next two eigenvalues are moderately large, resulting in the additional irregular structure observed in the simulation. Note that the last eigenvalue with positive real part has just barely crossed the imaginary axis, which again is apparent from Fig. 16. The modified platelet model also considered delays near R = 1 2 . The simulation in Fig. 2 showed the stability of the solution at R = 1 2 . For R = 1 2 , there are only two families of D-partitioning curves, leaving a fairly large region of stability. This is apparent in the leftmost graph of Fig. 17, where the point for the model parameters is clearly inside the region of stability. When the delay is decreased to R = 0.48, there is a transition, A * 1 = 24.5. Thus, when A = 100, the two family structure has broken down, and the D-partitioning curves form a very different pattern. The result is shown in the middle graph of Fig. 17, where the D-partitioning curves Γ 9 , Γ 11 , and Γ 13 are visible between the model parameter point and the region of stability. Once again, it is easy to compute the eigenvalues with positive real part for this case giving: The leading eigenvalue, λ 1 , has a frequency of 64.71, which suggests a period of 0.09710. This is consistent with the observed period of oscillation in Fig. 2. We observe that this oscillatory behavior is irregular, which reflects a strong contribution from λ 2 and λ 3 , which have higher and lower frequencies, respectively. For R = 0.51, the rightmost graph of Fig. 17 shows the D-partitioning curves Γ 10 and Γ 12 between the model point and the region of stability. When the eigenvalues are computed, we obtain: λ 1 = 0.03415 ± 60.02 i λ 2 = 0.02930 ± 72.28 i.
The frequency of the leading eigenvalue is 60.02, which yields a period of 0.1047. Again, this is consistent with the observed oscillations in Fig. 2. 8. Discussion and conclusion. Delay differential equations (DDEs) with multiple time delays are used in a wide array of applications. When studying the stability of two delay models, our results show the high sensitivity of the DDE for certain delays (rationally dependent) for some ranges of the parameters. In particular, the DDE (1) shows unusually large regions of stability for R = 1 n , when n is a small integer. For example, when R = 1 2 , the region of stability doubles over the Minimum Region of Stability (MRS), which is independent of the delay. Our geometric approach allows a systematic method for visualizing the region of stability and provides a simplified understanding for how the stability region evolves. In our motivating example, we demonstrated how easily the leading eigenvalues could be found, which helped explain some of the observed behavior in the nonlinear problem.
The characteristic equation for (1) is an exponential polynomial, which is deceptively complex to analyze. For rational delays, R, this characteristic equation organizes into families of D-partitioning curves, which enter or leave the boundary of stability in only a limited number of ways. The most significant changes occur at values of the parameter, A, which we defined as transitions, A * j . At a transition, a new D-partitioning curve (higher frequency eigenvalue) enters the boundary of the stability surface and does so by distorting the boundary outward, enlarging the stability surface. As R → 1 n − , the transition A * n−1 → +∞, moving the accompanying distortion further away and maintaining an increased region of stability for larger values of A. We numerically showed that for any R < 1 n , but close, the boundary of the stability region in the BC-cross section at A * n−1 reduces primarily to four simple curves, which allowed us to numerically compute the size of this increased stability region. One interesting phenomenon that can occur at a transition is a "stability spur," which can lead to interesting disconnected regions of stability in the BC-cross sectional parameter space.
The evolution of the stability region had limited, yet very orderly ways of changing for R near 1 n . This is the quasi-periodicity of the families of D-partitioning curves, which intersect mostly through tangencies. The organization of the curves, as seen in many of the figures, produced clear patterns that could be carefully analyzed for R = 1 n , resulting in the observed larger regions of stability. Because any rational delay has a finite set of families of curves, as given in Def. 5.4, we conjecture that this finite set is incapable of asymptotically approaching the MRS. Thus, any rational delay results in a larger region of stability for large A. Irrational delays would not have this organization of the D-partitioning curves, which would allow some curve to approach any point along the MRS for large A. Because R = 1 2 has only two families of curves, this delay should result in the largest region of stability for large A.
In this paper we analytically proved some results to support our claims. Our analytic results are limited to convergence of "degeneracy lines" to the MRS and the ordering of sequences of particular points of interest. These analytic results prove some key features of the enlarged stability region, including the length of distortion for R = 1 n and some symmetry properties. However, significantly more work is needed to prove results about the complete D-partitioning curves on the boundary and some asymptotic properties. Extensive numerical studies and the intriguing ordering of the families of curves suggest that the entire D-partitioning curves follow the ordering of the pointwise sequences. Still more analytic work is needed around the singularities that occur in the characteristic equation for R = 1 n .

JOSEPH M. MAHAFFY AND TIMOTHY C. BUSKEN
Furthermore, other rational delays, like R = 2 5 , show similar increases in their regions of stability, but we have not investigated the details on how these rationally dependent delays produce larger regions of stability. We have produced a framework for future studies of the DDE (1) and have several MatLab programs for further geometric investigations. (Readers can access programs to generate the D-partitioning curves and the 3D figures for additional viewing angles at http://www-rohan.sdsu.edu/~jmahaffy/two_delay_supplement.html.) Understanding the stability properties of DDE (1) is very important for a number of applications with time delays. Our results show that selecting delays of R = 1 n for n small in a model could give the investigator stability that is easily lost with only a very small change in the delay. The ultra-sensitivity in the model can be explained by our results. This two delay problem is very complex, but our geometric techniques provide a valuable tool for future stability analysis of delay models.