Samuel Bernard

Fabien Crauste

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 [START_REF] Hutchinson | Circular causal systems in ecology[END_REF][START_REF] Mackey | Unified hypothesis of the origin of aplastic anaemia and periodic hematopoiesis[END_REF], signal propagation along neuronal axons [START_REF] Campbell | Time delays in neural systems[END_REF], and post-translational protein modifications [START_REF] Bernard | Modelling Transcriptional Feedback Loops: The Role of Gro/TLE1 in Hes1 Oscillations[END_REF][START_REF] Monk | Oscillatory expression of Hes1, p53, and NF-κB driven by transcriptional time delays[END_REF]. Introduction of a discrete delay in an ordinary differential equation can destabilize steady states and generate complex dynamics, from limit cycles to chaos [START_REF] Kuang | Delay differential equations: With applications in population dynamics[END_REF]. 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 [START_REF] Campbell | Approximating the Stability Region for a Differential Equation with a Distributed Delay[END_REF].

The delayed feedback differential equation of the form

ẋ = F x, ∞ 0 x(t -τ )dη(τ )
Work supported by ANR grant ProCell ANR-09-JCJC-0100-01.

1 is a model paradigm in biology and physics [START_REF] Adimy | A mathematical study of the hematopoiesis process with applications to chronic myelogenous leukemia[END_REF][START_REF] Atay | Distributed delays facilitate amplitude death of coupled oscillators[END_REF][START_REF] Eurich | Distributed delays stabilize ecological feedback systems[END_REF][START_REF] Meyer | Distributed delays stabilize neural feedback systems[END_REF][START_REF] Monk | Oscillatory expression of Hes1, p53, and NF-κB driven by transcriptional time delays[END_REF][START_REF] Rateitschak | Intracellular delay limits cyclic changes in gene expression[END_REF]. 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 states x, ẋ = -ax -b ∞ 0

x(t -τ )dη(τ ) [START_REF] Adimy | A mathematical study of the hematopoiesis process with applications to chronic myelogenous leukemia[END_REF] 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 [START_REF] Bellman | Differential-difference equations[END_REF] and [START_REF] Hale | Introduction to functional differential equations[END_REF]. Additional applications can be found in previously mentioned references and in [START_REF] Erneux | Applied delay differential equations[END_REF][START_REF] Kuang | Delay differential equations: With applications in population dynamics[END_REF].

Stability analysis of Eq. ( 1), when the distribution function η differs from the Dirac distribution, has been the subject of several works. In 1989, Boese [START_REF] Boese | The stability chart for the linearized cushing equation with a discrete delay and gamma-distributed delays[END_REF] analyzed the stability of (1) for a Gamma distribution, and determined rather technical sufficient conditions for its asymptotic stability. Kuang [START_REF]Nonoccurrence of stability switching in systems of differential equations with distributed delays[END_REF], 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 [START_REF] Bernard | Sufficient conditions for stability of linear differential equations with distributed delay[END_REF] considered [START_REF] Adimy | A mathematical study of the hematopoiesis process with applications to chronic myelogenous leukemia[END_REF] 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 [START_REF] Adimy | A mathematical study of the hematopoiesis process with applications to chronic myelogenous leukemia[END_REF]. Atay [START_REF]Delayed feedback control near Hopf bifurcation[END_REF] 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 [START_REF] Adimy | A mathematical study of the hematopoiesis process with applications to chronic myelogenous leukemia[END_REF], he showed that if the delay has a destabilizing effect, then the discrete delay is locally the most destabilizing delay distribution.

Huang and Vandewalle [START_REF] Huang | An analysis of delay-dependent stability for ordinary and partial differential equations with fixed and distributed delays[END_REF] and Tang [START_REF] Tang | Asymptotic behavior of a differential equation with distributed delays[END_REF] also analyzed the stability of equations similar to [START_REF] Adimy | A mathematical study of the hematopoiesis process with applications to chronic myelogenous leukemia[END_REF]. 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 [START_REF] Tang | Asymptotic behavior of a differential equation with distributed delays[END_REF], 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 [START_REF] Adimy | A mathematical study of the hematopoiesis process with applications to chronic myelogenous leukemia[END_REF] and Crauste [START_REF] Crauste | Stability and hopf bifurcation for a first-order delay differential equation with distributed delay[END_REF] obtained sufficient conditions for the existence of a Hopf bifurcation when the delay density function is decreasing. In [START_REF] Ozbay | Stability analysis of systems with distributed delays and application to hematopoietic cell maturation dynamics, in Decision and Control[END_REF], 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 [START_REF] Solomon | New stability conditions for systems with distributed delays[END_REF]. Berezansky and Braverman recently obtained sufficient conditions for the stability of non-autonomous differential equations with distributed delay [START_REF] Berezansky | Stability of linear differential equations with a distributed delay[END_REF][START_REF]Stability of equations with a distributed delay, monotone production and nonlinear mortality[END_REF].

