CHARACTERISTIC ROOTS FOR TWO-LAG LINEAR DELAY DIFFERENTIAL EQUATIONS

. We consider the class of two-lag linear delay diﬀerential equations and develop a series expansion to solve for the roots of the nonlinear char- acteristic equation. The expansion draws on results from complex analysis, combinatorics, special functions, and classical analysis for diﬀerential equa- tions. Supporting numerical results are presented along with application of our method to study the stability of a two-lag model from ecology.


1.
Introduction. Delay Differential Equations (DDEs) naturally arise in many fields of science and engineering 1 and accordingly have been the focus of extensive study over the years. The classic books by Bellman and Cooke [10] and El'sgol'ts and Norkin [26] describe the significant mathematical results up to the 1960's, while subsequent books by Hale and Verduyn Lunel [33] and Diekmann et al., [25] cover the later semigroup-based frameworks. A specific area in which delays arise frequently is biological modeling, e.g., due to gestation periods [13] or biochemical delays [38], and there are several good texts on population modeling and analysis with DDEs in biology [22,30,40,45] including a recent introductory textbook by Smith [59]. While the focus of the illustrating example in this work is from biology, DDEs also arise in other areas including spindle speed milling [44] and compound optical resonators [54]. Lastly, we direct the interested reader to several recent review articles [31,51] and edited works [6,43,58] which thoroughly cover the stateof-the-art advances in the study and application of DDEs.
1.1. DDE stability analyses. Consider the first order, linear, single-lag DDĖ where φ(t) is the initial history function (which need not satisfy (1)). The characteristic equation of (1) is a transcendental equation with an infinite number of roots {s j } in the complex plane (see Figure 1 for an example). Thus the solutions are constructed as an infinite series and we note that these basis functions e sj t form a Schauder basis for L p functions over a finite domain. The coefficients C j can be computed via evaluating the history function φ(t) and the basis functions e sj t at 2N + 1 timepoints in [−τ, 0] and then linearly solving for 2N + 1 of the C j 's. The transcendental equation in (3) belongs to a well-known class of equations known as exponential polynomials [9,49,52] and many researchers in DDEs have (understandably) chosen to study them [5,19]. Moreover, they arise quite naturally in the Evans functions used to determine stability of traveling waves for PDEs with time delays [55]. In the broader mathematics community, research into this class of algebraic equations is very active and extends to numerical schemes for computing the roots (see papers by Jarlebring such as [34]) and applications to fields afar as quantum computing [56]. Regardless of their origin, the computation and estimation of properties of the roots to equations such as (3) is clearly important in many ares of applied mathematics.
A primary goal in the study of linear (and linearized) DDEs is to identify the value of the the delay τ for which equilibria become unstable, i.e., when the real part of the principal root Re(s 0 ) becomes positive. A traditional strategy for proving that such a τ exists is to separate (3) into real and imaginary equations and prove that for the τ * such that s 0 (τ * ) is purely imaginary, the root passes through to the positive real half of the plane, i.e.,

