Optimal linear stability condition for scalar differential equations with distributed delay

Linear scalar differential equations with distributed delays appear in the study of the local stability of nonlinear differential equations with feedback, which are common in biology and physics. Negative feedback loops tend to promote oscillations around steady states, and their stability depends on the particular shape of the delay distribution. Since in applications the mean delay is often the only reliable information available about the distribution, it is desirable to find conditions for stability that are independent from the shape of the distribution. We show here that for a given mean delay, the linear equation with distributed delay is asymptotically stable if the associated differential equation with a discrete delay is asymptotically stable. We illustrate this criterion on a compartment model of hematopoietic cell dynamics to obtain sufficient conditions for stability.


Introduction
Models of self-regulating systems often include discrete delays in the feedback loop to account for the finite time required to perform essential steps before the loop is closed. Such mathematical simplifications are especially welcome in biological applications, where knowledge about the loop steps is usually sparse. This includes maturation and growth times needed to reach reproductive age in a population [31,40], signal propagation along neuronal axons [17], and post-translational protein modifications [15,44]. Introduction of a discrete delay in an ordinary differential equation can destabilize steady states and generate complex dynamics, from limit cycles to chaos [36]. Although the linear stability properties of scalar equations with single discrete delays are fairly well characterized, lumping intermediate steps into a delayed term can produce broad and atypical delay distributions that deviate from discrete delays, and it is still not clear how that affects the stability of the equation [18].
is a model paradigm in biology and physics [1,5,25,42,44,46]. The first argument of F is the instantaneous part of the loop and the second one, the delayed or retarded part, which closes the feedback loop. The integral is taken in the Riemann-Stieltjes sense. The function η is a cumulative probability distribution function, it can be continuous, discrete, or a mixture of continuous and discrete elements. In most cases, the stability of the above equation is related to its linearized equation about one of its steady statesx,ẋ where the constants a and b ∈ R are the negatives of the derivatives of the instantaneous and the delayed parts of F at x =x, Eq. (1) is also called a linear retarded functional differential equation. Basic theory for delay differential equations and functional differential equations can be found in [9] and [28]. Additional applications can be found in previously mentioned references and in [24,36].
Stability analysis of Eq. (1), when the distribution function η differs from the Dirac distribution, has been the subject of several works. In 1989, Boese [16] analyzed the stability of (1) for a Gamma distribution, and determined rather technical sufficient conditions for its asymptotic stability. Kuang [37], in 1994, considered a system of two differential equations with continuous distributed delay, possibly infinite. He focused on the existence of pure imaginary eigenvalues, and determined conditions for their nonexistence, obtaining sufficient conditions for the asymptotic stability of his system. In 2001, Bernard et al [13] considered (1) and determined sufficient conditions for its stability, mainly in the case where the distribution is symmetric about its mean. They then conjectured that the single Dirac measure would be the most destabilizing distribution of delays for (1). Atay [6] recently gave arguments in that direction. He focused on the stability of delay differential equations near a Hopf bifurcation, and for linear delay differential equations, such as (1), he showed that if the delay has a destabilizing effect, then the discrete delay is locally the most destabilizing delay distribution.
Huang and Vandewalle [30] and Tang [50] also analyzed the stability of equations similar to (1). The first authors were interested in the numerical stability of differential equations with distributed delay, but they proposed an interesting geometrical approach to determine conditions for the stability of (1) for a special delay distribution. Unfortunately, their method cannot be generalized to general distributions. In [50], Tang determined sufficient stability conditions for very general differential equations with distributed delay, but his results are very technical and not easy to handle in practice. Adimy et al [1] and Crauste [23] obtained sufficient conditions for the existence of a Hopf bifurcation when the delay density function is decreasing. In [45], Ozbay et al. investigated the stability of linear systems of equations with distributed delays, and applied their results to a model of hematopoietic stem cell dynamics. Considering an exponential distribution of delays, they obtained necessary and sufficient conditions for the stability using the small gain theorem and Nyquist stability criterion. Solomon and Fridman, using Linear Matrix Inequalities, also established sufficient conditions for exponential stability of systems with infinite distributed delays [47]. Berezansky and Braverman recently obtained sufficient conditions for the stability of non-autonomous differential equations with distributed delay [11,12].
Finally, let us mention the work of Anderson [2,3], who focused on the stability of some delay differential equations, called regulator models, which are a particular form of (1). The theory developed by Anderson [2,3] focuses on the properties of the probability distribution η. Although the results of Anderson are only valid for some class of probability measures, they stress the importance of the shape of the delay distribution. Moreover, Anderson mentions that "the more concentrated the probability measure, the worse the stability property of the model" [3].
Although it has been observed that in general a greater relative variance provides a greater stability, a property linked to geometrical features of the delay distribution [2], there are counter-examples to this principle. Yet, as mentioned above, it has been conjectured that among distributions with a given mean, the discrete delay is the least stable one [6,13]. If this were true, a theorem due to Hayes [29] would provide a sufficient condition for the stability of the trivial solution of delay differential equations independently from the shape of the delay distribution. This conjecture has been proved by Krisztin using Lyapunov-Razumikhin functions when there is no instantaneous part [35], and by different authors for distributions that are symmetric about their means [6,13,33,43]. It is possible to lump the non-delayed term into the delay distribution and use the condition found in [35], but the resulting stability condition is not optimal. Here we prove that the conjecture is true for all delay distributions with exponential tails. That is, for a given mean delay, the scalar linear differential equation with a distributed delay is asymptotically stable provided that the corresponding equation with a single discrete delay is asymptotically stable. This sufficient condition for stability is optimal in the sense that if it is not satisfied, we can find a distribution with distributed delay for which the equation is not stable. To illustrate this general result, we consider a compartment model of hematopoiesis that can be expressed as a scalar differential equation with an arbitrarily complex delay distribution, and we obtain a simple stability condition.
In section 2, we provide definitions and, in Section 3, we set the stage for the main stability results. In section 4, we show that a distribution of discrete delays is necessarily stable when the discrete distribution with a single delay equal to the mean is stable. In section 5, we present the generalization to any distribution, hence showing that distributions with distributed delays provide more stability than the discrete distribution with the same mean. Section 6 is devoted to the presentation of a model for hematopoiesis and the illustration of the stability problem.
The notation x t (σ, φ) is also used to refer to the solution of Eq. (2) associated with the initial conditions σ and φ.
is called the characteristic function of the linear equation (2). The equation D(λ) = 0 is called the characteristic equation of (2).
The following theorem [26,48] gives a necessary and sufficient condition for the asymptotic stability of x = 0.
Theorem 2.4. Suppose that there exists ν > 0 such that the following inequality is satisfied: The solution x = 0 of Eq.  Inequality (4) implies that the mean delay value is finite, We assume in the following that inequality (4) is always satisfied. For more details concerning retarded functional differential equations with infinite delays, see [26,27]. When η represents a single discrete delay (η a heaviside function), the asymptotic stability of the zero solution of Eq. (2) is fully determined by the following theorem, originally due to Hayes [29].
More generally, the following statements always hold for any delay distribution: (i) When a ≤ −b, the characteristic equation of Eq.
(2) has a positive real root.
(ii) When a ≥ |b| and a > −b, the characteristic equation of Eq.
(2) has no root with positive real part.
Therefore, the stability of the solution x = 0 depends on the delay distribution only in the parameter space region b > |a| and, from now on, we restrict the stability analysis to that region. Assuming b > 0 and making the change of timescale t → bt, we have a → a/b, b → 1 and η(τ ) → η(bτ ). Eq. (2) can be rewritten aṡ The delay distributions affect the stability of Eq. (5) when a ∈]0, 1[. The characteristic equation is called stable if all roots have (λ) < 0 [48]. To emphasize the relation between the stability and the delay distribution, we give a similar definition for the delay distribution.
The integral term in Eq. (6) is the Laplace transform L of the distribution η. Along the imaginary axis λ = iω, the Laplace transform can be expressed as (Lη)(iω) = C(ω) − iS(ω), where The strategy for determining the stability of distributed delays is the following. We use a geometric argument to bound the roots of characteric equation (6) by the roots of the characteristic equation for a single discrete delay. More precisely, we will show that if the leading roots associated to the discrete delay are a pair of imaginary roots, then all the roots associated to the distribution of delays have negative real parts. We first state, in Section 3, a criterion for stability: if S(ω) < ω whenever C(ω) = −a, then the distribution is stable. We then show in Theorem 4.5 that a distribution of n discrete delays is more stable than a certain distribution with two delays (in the sense that S(ω) ≤ S * (ω), where the distribution with n delays is denoted by η and the "special" distribution with two delays by η * ). We construct this most "unstable" distribution and determine that only one of the delays is positive, so that its stability can be determined using Theorem 2.5. We then generalize for any distribution of delays in Section 5.

General Stability Criteria
Assume a ∈] − 1, 1[, and let η be a distribution with mean E. We consider the family of distributions, scaled with the parameter ρ ≥ 0, where H(τ ) is the step or heaviside function at 0, corresponding to a single discrete delay vanishing at τ = 0. The distribution η ρ has a mean ρE ≥ 0. The notation D ρ is used to refer to the characteristic equation associated with the scaled distribution η ρ . The characteristic equation for the distribution η 0 is D 0 (λ) := λ + a + 1 = 0.
The next proposition provides a necessary condition for instability. It is a direct consequence of Theorem 2.19 in [48].
Proof. Suppose that the distribution η is unstable, i.e. that the characteristic equation has roots λ with (λ) ≥ 0. Consider the family of scaled distributions η ρ . The roots of the characteristic equation D ρ = 0 depend continuously on the parameter ρ and roots with positive real parts can only appear by crossing the imaginary axis. The scaled distribution η ρ is stable for ρ = 0 (the only root is λ = −(a + 1) < 0) and unstable for ρ = 1. Hence there exists a critical value 0 < ρ ≤ 1 at which η ρ loses its stability, and this happens when the characteristic equation D ρ (λ) = 0 has a pair of imaginary roots λ = ±iω, with ω ≥ 0. Splitting the characteristic equation in real and imaginary parts, we have Since −ω satisfies the above system, we only look from now on and throughout this manuscript to positive values of ω. The upper bound on ω, ω c = √ 1 − a 2 , is obtained by applying Cauchy-Schwartz inequality, Finally, setting ω s := ρω, we obtain 0 < ω s ≤ ω ≤ ω c and This completes the proof.
Proposition 1 provides a sufficient condition for asymptotic stability, stated in the following corollary. Corollary 1. The distribution η is stable if one of the two following conditions is satisfied: The condition S(ω) < ω is not necessary for stability, as there are cases where S(ω) ≥ ω even though the distribution is stable. This happens when an unstable distribution switches back to stability as E is further increased (see [10] or [16]).

Stability of a distribution of discrete delays
In this section, we show that a distribution with n discrete delays and mean E is more stable than the distribution with a single discrete delay E. It is convenient to represent distributions of discrete delays by their densities. We denote a density of n discrete delays τ i ≥ 0, and weights p i > 0, i = 1, ..., n, n ≥ 1, as Following Corollary 1, for f n to be stable, it is enough to show that S n (ω s ) < ω s whenever C n (ω s ) = −a, ω s ≤ ω c . We now show that among all distributions satisfying C n (ω s ) = −a for a fixed ω s , there exists a density f * that maximizes the values of S n (ω s ). This density f * has only one positive delay, making it easy to show that S * (ω s ) < ω s . This would imply that all discrete delay distributions are stable.
Definition 4.1. We define the constants c ≈ 0.7246 and θ c ≈ 2.3311, where c is the smallest positive value such that cos(θ) ≥ 1 − cθ for all θ > 0, found by solving the two equations c = sin(θ) and 1 − θ sin(θ) = cos(θ) for c > 0, θ > 0, and θ c is the positive value for which cos(θ) = 1 − cθ. We define the convex function The following lemmas show how to find the distribution that maximizes S n (ω s ) for n = 2.
with ω c = √ 1 − a 2 . Suppose that there exists ω s ∈ [0, ω c ] and a density f 2 with mean E, such that Then ω s E < θ c and cos(ω s E) > −a.
Consequently, cos(ω s E) > −a and, using (11), we then deduce that cos(ω s E) > C 2 (ω s ). Furthermore, we have The first inequality comes from the definitions of c and g (see Definition 4.1): The second inequality is the convexity property of g. Thus, we deduce cos(ω s E) > g(ω s E). Since g(x) = cos(x) for x ≥ θ c , this means that ω s E < θ c .
and a density f 2 with mean E, such that equality (11) is satisfied. Then there exists a unique density f * with two discrete delays τ * 1 and τ * 2 , mean E, such that τ * Proof. Suppose there exists a density f 2 with two discrete delays τ 1 and τ 2 , weights p 1 and p 2 , mean E, satisfying τ 1 = 0 and τ 2 > 0. Necessarily, p 2 τ 2 = E (so f 2 has mean E). We are going to show that By using p 1 = 1 − p 2 and p 2 = E/τ 2 , Eq. (13) is equivalent to To see that inequality (15) is indeed satisfied, one can note that, using (11), (15) holds true. Consequently Eq. (14) has at least one solution satisfying 0 ≤ ω s τ 2 ≤ π. Moreover, since θ c is a tangent point (see Definition 4.1), there is exactly one solution satisfying 0 ≤ ω s τ 2 < θ c < π. Denote by τ * 2 the smallest value of τ 2 that solves Eq. (14), and define From the definition of τ * 2 , f * exists and is unique. It remains to show that f * is a well-defined density, that is p * 2 ∈ [0, 1]. Since τ * 2 is the smallest and unique positive solution in the interval [0, θ c ] of (14), the sign of cos(x) − (1 − (1 + a)x/(ω s E)) determines whether x is smaller or larger than τ * 2 in the interval [0, θ c ]. From Lemma 4.2, ω s E < θ c and cos(ω s E) > −a, or formulated equivalently, cos(ω s E) we obtain the result 0 < p * 2 < 1, which shows that f * is a well-defined density.
Lemma 4.4. Assume a ∈ ]−1, 1[ and E > 0 satisfies (10). Suppose that there exists ω s ∈ [0, ω c ] and a density f 2 with mean E, such that equality (11) is satisfied. Then for any density f 2 with mean E and satisfying Eq. (11), we have where the density f * is defined in Lemma 4.3.
Proof. We recast the problem in a slightly different way. Consider a density with two discrete delays τ 1 and τ 2 and mean E, such that C 2 (ω s ) = −a. Writing u = ω s τ 1 , v = ω s τ 2 and T = ω s E, we can express the weights p i in terms of u and v: By convention, 0 ≤ u < T < v. We consider C 2 (ω s ) and S 2 (ω s ) as functions of u and v; C, S : The subscripts 2 have been dropped to ease the reading. Equation (17) Since 0 ≤ u < T < π, cos(u) is decreasing, cos(u) > cos(T ), and cos(T ) > −a (From Lemma 4.2, we know that cos(ω s E) > −a), so cos(u) + a > 0. Eq. (18) writes ) and α(u) is increasing for u ∈ [0, T ). The slope of the right hand side of (19) is negative, cos(v i ) is decreasing with solutions v i of (19) ( Figure 1A). One may note that the points C(u, v i ), S(u, v i ) are at the intersection of the chord i between the unit circle points cos(u), sin(u) and cos(v i ), sin(v i ) and the vertical secant at −a. From (16) with C(u, v i ) = −a, it is easy to see that cos(v i ) < −a since cos(u) > −a. By displaying the above mentioned chords and the vertical secant on a unit circle ( Figure 1B v 1 (0) = ω s τ * 2 . Therefore, we need to show that S(0, v 1 (0)) maximizes S(u, v 1 (u)). The total derivative of S with respect to u is If ∂S/∂v < 0, the total derivative is strictly negative if and only if The partial derivative with respect to v is One can see that v 1 always satisfies v 1 (u) ≤ π. Indeed, if one assumes by contradiction v 1 (u) > π, then first Eq. (19) has no root on the interval [T, π], and second, since It follows that for all v > π, and Eq. (19) has no root, yielding a contradiction. The sine function is strictly concave on the interval [0, π] and this implies that This shows that ∂S/∂v < 0. Now, Inequality (20) can be re-expressed as It can be verified that this inequality is satisfied for v −u = z ∈ (0, π]. The left-hand side vanishes when z → 0, and the derivative is strictly positive for 0 < z ≤ π: The last inequality is obtained with inequality (21). Therefore, dS/du < 0 and S is maximized for u = ω s τ * 1 = 0 and v 1 (0) = ω s τ * 2 < π. Now that we established the existence of a density f * with two delays, one equal to zero the other one positive, and mean E which maximizes the quantity S 2 (ω s ), we prove in the next theorem the stability of all densities with n discrete delays and mean E satisfying (10).
Theorem 4.5. Assume a ∈ ]−1, 1[ and E > 0 satisfies inequality (10). Let f n be a discrete density with n ≥ 1 delays and mean E, then the density f n is stable.
In the second step, we reduce the number of strictly positive delays. All pairs of delay τ i < τ j for which the inequality holds are iteratively replaced by one positive and one vanishing delay, as done in Lemma 4.3. We note that inequality (22) reduces to C(ω s ) = p 1 cos(ω s τ 1 ) + p 2 cos(ω s τ 2 ) ≤ cos(ω s E) for a two discrete delay distribution, with delays τ 1 and τ 2 satisfying (9). This transformation preserves the values of mean E and C (ω s ), and increases the value of S (ω s ). This step is repeated until one of the two situations occurs: (i) There remains one density f * with exactly one delay τ * 1 = 0 and one delay τ * 2 > 0. Then the inequality S * (ω s ) < ω s follows from the first part of the proof. Therefore, S(ω s ) ≤ S (ω s ) ≤ S * (ω s ) < ω s , and, by Corollary 1 implies that f is stable. (ii) There remains a densityf with one delayτ 1 = 0 and two or more delaysτ k > 0, k = 2, . . . , m, m ≥ 3, such that for each pair i = j ∈ 2, . . . , m. Since Because the sine function is concave on the interval [0, π], any averaging of delays can only increase the value of S. We now have a density f with τ 1 = 0 and τ 2 > 0, p 1 =p 1 and p 2 = 1 −p 1 , C (ω s ) <C(ω s ) (from inequality (23)), E =Ē ≤ E, and S (ω s ) ≥S(ω s ). We now replace τ 2 with a delay τ * 2 < τ 2 , so as to obtain a density f * with C * (ω s ) =C(ω s ) = −a, and E * = E .
Consequently, this last change of delay has the effect of increasing the value S * (ω s ) ≥ S (ω s ), while maintaining the condition C * (ω s ) = −a. Since the mean E * of density f * satisfies inequality (10), we have S * (ω s ) < ω s as shown for the case n = 2. Therefore S(ω s ) ≤ S (ω s ) ≤S(ω s ) ≤ S (ω s ) ≤ S * (ω s ) < ω s . Corollary 1 implies that f is stable.

Stability of a general distribution of delays
We now show that the stability of discrete delays implies the stability of general distributions. First we need to bound the roots of the characteristic equation for general distributed delays.
Lemma 5.1. Assume a ∈ ]−1, 1[ and E > 0 satisfies inequality (10). Let η be a delay distribution with mean E and characteristic equation D(λ) = 0. There exists a sequence of distributions {η n } n≥1 with mean E, such that η n converges weakly to η as n → ∞, and λ is a root of the characteristic equation if and only if there exists a sequence of characteristic roots λ n for η n such that lim n→∞ λ n = λ. If {µ n } n≥1 is a sequence of real parts of characteristic roots λ n for η n , D n (λ n ) = 0, then lim sup n→∞ µ n < 0.
Proof. Existence of a sequence {η n } n≥1 of distributions with n delays and mean E, such that η n converges weakly to η as n → ∞ is rather straightforward, this sequence can be built explicitly. We do not detail this part here.
Consider λ n = µ n + iω n a root of the characterisitic equation for η n . The mean E satisfies inequality (10), so µ n < 0. Then, as n → ∞ by weak convergence. Thus any converging sub-sequence of roots converges to a root for η. The same way, if λ is a root for η, as n → ∞. Convergence is guaranteed by inequality (4). Thus each root λ n lies close to a corresponding root λ, and µ = lim sup n→∞ µ n , with µ n real part of a characteristic root λ n , is the real part of a characteristic root for η. Since µ n < 0, we have that µ is non-positive. Suppose µ = 0 and consider the scaled distribution η a,ρ (τ ) defined by (7), and the associated real parts µ a,ρ , where the subscript a is there to emphasize the dependence of the stability on the parameter a in the characteristic equation. Then, by continuity, there exists (ā, ρ) in an εneighborhood of the point (a, 1) for which µā ,ρ > 0. For sufficiently small ε > 0, inequality (10) is still satisfied: Additionally, the scaled discrete distributions η n,ā,ρ converge weakly to ηā ,ρ , so that the real parts µ n,ā,ρ of the roots converging to µā ,ρ become eventually positive. That is, there is N > 1 such that η n,ā,ρ is unstable for all n > N , a contradiction to Theorem 4.5, since inequality (10) still holds. Therefore µ < 0.
Proof. Consider a sequence of distributions with n delays {η n } n≥1 where η n converges weakly to η. By Lemma 5.1, the roots of the characteristic equation of η have strictly negative real parts. Therefore η is stable.
The results obtained above provide the most complete picture of the stability of Eq. (2) when the only information about the distribution of delays is the mean. These results are summarized in the following theorem and illustrated in Fig. 2.
The zero solution of Eq. (2) may not be asymptotically stable (depending on the particular distribution) if b > |a| and The zero solution of Eq. (2) is unstable if a ≤ −b.  (3): a delay-independent stability region (light grey, (1)), delimited by the condition a ≥ |b|; a discrete-delay stability region (conditionally stable, light-grey, (2)), delimited by condition 24; and a distributed-delay-dependent stability region (white, (3)). The instability region is composed of a distributed-delay-dependent instability region (conditionally stable, white, (4)) and a delay-independent instability region (unstable, dark grey, (5)), delimited by the curve b = −a. The discrete and distributed delay stability boundaries intersect at point (a = −1/E, b = 1/E). The arrow pointing leftward shows that there exists a region, for b > 1/E, where a stable steady state can become unstable through a decrease of the value of a, independently of the shape of the delay distribution. The distributed delay is f (τ ) = 0.8δ(τ − 0.625) + 0.2δ(τ − 3.5), with mean delay E = 1.2 (parameters as in Figure 1).

Compartment Model of Hematopoiesis
Circulating blood cells are continuously renewed by a hierarchical structure of cells maintained by hematopoietic stem cells (HSCs). Hematopoiesis consists in a complex set of feedback loops that control blood cell production. HSCs can either self-renew or differentiate to one of the three main blood cell lineages: white blood cells, platelets and red blood cells. Through successive division and differentiation stages, HSCs become progenitors (immature cells), precursors (differentiated cells), and then fully mature cells. At every stage of this hierarchy, feedback loops regulate cell differentiation, proliferation, and death. The process of red blood cell production is tightly controlled by erythropoietin, a growth factor released by the kidneys when blood oxygen is low, and whose action inhibits cell death [34]. Platelet production and white blood cell production processes are also controlled by growth factors (thrombopoietin [32] and G-CSF [7], respectively). It is usually thought that mature blood cells act negatively, through growth factor release, on precursors, progenitors and HSCs dynamics [19,20].
From a modeling viewpoint, the hierarchical structure of hematopoiesis can be described by a finite system of differential equations, each equation describing the dynamics of one cell generation [8,14,19,20,40,49]. Such a view is largely accepted, both by modelers and biologists, even though mechanisms involved in cell differentiation processes are complex and there is no reason to believe that cells always go through a forward differentiation process.
In 2005, Colijn and Mackey [19,20] proposed a compartment model of hematopoiesis, based on previous models of hematopoietic stem cell dynamics [40], white blood cell dynamics [14], platelet dynamics [4] and red blood cell dynamics [8]. This model consists in a system of 4 differential equations with discrete delays. Each equation describes the number of either HSCs, red blood cells, white cells or platelets. Cells spend a finite amount of time in each of these compartments during which they mature and divide. Delays account for cell stage durations. Colijn and Mackey's model [19,20] has been further justified and numerically analyzed by Colijn and Mackey [21] and Lei and Mackey [38], who showed that it exhibits multiple steady states. Stability analysis of this model is made difficult by the presence of several discrete delays. A simpler model, based on ordinary differential equations, can then be considered, similar to the one by Stiehl and Marciniak-Czochra [49]. However, even in this case, the structure of the system with several compartments induces a natural delay, and the stability analysis is not straightforward.
We consider a compartment model of hematopoiesis that encompasses the main dynamical properties existing hematopoiesis models, and focus on stability conditions for this system. The compartment model can be expressed as a single equation with a general distributed delay. We showed that among all delay distributions with a given fixed mean, the distribution with a single discrete delay (that is, the delay equals the mean) is the most unstable one. Consequently we can provide a condition for the stability of the hematopoiesis model by determining when the equation with a single delay is stable.
Let denote by x(t) the number of HSCs at time t, and by z i (t), i = 1, 2, 3, the densities of circulating platelets, white cells, and red blood cells, respectively. We assume that x produces the quantities z i through a linear chain process, describing the compartmental structure of each hematopoietic lineage. The number of mature cells z i act on a negative feedback loop that represses the production of x. The disappearance rate of HSCs, α, is assumed constant. The HSC production rate P is a function that depends on x and a weighted average z of the repressors z i . Namely The HSC number x is governed by the equation Each mature cell number z i (t), i = 1, 2, 3, is assumed to be the product of a linear chain of differential equations of the type In the i-th hematopoietic lineage, the cell number in generation j-th is denoted by y This situation hypothesizes that each compartment in each hematopoietic lineage depends only on the previous compartment and, except for the source term β i x, lineages are independent from each other. This system is an instance of a nonlinear system with a linear subsystem [22,39]. For each lineage i, thanks to the usual chain trick in System (26), the repressors z i can be expressed in terms of the history of x convoluted by a Gamma distribution, When one focuses only on one hematopoietic lineage, and z = z i (p j = 0 for j = i), Eq. (25) can be expressed as a distributed delay equation with a Gamma distribution with mean E i = q i /β i and variance V i = q i /β 2 i . Two limiting cases are useful to consider. When q i = 1, mature cells are produced directly from HSCs, and the Gamma distribution becomes an exponential distribution with parameter β i . When E i = q i /β i is made constant and q i → ∞, the Gamma distribution converges to a Dirac mass at E i .
In addition to these three standard delay distributions, more general delay distributions are obtained by considering the above-mentioned linear parallel negative feedback loops. From System (26), the weighted repressor z(t) remains a delayed version of x(t), where the density of the distributed delay is a weighted average of Gamma densities, The delay has a mean E p = 3 i=1 p i q i /β i . In the limiting case where the length q i of each loop becomes infinite while keeping the ratio q i /β i constant, the distribution becomes a combination of discrete delays. Therefore, by a suitable choice and number of parallel negative feedback loops, one can obtain an arbitrary complex distribution of delays.
After expressing the repressor z as a function of the history of x in (27), one can then write the following equation for x, from (25) and (27), The dynamics of System (25)-(26) is entirely contained in (28). Although the production term depends continuously on the history of x, the initial conditions need only to be known at a finite number of locations. Analyzing the stability of Eq. (28) is however as difficult as the stability of the System (25)- (26).
As a nonlinear production term P , we consider the case of a mixed feedback loop, observed when a repressor (mature cells) and an activator (immature cells) are competing. The nonlinear term in equation (28) is then (29) P (x, z) = k 0 x r 1 + z h .
The parameter r is related to the degree of cooperativity of the positive loop. For r > 1 the positive loop is positively cooperative and multiple stable steady states are possible. When r = 1 the positive loop is neutrally cooperative and at most one positive steady state exists. For 0 ≤ r < 1, the positive loop is negatively cooperative and there is a single positive steady state. When r = 0, the dependence on x of the production rate P is lost. To ensure solutions are bounded, we set r ≤ h. The parameter h is the Hill coefficient describing the degree of cooperativity of the negative loop. The higher the value of h, the steeper the negative control. We assume h > 1. With these conditions, there is always at least one steady statē x ≥ 0. Eq. (28) linearized around a positive steady statex > 0 iṡ For positive cooperativity (1 < r ≤ h), there is a stable steady statex 0 = 0. In addition, there are either zero, one or two positive steady states given by the roots of the equation αx h − k 0x r−1 + α = 0. In terms of Eq. (1), a = α(1 − r) < 0 and b(x) = α 2 hx h−r+1 /k 0 > 0. The smaller positive steady statex 1 satisfies a ≤ −b(x 1 ) and, by Theorem 5.3, is always unstable. The larger steady statex 2 satisfies a > −b(x 2 ) and the sufficient condition on stability of Theorem 5.3 can be applied in the following proposition.
Proposition 2 (positive cooperativity). Assume P is given by (29) and 1 < r ≤ h (mixed feedback loop with positive cooperativity). When they exist and are distinct, the smaller positive steady statex 1 of (28) is unstable, and the larger positive steady statex 2 is linearly asymptotically stable if Whenx 1 =x 2 , the positive steady state is unstable. The zero steady statex 0 = 0 is always linearly stable.
For negative cooperativity (0 ≤ r < 1), there exists a steady statex 0 = 0 only if r > 0, in which case it is unstable. In addition, there is a unique positive steady state given by the root of the equation α(1 +x h )x 1−r = k 0 . The linear equation is given by equation (30), and the instantaneous coefficient is For neutral cooperativity (r = 1), there is a steady statex 0 = 0, whose stability depends on the existence of a positive steady state. There exists a positive steady statex = ((k 0 − α)/α) 1/h only if k 0 > α, and in this case a = 0 and b = αh(k 0 − α)/k 0 > 0. Theorem 5.3 can be applied in the following proposition to determine stability .
A model with neutral cooperativity has been considered before by Mackey and Glass [41] in the context of blood cell production. Neutral cooperativity arises when HSCs proliferate at a rate proportional to their number. In this situation, the steady state can be solved explicitly and the stability condition is relatively simple to state. The existence condition defines whether stem cells reproduce quickly enough to maintain their population (k 0 > α) or not. The original Mackey-Glass equation contained a single discrete delay at E p . Replacing the discrete delay by a general delay distribution cannot make the positive steady state unstable, as illustrated in Fig. 3.

Conclusion
We have shown that for a given mean delay, the scalar linear differential equation with a distributed delay is asymptotically stable provided that the corresponding equation with a single discrete delay is asymptotically stable. Hence, linear systems with a discrete delay are "more" unstable than linear systems with distributed delay. This result provides a sufficient condition for the stability of a large class of linear systems, as instanced by a model of hematopoiesis with parallel lineages.
Quite often the aim of the modeling is not to reproduce stability but rather instability, via periodic oscillations. Pathological cases in hematopoiesis (blood diseases, leukemias) can for instance often be explained by the destabilization of the steady state which starts oscillating periodically. Our result shows that it is more difficult to reproduce periodic oscillations, observed experimentally, with a distributed delay than with a discrete delay.  is an average of three Gamma densities (grey lines) with mean delay E = q/β = {1, 2, 3}. The discrete delay is the mean delay E p = 3 i=1 p i E i = 2 (dashed ). (B) Stability chart of the positive steady statex = 1, for varying r and h, and the stability condition is given in Proposition 3. The distributed delay is stable at points i, ii and iv, while the discrete delay is stable at points i and iv. Color coding is as in Fig. 2. (C, D) Time series of the system with distributed (solid ) or discrete delay (dashed ). (C) Neutral cooperativity, increasing Hill coefficient: r = 1 and h = 1.5 (i), 1.9 (ii) and 3.0 (iii). (D) Constant Hill coefficient, increasing cooperativity: h = 1.9 and r = 0 (iv), 1 (v) and 1.3 (vi).