Finally, let us mention the work of Anderson [START_REF] Anderson | Geometric and probabilistic stability criteria for delay systems[END_REF][START_REF]Intrinsic parameters and stability of differential-delay equations[END_REF], who focused on the stability of some delay differential equations, called regulator models, which are a particular form of [START_REF] Adimy | A mathematical study of the hematopoiesis process with applications to chronic myelogenous leukemia[END_REF]. The theory developed by Anderson [START_REF] Anderson | Geometric and probabilistic stability criteria for delay systems[END_REF][START_REF]Intrinsic parameters and stability of differential-delay equations[END_REF] 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" [START_REF]Intrinsic parameters and stability of differential-delay equations[END_REF].

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 [START_REF] Anderson | Geometric and probabilistic stability criteria for delay systems[END_REF], 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 [START_REF]Delayed feedback control near Hopf bifurcation[END_REF][START_REF] Bernard | Sufficient conditions for stability of linear differential equations with distributed delay[END_REF]. If this were true, a theorem due to Hayes [START_REF] Hayes | Roots of the transcendental equation associated with a certain differencedifferential equation[END_REF] 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 [START_REF] Krisztin | Stability for functional differential equations and some variational problems[END_REF], and by different authors for distributions that are symmetric about their means [START_REF]Delayed feedback control near Hopf bifurcation[END_REF][START_REF] Bernard | Sufficient conditions for stability of linear differential equations with distributed delay[END_REF][START_REF] Kiss | Stability implications of delay distribution for first-order and second-order systems[END_REF][START_REF] Miyazaki | Characteristic equation and asymptotic behavior of delay-differential equation[END_REF]. It is possible to lump the non-delayed term into the delay distribution and use the condition found in [START_REF] Krisztin | Stability for functional differential equations and some variational problems[END_REF], 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.

Definitions

We consider the linear retarded functional differential equation

ẋ = -ax -b ∞ 0 x(t -τ )dη(τ ) (2)
with real constants a and b. We assume that η is a cumulative probability distribution function: η : [0, ∞) → [0, 1] is monotone nondecreasing, right-continuous, η(τ ) = 0 for τ < 0 and η(+∞) = 1. The corresponding probability density functional f (τ ) is given by the generalized derivative dη(τ ) = f (τ )dτ . The following definitions and Theorem 2.4 follow from Stépán [START_REF] Stépán | Retarded dynamical systems: stability and characteristic functions[END_REF].

Let B be the vector space of continuous and bounded functions on

[-∞, 0] → R. With the norm ||φ|| = sup θ∈[-∞,0] |φ(θ)|, φ ∈ B, B is a Banach space.
Definition 2.1. The function x : R → R is a solution of Eq. ( 2) with the initial condition

(3) x σ = φ, σ ∈ R, φ ∈ B,
if there exists a scalar δ > 0 such that x t ≡ x(t + θ) ∈ B for θ ∈ [-∞, 0] and x satisfies Eqs. ( 2) and ( 3) for all t ∈ [σ, σ + δ).

The notation x t (σ, φ) is also used to refer to the solution of Eq. ( 2) associated with the initial conditions σ and φ.

Definition 2.2. The trivial solution x = 0 of Eq. ( 2) is stable if for every σ ∈ R and ε > 0 there exists δ = δ(ε) such that ||x t (σ, φ)|| < ǫ for any t ≥ σ and for any function φ ∈ B satisfying ||φ|| < δ. The trivial solution x = 0 is called asymptotically stable if it is stable, and for every σ ∈ R there exists ∆ = ∆(σ) such that lim t→∞ ||x t (σ, φ)|| = 0 for any φ ∈ B satisfying ||φ|| < ∆.

D(λ) = λ + a + b ∞ 0 e -λτ dη(τ ),
is called the characteristic function of the linear equation [START_REF] Anderson | Geometric and probabilistic stability criteria for delay systems[END_REF]. The equation D(λ) = 0 is called the characteristic equation of [START_REF] Anderson | Geometric and probabilistic stability criteria for delay systems[END_REF].

The following theorem [START_REF] Hale | Functional differential equations with infinite delays[END_REF][START_REF] Stépán | Retarded dynamical systems: stability and characteristic functions[END_REF] 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:

∞ 0 e ντ dη(τ ) < ∞. ( 4
)
The solution x = 0 of Eq. ( 2) is (exponentially) asymptotically stable if and only if all roots of the characteristic equation D(λ) = 0 have ℜ(λ) < 0. Theorem 2.4 is equivalent to the statement that solutions of Eq. (2) of the form x(t) = ∞ i=1 ν i (t)e λit where λ i are the roots of the characteristic equation and ν i (t) polynomials, are enough to determine the stability of x = 0. Other solutions, the small solutions, decay faster than any exponential; hence the exponential stability.

Inequality (4) implies that the mean delay value is finite,

E := ∞ 0 τ dη(τ ) < ∞.
We assume in the following that inequality (4) is always satisfied. For more details concerning retarded functional differential equations with infinite delays, see [START_REF] Hale | Functional differential equations with infinite delays[END_REF][START_REF] Hale | Phase space for retarded equations with infinite delay[END_REF]. 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 [START_REF] Hayes | Roots of the transcendental equation associated with a certain differencedifferential equation[END_REF].

E < arccos(-a/b) √ b 2 -a 2 .
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 as