1.2.
Relationship between DDEs and Lambert W function. In their 2003 article, Asl and Ulsoy [3] connected the roots of (3) with solutions to a specific exponential polynomial known as the Lambert W equation [20,21]. Briefly: to solve for the roots to (3), rewrite the equation in the form and define a new function W such that W (βτ e −ατ ) = τ (s − α). The equation to be solved is now which is in the same form as W (x)e W (x) = x, the well-studied Lambert W equation. The roots to (3) are thus The principal solution to the Lambert W function can be computed using the Lagrange inversion theorem (see [16]). All other branches have been computed as well with coefficients expressed in terms of the unsigned Stirling numbers of the first kind [1]. This expansion to compute the different roots for the Lambert W function was presented by de Bruijn [23], with Comtet [17] rewriting de Bruijn's expansion in terms of the Stirling numbers and Corless et al. [20] correcting a typographic error in equation (2.4.4) in Comtet [17]. Lastly, we note the lambertw function in Matlab, Lam-bertW in Maple, and ProductLog in Mathematica allow for the numerically accurate computation of W j (x) for any branch j. Since Asl and Ulsoy's original article, the connection between DDE characteristic roots and solutions to Lambert W has been substantially extended and expanded and we direct the interested reader to the monograph by Yi, Nelson, and Ulsoy [65]. These results, however, do not extend to multi-lag DDEs as (to the best of our knowledge) the nonlinear characteristic equation can be cast as a Lambert W equation only when the DDE has a single lag. Several researchers have studied the two-lag case [11,15,32,42,46,50,54,57,60,62,63], including those that focus specifically on diabetes [38,41]. However, the majority of these efforts have focused on the geometry of the stability regions and the specific value of the delay at the Hopf bifurcation that commonly occurs as the delay is increased.
In this work, we develop a novel series expansion for all the roots to the nonlinear characteristic equation corresponding to a two-lag DDE. While we were inspired by the original derivation of (5), the development presented below is a bit more involved than the Lambert W derivation in de Bruijn [23] partially due to our goal of making the truncated sums computationally efficient to evaluate.
We note that our asymptotic expansion approach here differs from much of the literature in that we are interested in computing the roots regardless of the delay value. The other approaches, including D-composition and τ -composition [54,61], spectral elements [37], and semi-discretization [44] focus on identifying the value of the delay at the Hopf bifurcation. We do note that some of these methods have been applied to specifically study two-lag DDEs [54,57,63]. However, the result in those articles involve collections of stability criteria and do not generate the complete set of roots. In our future efforts we plan to pursue a rigorous comparison with these existing methods.
Section 2 presents the main result of the paper, deriving the formula 2 for the roots. Section 3 provides numerical results supporting our claims of accuracy and convergence and including an application to a two-lag model from ecology. Lastly, Section 4 summarizes the results of the work and discusses the natural next steps for future work.
2. Two-lag development. In this section, we develop our series expansion for the roots to the nonlinear characteristic equation arising from the two-lag DDE. Note that for the multi-delay case, the conventional approach focuses on developing geometric bounds for the roots (see [32] for a specific example and [10, §12] for a good summary and overview), whereas here, we are explicitly computing the roots.
Consider a DDE with two lags τ 1 and τ 2 with corresponding coefficients α, β, and γ. We will rescale and rearrange terms, defining several new variables and parameters in order to reduce the equation to its essential mathematical features. We rescale time by the first lag t = t/τ 1 , the parameters so that After a Laplace transform, the corresponding nonlinear characteristic equation for s ∈ C is thus and multiply through by e λτ , yielding the much cleaner equation Inspired by the original expansion in [23], we note that for |γ 2 | near to or far from zero, we can write λ as a small perturbation u of ln γ 2 λ = (ln γ 2 + u)/τ , and then safely assume that |u| |ln γ 2 | . (8) Substituting this form of λ into equation (7) we find that To proceed, we will also assume that in addition to (8). This restriction on the coefficients is a result of the chosen rearrangement of (6) into (7). There is nothing particularly special about (10) and choosing to recast (6) differently would simply necessitate different assumptions on the parameters in order for an expansion to converge. In our experience, this is a reasonable assumption encountered in several models (see the examples in Section 3).
Moving forward, we can solve (9) for u as which means that u = ln τ + ln( 1 /ln γ2) + v (11) where our attention now focuses on solving for a new variable v. As discussed by Corless et al. [20], the nested logarithms in (11) need not select the same branch.
As Corless et al. do, we will let the outer logarithm be the principal branch and denote Ln as the inner logarithm for which we have not yet chosen a branch, i.e., By substituting u back into (9), we can cancel several terms and simplify our problem to one of finding the roots v to the equation then our efforts will focus on identifying the roots of Note that the motivation for adding and subtracting c will become self-evident in the proof for the following lemma.
Lemma 2.1. If there exists positive numbers a and b such that |σ| < a and |µ| < a then (12) has a single root inside the region |v| ≤ b.
By Cauchy's theorem, we can then compute v using the well-known formula from complex analysis for b = π min{1, τ }. The 1/f (ζ) can be expanded into the following (absolutely and uniformly convergent) power series with some coefficients a m . The substitution of (14) into (13) and evaluation of the contour integral generates a double series in m and k. Note, however, that the m = 0 term is zero because the integrand has the same number of roots (at ζ = 0) in the numerator and the denominator. The doubly infinite and absolutely convergent power series for v in σ and µ is thus where the a km are coefficients independent of σ and µ. The coefficients a km are computed by first computing a power series expansion in σ and then applying the Lagrange Inversion Theorem [1] to obtain v as a function of µ. To facilitate rearrangement of this series into the form of (15), we denote and note that the derivative of f (w) with respect to σ is −w. Equation (15) can thus be written (after the power series expansion in σ around zero) as where mk is a rising factorial 3 , and h m,k = lim w→0 h m,k (w) with