ẋ = -ax - ∞ 0 x(t -τ )dη(τ ). ( 5
)
The delay distributions affect the stability of Eq. ( 5) when a ∈]0, 1[. The characteristic equation is called stable if all roots have ℜ(λ) < 0 [START_REF] Stépán | Retarded dynamical systems: stability and characteristic functions[END_REF]. To emphasize the relation between the stability and the delay distribution, we give a similar definition for the delay distribution.

Definition 2.6. The delay distribution η (or the density f ) is called stable if all roots of the characteristic equation of Eq. ( 2), or Eq. ( 5), have ℜ(λ) < 0. The delay distribution η (or the density f ) is called unstable if there exists a characteristic root with ℜ(λ) > 0.

According to Theorem 2.4 and using Definition 2.6, the solution x = 0 of Eq. ( 5) is asymptotically stable if and only if the delay distribution is stable. The characteristic equation of Eq. ( 5) is

D(λ) = λ + a + ∞ 0 e -λτ dη(τ ) = 0. ( 6
)
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

C(ω) = ∞ 0 cos(ωτ )dη(τ ), S(ω) = ∞ 0 sin(ωτ )dη(τ ).
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,

η ρ (τ ) = η(τ /ρ), ρ > 0, H(τ ), ρ = 0, (7)
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 [START_REF] Stépán | Retarded dynamical systems: stability and characteristic functions[END_REF].

Proposition 1. If the distribution η is unstable, then there exists ω s ∈ (0, ω c ], ω c = √ 1 -a 2 , such that C(ω s ) = -a and S(ω s ) ≥ ω s .
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

     ℜ(D ρ (iω)) = ∞ 0 cos(ωτ )dη ρ (τ ) + a = 0, ℑ(D ρ (iω)) = ω - ∞ 0 sin(ωτ )dη ρ (τ ) = 0. (8)
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,

a 2 + ω 2 = ∞ 0 cos(ωτ )dη ρ (τ ) 2 + ∞ 0 sin(ωτ )dη ρ (τ ) 2 ≤ 1.
Rewriting [START_REF] Bélair | Age-structured and two-delay models for erythropoiesis[END_REF] in term of η, we have

∞ 0 cos(ωρτ )dη(τ ) = -a, ∞ 0 sin(ωρτ )dη(τ ) = ω.
Finally, setting ω s := ρω, we obtain 0 < ω s ≤ ω ≤ ω c and

C(ω s ) = ∞ 0 cos(ω s τ )dη(τ ) = -a, S(ω s ) = ∞ 0 sin(ω s τ )dη(τ ) = ω ≥ ω s .
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:

(

i) C(ω) > -a for all ω ∈ [0, ω c ], (ii) C(ω) = -a, for ω ∈ ]0, ω c ], implies that S(ω) < ω.
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 [START_REF] Beretta | Geometric stability switch criteria in delay differential systems with delay dependent parameters[END_REF] or [START_REF] Boese | The stability chart for the linearized cushing equation with a discrete delay and gamma-distributed delays[END_REF]).

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

f n (τ ) = n i=1 p i δ(τ -τ i )
where δ(τ -τ i ) is a Dirac mass at τ i , and ( 9)

n i=1 p i τ i = E,
and

n i=1 p i = 1.
The characteristic equation associated with the density

f n is D n (λ) = λ + a + n i=1 p i e -λτi = 0. Likewise, we denote C n (ω) = n i=1 p i cos(ωτ i ), S n (ω) = n i=1 p i sin(ωτ i ).
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

g(x) : [0, π] → [-1, 1] by g(x) = 1 -cx, 0 ≤ x < θ c , cos(x), θ c ≤ x ≤ π. Convexity implies g(px 1 + (1 -p)x 2 ) ≤ pg(x 1 ) + (1 -p)g(x 2
), for p ∈ [0, 1], and

x 1 , x 2 ∈ [0, π].
In addition, we have g(x) ≤ cos(x).

The following lemmas show how to find the distribution that maximizes S n (ω s ) for n = 2.

Lemma 4.2. Assume a ∈ ]-1, 1[ and E > 0 satisfies E < arccos (-a) ω c , ( 10
)
with ω c = √ 1 -a 2 .
Suppose that there exists ω s ∈ [0, ω c ] and a density f 2 with mean E, such that

(11) C 2 (ω s ) := p 1 cos(ω s τ 1 ) + p 2 cos(ω s τ 2 ) = -a.
Then ω s E < θ c and cos(ω s E) > -a.

Proof. From inequality [START_REF] Beretta | Geometric stability switch criteria in delay differential systems with delay dependent parameters[END_REF], one gets ω c E < arccos(-a) < π, so cos(ω c E) > -a. Moreover, since ω c ≤ 1, the inequality ω s ≤ ω c implies cos(ω s E) ≥ cos(ω c E). Consequently, cos(ω s E) > -a and, using [START_REF] Berezansky | Stability of linear differential equations with a distributed delay[END_REF], we then deduce that cos(ω s E) > C 2 (ω s ). Furthermore, we have

C 2 (ω s ) := p 1 cos(ω s τ 1 ) + p 2 cos(ω s τ 2 ) ≥ p 1 g(ω s τ 1 ) + p 2 g(ω s τ 2 ) ≥ g(ω s E).
The first inequality comes from the definitions of c and g (see Definition 4.1): cos(x) ≥ 1 -cx for x ≥ 0. 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 . Lemma 4.3. Assume a ∈ ]-1, 1
[ and E > 0 satisfies [START_REF] Beretta | Geometric stability switch criteria in delay differential systems with delay dependent parameters[END_REF]. Suppose that there exists ω s ∈ [0, ω c ] and a density f 2 with mean E, such that equality [START_REF] Berezansky | Stability of linear differential equations with a distributed delay[END_REF] is satisfied.

Then there exists a unique density f * with two discrete delays τ * 1 and τ * 2 , mean E, such that τ * 1 = 0 and 0 < ω s τ * 2 ≤ θ c ≤ π, and satisfying

C * 2 (ω s ) = -a. ( 12
)
Proof. Suppose there exists a density f ′ 2 with two discrete delays τ ′ 1 and τ ′ 2 , weights 13) is equivalent to

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 (13) C ′ 2 (ω s ) = -a. By using p ′ 1 = 1 -p ′ 2 and p ′ 2 = E/τ ′ 2 , Eq. (
(14) cos(ω s τ ′ 2 ) = 1 - 1 + a ω s E ω s τ ′ 2 .
From the definitions of the constant c and the function g (Definition 4.1), the equation cos(

x) = 1 -(1 + a)x/ω s E has positive solutions in [0, π] if and only if cos(x) ≥ 1 -cx, that is (15) c ≥ 1 + a ω s E .
To see that inequality [START_REF] Bernard | Modelling Transcriptional Feedback Loops: The Role of Gro/TLE1 in Hes1 Oscillations[END_REF] is indeed satisfied, one can note that, using [START_REF] Berezansky | Stability of linear differential equations with a distributed delay[END_REF],

-a = C 2 (ω s ) = 2 i=1 p i cos(ω s τ i ) ≥ 2 i=1 p i (1 -cω s τ i ) = 1 -cω s E, so -a ≥ 1 -cω s E.
Thus (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

f * = 2 i=1 p * i δ(τ -τ * i ), with p * 2 = E/τ * 2 , p * 1 = 1 -p * 2 ,
E) > 1 -(1 + a)/(ω s E)ω s E. Thus, cos(ω s E)-(1-(1+a)ω s E/(ω s E)) > 0, which implies that ω s E < ω s τ * 2 . Since E = p * 2 τ * 2
, 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 [START_REF] Beretta | Geometric stability switch criteria in delay differential systems with delay dependent parameters[END_REF]. Suppose that there exists ω s ∈ [0, ω c ] and a density f 2 with mean E, such that equality [START_REF] Berezansky | Stability of linear differential equations with a distributed delay[END_REF] is satisfied. Then for any density f 2 with mean E and satisfying Eq. ( 11), we have

S 2 (ω s ) ≤ S * (ω s ),
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:

p 1 = v -T v -u and p 2 = T -u v -u .
By convention, 0 ≤ u < T < v. We consider C 2 (ω s ) and S 2 (ω s ) as functions of u and v; C, S : [0,

T ) × (T, ∞) → [-1, 1] with C(u, v) = v -T v -u cos(u) + T -u v -u cos(v), (16)
S(u, v) = v -T v -u sin(u) + T -u v -u sin(v). ( 17
)
The subscripts 2 have been dropped to ease the reading. Equation ( 17) is to be maximized with the restriction C(u, v) = -a. The equation C(u, v) = -a defines a one-parameter family of solutions. For u fixed, C(u, v) = -a has at most a finite number of solutions v, with v > T . As shown in Lemma 4.3, there is always a solution when u = 0 (equivalent to τ 1 = 0). Fixing u and labeling the solutions in increasing order v i , i = 1, ..., r, r ≥ 1, we have S(u, v 1 ) ≥ S(u, v i ) for all i. Indeed, the equality C(u, v) = -a can be rewritten as

(18) - a + cos(u) T -u v + au + T cos(u) T -u = cos(v).
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

(19) cos(v) = α(u) -β(u)v,
where β(u) > 0 for all u ∈ [0, T ), α(0) = 1 (in case u = 0, Eq. ( 19) reduces to Eq. ( 14)) 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), it follows that all the chords i, i > 1, lie below chord 1, and thus S(u, v 1 ) ≥ S(u, v i ), i ≥ 1.

It is therefore enough to look, for each u, at the smallest solution v 1 of the equation C(u, v) = -a. The solution, which exists for u ∈ [0, T ), can be parameterized by u, with