Straightforward application of the product rule yields that
The lth derivative of f 2 (w) −(m+k) is computed using Faà di Bruno's formula [36] for higher derivatives of composite functions and thus and B l,p are the partial Bell polynomials [2,8,9,64]. Recall that for an arbitrary sequence {x i } and l, p ∈ N, the partial Bell polynomials are In what follows, the argument for the Bell polynomial frequently takes the form {f To improve the clarity of the formulas, we denote ; n ∈ N, to be a finite sequence of higher derivatives of f 2 .
Formally taking the limit for lim w→0 h m,k (w) yields The evaluation of most terms in (17) is computationally straightforward and efficient, except the term involving B (q) l,p . For example, the sum in (16) for m = 10 and k = 100 can take around one minute on an Intel i7-2620M cpu when using Mathematica to analytically find the qth derivative of B l,p and then sum the terms. In [48], Mishkov derived a formula for B (q) l,p , but to implement it would consists of solving a cascade of linear Diophantine equations 4 to obtain the correct indices for the sums.
We will use a recursive property of B (k) l,p to avoid having to solve systems of Diophantine equations. Note that for an arbitrary sequence and thus in the context of (17) we have that recursively allows (eventual) computation of the derivatives by evaluation of the partial Bell polynomials themselves. 5 The use of Mathematica's memoization features for recursive sums substantially sped up the computation as well, naturally at the expense of increased memory usage. 6 We now have that the principal branch of the roots of (7) is where We note that for i ≥ 1 f . and thus h m,k is defined as where and recalling that B (q) l,p is defined recursively in (18). The full description of the jth branch root for the characteristic equation in (6) is therefore where ln j (z) = |z| + arg(z) + 2πij, v is defined in (20), and the subscript on the v indicates the chosen branch for the logarithms in σ, c, and µ.
For the convenience of the reader, the Mathematica code is available in the github repository for the Mathematical Biology Group at the University of Colorado (https://github.com/MathBioCU). This code solves for s j as a function of α, β, γ, τ 1 , τ 2 , j, m, and k. Lastly, upon request, the author will supply versions of the code with arguments α 1 , β 2 , γ 2 , τ , j, m, and k as well as versions which obtain significant speedups by using Mathematica's memoization, compilation, and parallelization features.

Numerical results.
In what follows, we denoteŝ j as the actual root. Where relevant, the s j values were computed using the damped Newton's method implemented in Mathematica's FindRoot applied to equation (6).
3.1. Convergence. In this section, we provide numerical results for the estimation of the roots. Consider the two lag DDĖ x(t) = φ(t) ; t ≤ 0 , system from above with transfer function We note that the accuracy of our asymptotic expansion rests on the assumption that (10) is small and so we arbitrarily choose which means that the parameters in (22) are The parameters for the two-lag DDE were those from (25) .
In the illustrating computations of s j below, we choose to truncate the series at m = 8 and k = 1000 as this provide a good balance between accuracy (as measured by comparison with the results from performing a root finding algorithm) and computation time. 7 We obtain the following roots as depicted in Figure 1. For comparison, Figure 2 depicts a contour plot of the transfer function in (23) (lighter colors are larger values) and we note that our computed roots match exactly to eight decimals with the singularities of the transfer function.

3.2.
Example. In modeling the life cycle of the Australian blowfly Lucila cuprina, Braddock and van den Driessche [13] developed the following two-lag DDE The positive equilibrium is at x * = 1/(a 1 + a 2 ) and a linearization about this point (shifting to y = x − x * ) yields the equation Theorem 6 in Ruan [53] summarizes the results of [13] nicely, providing 3 different sufficiency conditions for the existence of a value of τ 2 which will incur a Hopf bifurcation. We present here only the last of the possible sufficiency conditions in the following theorem from [53], since these restrictions also satisfy our conditions for accurate series expansion.
Theorem 3.1. If a 1 = a 2 and τ 1 > (2x * ra 1 ) −1 then there is a τ 0 2 > 0 such that when τ 2 = τ 0 2 the two-lag DDE (26) exhibits a Hopf bifurcation at the equilibrium x * = 1/(a 1 + a 2 ). The horizontal axis is the real part of s and the vertical axis is the imaginary part of s. Note how the peaks (lighter regions) coincide with the computed roots in Figure 1. The parameters for the twolag DDE were those from (25).
To illustrate the match between our expansion and this theory, let r = 1, a 1 = a 2 = 1, τ 1 = 10. Figure 3 depicts the principal root as a function of τ 2 and we note the root at τ 2 = 0.379414. We also note that the relation which must be small (10) is on the order of 10 −17 for this set of parameters. Figure 4 depicts asymptotically stable oscillations (solid curve for τ 2 = 0.2) and unstable oscillations (dashed curve for τ 2 = 1.3). To better illustrate the stable and unstable behaviors, we arbitrarily choose an initial history of φ(t) = 1 ; t ≤ 0 for the stable solution and an initial history of φ(t) = 1 /10 ; t ≤ 0 for the unstable solution.
4. Conclusion and discussion. In this work we have developed an asymptotic expansion for the roots to a transcendental characteristic equation associated with two-lag DDEs. We have provided numerical evidence supporting the accuracy of our expansion. We have also applied the expansion to a two-lag DDE model of the Australian blowfly, illustrating the Hopf bifurcation transition from stable to unstable oscillations.
There are many intriguing directions one could take with this work. The natural first steps are to provide a rigorous investigation of truncation criteria for the series computation as well as a comparison with other published approaches.   3) solutions to the linearized two-lag blowfly model in (27). For the stable solution, the initial history is φ(t) = 1 for t ≤ 0 and for the unstable solution, the initial history is φ(t) = 0.1 for t ≤ 0.

DAVID M. BORTZ
Concerning systems of single-lag DDEs, it is known that the approach presented in Asl and Ulsoy's original article [3] is only valid when the coefficient matrices in front of the state and lagged state commute. This limitation has been discussed in the literature [4,67] with Jarlebring [35] also noting that the matrices need only be triangularizable. Fortunately, there is an algorithm for working around this problem and the monograph [65] has a very clear description of how to solve single-lag DDE systems using a matrix Lambert W function. And a reasonable extension would be to develop the framework to identify characteristic roots for two-lag systems of DDEs. Indeed, we have made preliminary efforts in this direction as it is a natural extension to some of the author's previous work [7].
Another direction for future work would be to consider more general multi-lag DDEs. We were surprised to discover that the restrictions placed on the coefficients in the two-lag DDE were satisfied for many of the DDEs we considered. We therefore plan to investigate if the coefficient restrictions for an asymptotic expansion in the three-lag (and higher) DDEs are as generous.

5.
Disclaimer. The views expressed are those of the author and do not reflect the official policy or position of the Department of Defense or the U.S. Government. This is in accordance with DoDI 5230.29, January 8, 2009.