v 1 = v 1 (u) = min{v|C(u, v) = -a}. At u = 0, the solution A 0 0.5 1 1.5 2 2.5 3 -1.0 -0.8 -0.6 -0.4 -0.2 v 2 v 1 T u v B -1 0 1 -1.0 -0.8 -0.6 -0.4 -0.2 u v 2 v 1 cos sin Figure 1.
Two solutions in v of Eq. ( 19) for fixed u, with density f 2 (τ ) = 0.8δ(τ -0.625) + 0.2δ(τ -3.5) and a = -0.5 (parameter chosen to satisfy the conditions of Lemma 4.2). Then ω c = √ 1 -a 2 = 0.8660, E = 1.2 < arccos(-a)/ω c = 1.2092. Eq. ( 11) was solved for ω s = 0.8308 < ω c to yield T = ω s E = 0.9969 and u = ω s τ 1 = 0.5192). The solution v 2 = ω s τ 2 = 2.9078 corresponds to the density f 2 and v 1 = 1.3056, to the density f * = 0.3925δ(τ -0.625) + 0.6075δ(τ -1.5715). (A) Solutions along the cosine. (B) Solutions parametrized on the circle, illustrating that at the intersection of the secant at -a, the value of S(u, v 1 ) (•) is strictly larger than S(u, v 2 ) (⋄).

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

d du S(u, v 1 (u)) = ∂S ∂u + ∂S ∂v dv 1 du .
If ∂S/∂v < 0, the total derivative is strictly negative if and only if

(20) dv 1 du > - ∂S ∂u / ∂S ∂v .
The partial derivative with respect to v is

∂S ∂v = T -u v -u sin(u) -sin(v) v -u + cos(v) .
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 cos(T )

> -a = α(u) -β(u)T , one gets cos(v) > α(u) -β(u)v for v ∈ [T, π]. It follows that for all v > π, α(u) -β(u)v < α(u) -β(u)π < cos(π) = -1 ≤ cos(v),
and Eq. ( 19) has no root, yielding a contradiction. The sine function is strictly concave on the interval [0, π] and this implies that

(21) sin(u) < sin(v) + (u -v) d dv sin(v),
or equivalently that (sin(u) -sin(v))/(v -u) + cos(v) < 0, for all 0 ≤ u < v ≤ π. This shows that ∂S/∂v < 0. Now,

dv 1 du = v -T T -u cos(v) -cos(u) + (v -u) sin(u) cos(u) -cos(v) -(v -u) sin(v) , - ∂S ∂u / ∂S ∂v = v -T T -u sin(v) -sin(u) -(v -u) cos(u) sin(u) -sin(v) + (v -u) cos(v) .
Inequality ( 20) can be re-expressed as

(v -u) 2 -2 cos(v -u) -(v -u) sin(v -u) > 0.
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 ≤ π:

d dz 2 -2 cos(z) -z sin(z) = sin(z) -z cos(z) > 0.
The last inequality is obtained with inequality [START_REF]Bifurcation and bistability in a model of hematopoietic regulation[END_REF]. 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 [START_REF] Beretta | Geometric stability switch criteria in delay differential systems with delay dependent parameters[END_REF]. Theorem 4.5. Assume a ∈ ]-1, 1[ and E > 0 satisfies inequality [START_REF] Beretta | Geometric stability switch criteria in delay differential systems with delay dependent parameters[END_REF]. Let f n be a discrete density with n ≥ 1 delays and mean E, then the density f n is stable.

Proof. Case n = 1. Single delay distributions (n = 1) are stable by Theorem 2.5.

Case n = 2. Consider a density f 2 with two delays τ 1 < τ 2 . If C 2 (ω s ) > -a for every ω s ∈ [0, ω c ], Corollary 1 states that f 2 is stable. Suppose C 2 (ω s ) = -a for a value ω s ∈ [0, ω c ]. From Lemmas 4.3 and 4.4, there exists a density f * with τ * 1 = 0 and 0 < τ * 2 ≤ π/ω s such that C * (ω s ) = C 2 (ω s ) and S * (ω s ) ≥ S 2 (ω s ). Since S * maximizes the value of S 2 , if we are able to show that any distribution f * with a zero and a positive delay, and C * (ω s ) = -a, satisfies S * (ω s ) < ω s , then from Corollary 1 all distributions with two delays will be stable.

Let the density f * (τ ) = (1 -p)δ(τ ) + pδ(τ -τ * ) with p ∈ (0, 1] and τ * ∈ [E, π/ω s ]. We have C * (ω s ) = 1 -p + p cos(ω s τ * ) = -a. We must show that S * (ω s ) = p sin(ω s ) < ω s . Summing up the squares of the cosine and the sine, we then obtain p 2 = (-a + p -1) 2 + S * 2 (ω s ), so S * (ω s ) = p 2 -(-a + p -1) 2 . Since E satisfies inequality [START_REF] Beretta | Geometric stability switch criteria in delay differential systems with delay dependent parameters[END_REF], then τ * < arccos(-a)/p √ 1 -a 2 . From C * (ω s ) = -a we get ω s = arccos(-(a + 1 -p)p -1 )/τ * . Thus,

p 1 -a 2 arccos -(a + 1 -p)p -1 arccos(-a) < arccos -(a + 1 -p)p -1 τ * = ω s .
Since (a + 1 -p)p -1 ≥ a for p ∈ (0, 1] and a ∈ ]-1, 1[, we have the following inequality arccos(-a)

√ 1 -a 2 ≤ arccos -(a + 1 -p)p -1 1 -(a + 1 -p)p -1 2 , which implies p 1 -(a + 1 -p)p -1 2 ≤ p 1 -a 2 arccos -(a + 1 -p)p -1 arccos(-a) .
Thus,

S * (ω s ) = p 2 -(-a + p -1) 2 ≤ p 1 -a 2 arccos -(a + 1 -p)p -1 arccos(-a) < ω s .
This completes the proof for the case n = 2.

Case n > 2. For densities f with n > 2 delays, the strategy is also to find an upper bound for the value of S(ω s ) via a new distribution f * that keeps C(ω s ) = -a constant. If, for the new distribution, S(ω s ) ≤ S * (ω s ) < ω s holds true, then Corollary 1 can be applied. The construction of f * requires two or three steps. In the first step, all delays τ i > π/ω s are replaced by smaller delays τ ′ i < π/ω s , in order to use the concavity of the sine function on the interval [0, π] as done in the proof of Lemma 4.4, in the following way:

τ ′ i = τ i -2k i π/ω s if sin(ω s τ i ) ≥ 0, 2(k i + 2)π/ω s -τ i if sin(ω s τ i ) < 0,
where k i = max{j|2jπ/ω s ≤ τ i }. This transformation preserves C(ω s ): cos(ω s τ ′ i ) = cos(ω s τ i ), and ensures that S(ω s ) increases: sin(ω s τ ′ i ) = | sin(ω s τ i )|. That way, we obtain an associated delay density f ′ with C ′ (ω s ) = -a, S ′ (ω s ) ≥ S(ω s ), E ′ ≤ E and τ ′ i ≤ π/ω s . In the second step, we reduce the number of strictly positive delays. All pairs of delay τ ′ i < τ ′ j for which the inequality

p i cos(ω s τ ′ i ) + p j cos(ω s τ ′ j ) p i + p j ≤ cos ω s p i τ ′ i + p j τ ′ j p i + p j (22)
holds are iteratively replaced by one positive and one vanishing delay, as done in Lemma 4.3. We note that inequality [START_REF] Cooke | Discrete delay, distributed delay and stability switches[END_REF] 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 density f with one delay τ1 = 0 and two or more delays τk > 0, k = 2, . . . , m, m ≥ 3, such that pi cos(ω s τi ) + pj cos(ω s τj ) pi + pj > cos ω s pi τi + pj τj pi + pj , 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 = p1 and p ′′ 2 = 1 -p1 , C ′′ (ω s ) < C(ω s ) (from inequality [START_REF] Crauste | Stability and hopf bifurcation for a first-order delay differential equation with distributed delay[END_REF]), 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 ′′ . Indeed, this consists in finding (p * 2 , τ * 2 ) such that p * 2 τ * 2 = E ′′ = p ′′ 2 τ ′′ 2 , τ * 2 < τ ′′ 2 , and 1 -p * 2 + p * 2 cos(ω s τ * 2 ) = -a. Hence, this is equivalent to finding τ * 2 ∈ ]E ′′ , τ ′′ 2 [ such that χ(τ * 2 ) := 1 - p ′′ 2 τ ′′ 2 τ * 2 + p ′′ 2 τ ′′ 2 τ * 2 cos(ω s τ * 2 ) = -a. Since χ is continuous, with χ(E ′′ ) = cos(ω s E ′′ ) = cos(ω s Ē) ≥ -a, and χ(τ ′′ 2 ) = C ′′ 2 (ω s ) < -a,
= E ′′ = p ′′ 2 τ ′′ 2 , that p * 2 sin(ω s τ * 2 ) ≥ p ′′ 2 sin(ω s τ ′′ 2 ), or equivalently, S * (ω s ) ≥ S ′′ (ω s ).
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 [START_REF] Beretta | Geometric stability switch criteria in delay differential systems with delay dependent parameters[END_REF], 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 [START_REF] Beretta | Geometric stability switch criteria in delay differential systems with delay dependent parameters[END_REF]. 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 [START_REF] Beretta | Geometric stability switch criteria in delay differential systems with delay dependent parameters[END_REF], so µ n < 0. Then,

D(λ n ) = λ n + a + ∞ 0 e -λnτ dη(τ ) = λ n + a + ∞ 0 e -λnτ d[η(τ ) -η n (τ )] + ∞ 0 e -λnτ dη n (τ ) = ∞ 0 e -λnτ d[η(τ ) -η n (τ )] → 0,
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 η,

D n (λ) = λ + a + ∞ 0 e -λτ dη n (τ ) = λ + a + ∞ 0 e -λτ d[η n (τ ) -η(τ )] + ∞ 0 e -λτ dη(τ ) = ∞ 0 e -λτ d[η n (τ ) -η(τ )] → 0,
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 [START_REF] Basu | G-CSF: function and modes of action[END_REF], 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 [START_REF] Beretta | Geometric stability switch criteria in delay differential systems with delay dependent parameters[END_REF] is still satisfied:

ρE < arccos(-ā) √ 1 -ā2 .
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.

Theorem 5.2. Assume a ∈ ]-1, 1[ and E > 0 satisfies inequality [START_REF] Beretta | Geometric stability switch criteria in delay differential systems with delay dependent parameters[END_REF]. Let η be a delay distribution with mean E, then the distribution η is stable.

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.

Theorem 5.3. The zero solution of Eq. ( 2) is asymptotically stable if a > -b and a ≥ |b|, or if b > |a| and the mean E of η satisfies

(24) E < arccos(-a/b) √ b 2 -a 2 .
The zero solution of Eq. ( 2) may not be asymptotically stable (depending on the particular distribution) if b > |a| and

E ≥ arccos(-a/b) √ b 2 -a 2 .
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 [START_REF] Koury | Erythropoietin retards DNA breakdown and prevents programmed death in erythroid progenitor cells[END_REF]. Platelet production and white blood cell production processes are also controlled by growth factors (thrombopoietin [START_REF] Kaushansky | The molecular mechanisms that control thrombopoiesis[END_REF] and G-CSF [START_REF] Basu | G-CSF: function and modes of action[END_REF], respectively). It is usually thought that mature blood cells act negatively, through growth factor release, on precursors, progenitors and HSCs dynamics [START_REF] Colijn | A mathematical model of hematopoiesis -I. Periodic chronic myelogenous leukemia[END_REF][START_REF]A mathematical model of hematopoiesis -II. Cyclical neutropenia[END_REF].

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 [START_REF] Bélair | Age-structured and two-delay models for erythropoiesis[END_REF][START_REF] Bernard | Oscillations in cyclical neutropenia: New evidence based on mathematical modeling[END_REF][START_REF] Colijn | A mathematical model of hematopoiesis -I. Periodic chronic myelogenous leukemia[END_REF][START_REF]A mathematical model of hematopoiesis -II. Cyclical neutropenia[END_REF][START_REF] Mackey | Unified hypothesis of the origin of aplastic anaemia and periodic hematopoiesis[END_REF][START_REF] Stiehl | Characterization of stem cells using mathematical models of multistage cell lineages[END_REF]. 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 [START_REF] Colijn | A mathematical model of hematopoiesis -I. Periodic chronic myelogenous leukemia[END_REF][START_REF]A mathematical model of hematopoiesis -II. Cyclical neutropenia[END_REF] proposed a compartment model of hematopoiesis, based on previous models of hematopoietic stem cell dynamics [START_REF] Mackey | Unified hypothesis of the origin of aplastic anaemia and periodic hematopoiesis[END_REF], white blood cell dynamics [START_REF] Bernard | Oscillations in cyclical neutropenia: New evidence based on mathematical modeling[END_REF], platelet dynamics [START_REF] Apostu | Understanding Cyclical Thrombocytopenia: A mathematical modeling approach[END_REF] and red blood cell dynamics [START_REF] Bélair | Age-structured and two-delay models for erythropoiesis[END_REF]. 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 [START_REF] Colijn | A mathematical model of hematopoiesis -I. Periodic chronic myelogenous leukemia[END_REF][START_REF]A mathematical model of hematopoiesis -II. Cyclical neutropenia[END_REF] has been further justified and numerically analyzed by Colijn and Mackey [START_REF]Bifurcation and bistability in a model of hematopoietic regulation[END_REF] and Lei and Mackey [START_REF] Lei | Multistability in an age-structured model of hematopoisis: Cyclical neutropenia[END_REF], 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 [START_REF] Stiehl | Characterization of stem cells using mathematical models of multistage cell lineages[END_REF]. 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 z = 3 i=1 p i z i , where p i ≥ 0 and 3 i=1 p i = 1. The HSC number x is governed by the equation [START_REF] Eurich | Distributed delays stabilize ecological feedback systems[END_REF] ẋ = P (x, z) -αx.

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 ( 26)

     ẏ(1) i = β i x -y (1) i , ẏ(j) i = β i y (j-1) i -y (j) i , j = 1, . . . , q i -1, żi = β i y (qi-1) i -z i .
In the i-th hematopoietic lineage, the cell number in generation j-th is denoted by y (j) i , j = 1, . . . , q i -1. Mature cells z i form compartment q i , and immature cells x compartment 0. System ( 25)-( 26) describes a hierarchical structure with parallel negative feedback loops of length q i , with kinetic parameters β i , i = 1, 2, 3. 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 [START_REF] Cooke | Discrete delay, distributed delay and stability switches[END_REF][START_REF] Macdonald | Biological delay systems: linear stability theory[END_REF]. For each lineage i, thanks to the usual chain trick in System [START_REF] Hale | Functional differential equations with infinite delays[END_REF], the repressors z i can be expressed in terms of the history of x convoluted by a Gamma distribution,

z i (t) = ∞ 0 x(t -τ )g(τ, q i , β i )dτ, with g(τ, q i , β i ) = β qi i Γ(q i ) τ qi-1 e -βiτ .
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),

(27) z(t) = ∞ 0 x(t -τ )f p (τ )dτ,
where the density of the distributed delay is a weighted average of Gamma densities,

f p (τ ) = 3 i=1 p i g(τ, q i , β i ).
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 [START_REF] Hale | Phase space for retarded equations with infinite delay[END_REF], one can then write the following equation for x, from ( 25) and ( 27), [START_REF] Hale | Introduction to functional differential equations[END_REF] ẋ

= P x, ∞ 0 x(t -τ )f p (τ )dτ -αx.
The dynamics of System ( 25)-( 26) is entirely contained in [START_REF] Hale | Introduction to functional differential equations[END_REF]. 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)- [START_REF] Hale | Functional differential equations with infinite delays[END_REF]. 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 state x ≥ 0. Eq. ( 28) linearized around a positive steady state x > 0 is

ẋ = -α(1 -r)x - α 2 h k 0 xh-r+1 ∞ 0 x(t -τ )f p (τ )dτ. ( 30
)
For positive cooperativity (1 < r ≤ h), there is a stable steady state x0 = 0. In addition, there are either zero, one or two positive steady states given by the roots of the equation αx h -k 0 xr-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 state x1 satisfies a ≤ -b(x 1 ) and, by Theorem 5.3, is always unstable. The larger steady state x2 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 state x1 of ( 28) is unstable, and the larger positive steady state x2 is linearly asymptotically stable if [START_REF] Hutchinson | Circular causal systems in ecology[END_REF] E p := 2 .

3 i=1 p i q i β i < arccos (r-1)k0 αh(x2) h-r+1 α (αh(x 2 ) h-r+1 /k 0 ) 2 -(r -1)
When x1 = x2 , the positive steady state is unstable. The zero steady state x0 = 0 is always linearly stable.

For negative cooperativity (0 ≤ r < 1), there exists a steady state x0 = 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 + xh )x 1-r = k 0 . The linear equation is given by equation [START_REF] Huang | An analysis of delay-dependent stability for ordinary and partial differential equations with fixed and distributed delays[END_REF], and the instantaneous coefficient is a

= α(1 -r) > 0, the delayed coefficient is b = α 2 hx h-r+1 /k 0 > 0.
For neutral cooperativity (r = 1), there is a steady state x0 = 0, whose stability depends on the existence of a positive steady state. There exists a positive steady state x = ((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 .

Proposition 3 (neutral and negative cooperativity). Assume P is given by [START_REF] Hayes | Roots of the transcendental equation associated with a certain differencedifferential equation[END_REF]. When r = 1 (mixed feedback loop with neutral cooperativity), a unique positive steady state, x = ((k 0 -α)/α) 1/h , of [START_REF] Hale | Introduction to functional differential equations[END_REF] exists if k 0 > α. If it exists, it is linearly asymptotically stable if

E p := 3 i=1 p i q i β i < k 0 π 2αh(k 0 -α) .
The zero steady state x0 = 0 is stable if k 0 < α and unstable otherwise. When 0 ≤ r < 1 (negative cooperativity), a unique positive steady state x of (28) exists. It is linearly asymptotically stable if xh+1-r ≤ k 0 (1 -r)/(αh), or if xh+1-r > k 0 (1 -r)/(αh) and ( 31) holds true. The zero steady state x0 = 0 is always unstable.

A model with neutral cooperativity has been considered before by Mackey and Glass [START_REF] Mackey | Oscillation and chaos in physiological control systems[END_REF] 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. 28) with a distributed delay and with a discrete delay for varying values of r and h. Fixed parameters values are α = 1 and k 0 = 2 (so x = 1 is a steady state of [START_REF] Hale | Introduction to functional differential equations[END_REF], whatever the values of r and h), and q = {2, 20, 60}, β = {2, 10, 20}, p = {0.3, 0.4, 0.3}. (A) The distributed delay (shaded area) 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 state x = 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).

  a = -∂ ∂x F (x, y) x=y=x and b = -∂ ∂y F (x, y) x=y=x .

Definition 2 . 3 .

 23 The function D : C → C given by

Theorem 2 . 5 .

 25 Let f (τ ) = δ(τ -E) a Dirac mass at E. The zero solution of Eq. (2) is asymptotically stable if and only if a > -b and a ≥ |b|, or if b > |a| and

  ω s ) := m k=1 pk cos(ω s τk ) = -a ≤ cos(ω s Ē).The third step is to replace all positive delays τk , k = 2, . . . , m, with the single mean delay

Figure 2 .

 2 Figure 2. Stability chart of distributions of delay in the (a, b) plane, obtained from Theorem 5.3. The asymptotic stability region is composed of regions (1) to (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 Figure1).

Figure 3 .

 3 Figure 3. Comparison of solutions of Eq. (28) with a distributed delay and with a discrete delay for varying values of r and h. Fixed parameters values are α = 1 and k 0 = 2 (so x = 1 is a steady state of[START_REF] Hale | Introduction to functional differential equations[END_REF], whatever the values of r and h), and q = {2, 20, 60}, β = {2, 10, 20}, p = {0.3, 0.4, 0.3}. (A) The distributed delay (shaded area) 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 =

  there is at least one τ * 2 ∈ ]E ′′ , τ ′′

	2 [ satisfying the above conditions, 2 and the function sin(x)/x is decreasing 2 < τ ′′ 2 . Moreover, since τ * 2 /τ * 2 τ ′′ 2 := p ′′ with p * on (0, π), one obtains, using p * 2 τ * 2

Acknowledgments

This work has been supported by ANR grant ProCell ANR-09-JCJC-0100-01.