Stability analysis of an equation with two delays and application to the production of platelets

. We analyze the stability of a diﬀerential equation with two delays originating from a model for a population divided into two subpopulations, immature and mature, and we apply this analysis to a model for platelet production. The dynamics of mature individuals is described by the following nonlinear diﬀerential equation with two delays: x (cid:48) ( t ) = − γx ( t )+ g ( x ( t − τ 1 )) − g ( x ( t − τ 1 − τ 2 )) e − γτ 2 . The method of D-decomposition is used to compute the stability regions for a given equilibrium. The centre manifold theory is used to investigate the steady-state bifurcation and the Hopf bifurcation. Similarly, analysis of the centre manifold associated with a double bifurcation is used to identify a set of parameters such that the solution is a torus in the pseudo-phase space. Finally, the results of the local stability analysis are used to study the impact of an increase of the death rate γ or of a decrease of the survival time τ 2 of platelets on the onset of oscillations. We show that the stability is lost through a small decrease of survival time (from 8.4 to 7 days), or through an important increase of the death rate (from 0.05 to 0.625 days − 1 ).


1.
Introduction. Differential equations with two delays arise when a system includes two "non-instantaneous" processes requiring a finite time to be completed. Take for example a population composed of immature individuals and mature individuals, such that immature individuals become mature after a time τ 1 and mature individuals die after having been mature for a time τ 2 . If we assume that at any time t ≥ 0, the rate of production of immature individuals is a positive function g of the total population x(t) of mature individuals, then the dynamics of x is described by the following differential equation: Adding a random destruction rate γ, Bélair et al. [2] formulated a model for the production of platelets whose dynamics are given by (1.1) The derivation of this equation from the structured PDEs describing the two populations can be found in [10]. The immature cells are megakaryocytes, whose production rate is a decreasing function of the platelet count, and which release between 1000 and 3000 platelets after a maturation time τ 1 . In the context of a disease called "cyclic thrombocytopenia" [11], Bélair et al. [2] found that in the case where the function g has a bell-shape given by g(x) := f 0 θ n x θ n +x n , increasing the death rate γ induces a de-stabilization of the positive equilibrium. Equation (1.1) with the same function g has also been studied by El-Morshedy et al. [10] as authors identified conditions for the global stability of either the trivial equilibrium or the positive equilibrium.
Notice that if x * is an equilibrium of (1.1), then the corresponding characteristic equation is given by λ + γ − g (x * )e −λτ1 + g (x * )e −γτ2 e −λ(τ2+τ1) = 0, (1.2) a particular case of the general form λ + A + Be −λr1 + Ce −λr2 = 0. (1.3) Stability with respect to A, B and C has been studied by Mahaffy et al. [21] in the case where r 1 > r 2 , and by Mahaffy & Busken [7] when r 1 = nr 2 with n ∈ N. Besse also studied this general form [4], although with respect to A and B and in the case where A, B ≥ 0 and −π/r 1 ≤ C ≤ 0. Bélair & Campbell [1] studied the case where A = 0 and B = 1: they identified the stability regions in the (A, r 2 ) plan, and were able to show that increasing r 1 leads to the separation of the stability region into multiple disconnected regions. Although Equation (1.2) constitutes a more complex case than the one treated by Mahaffy & Busken [7], as C depends on r 1 , r 2 and A, this specificity reduces the number of parameters by one, allowing for a more complete analysis of the role each parameter plays in shaping the regions of stability.
The assumption of a regulation of platelet count relying on variable megakaryocyte production was recently studied by Boullu et al. [6]. Using a non constant maturation rate for megakaryoblasts (cells from which megakaryocytes originate), this lead to a system of two delay-differential equations. The authors showed that increasing the death rate of megakaryoblasts could induce periodic solutions. But the use of random-only platelet death prevented from studying the effect of an increased destruction of platelets. The objective of this paper is to study in detail the local stability of the equilibria of (1.1). Then, (1.1) is used as a simplification of the model presented in [6] to study the impact of an increased destruction of platelets on the stability, through the increase of the random death rate γ as well as through the decrease of survival time τ 2 . Indeed, both the increase of γ and the decrease of τ 2 may be induced by the presence of platelet-specific antibodies [12,19,22].
In Section 2, we perform the linear stability analysis of (1.1). Since for B := g (x * )τ 1 = 0 the non-trivial equilibrium is locally asymptotically stable, we study the existence of purely imaginary eigenvalues as |B| increases from zero for different values of τ := τ 2 /τ 1 and A := γτ 1 . Following the D-decomposition approach [13], we obtain implicit expressions for τ , A and B allowing us to numerically plot the curves where λ = iω is a root of (1.2), in both the (τ, B) and (A, B) planes. In order to determine whether the existence of a pair of purely imaginary eigenvalues implies a change in the number of eigenvalues with positive real parts, we study the sign of the derivative of Re(λ) with respect to one of the parameters. This enables us to associate each region to a number of eigenvalues with positive real part. In order to determine the nature of the changes of dynamic associated with a loss of stability, we perform a centre manifold analysis for the steady-state bifurcation and for a single Hopf bifurcation in Section 3. Then, we notice that the crossing of stability curves indicates that complex behaviours, such as tori, are possible. To identify parameters associated with these complex behaviours we perform the centre manifold analysis of a double Hopf bifurcation in Section 4. Finally, the results of Section 2 are applied in Section 5 to (1.1) the equation for platelet count presented above.

2.1.
Curves associated with purely imaginary eigenvalues. In order to decrease the complexity further, we rescale time in units of one of the delays, setting s = t/τ 1 . Thus, if z is a solution of the linearization of (1.1) about one of its equilibrium x * , that is, (2.1) Then we define a function y as y(s) = z(sτ 1 ) and obtain where A = γτ 1 > 0, B = g (x * )τ 1 and τ = τ 2 /τ 1 > 0. Notice that γτ 2 = Aτ , such that the corresponding characteristic equation is written there is only one eigenvalue and it is real and negative, therefore the equilibrium is locally asymptotically stable. Because the number of eigenvalues λ of Equation (2.3) with Re(λ) > 0 can change only by a crossing of λ through the imaginary axis, the boundaries of the stability regions in the plane of parameters correspond to the values such that there exists a ω ∈ R + with λ = iω solution of (2.3). For such a λ, (2.3) becomes In the case ω = 0, (2.4) leads to so that λ = 0 is an eigenvalue if and only if B = A 1−e −Aτ . If λ = iω = 0, we need to separate real and imaginary parts: (2.5) For B = 0, we divide one equality by the other to obtain and when ω, A and τ are known, B is given by .
In order to determine whether the pair of purely imaginary eigenvalues corresponds to an increase or a decrease in the number of eigenvalues with positive real parts, we compute both d Re λ dA (A) and d Re λ dτ (τ ), then evaluate the transversality condition, that is whether the corresponding derivatives are positive or not.
and thus dλ dA Using the general decomposition x+iy u+iv = xu+yv u 2 +v 2 + i yu−xv u 2 +v 2 of any complex number, we obtain the sign of d Re λ dA as the sign of Region of stability

2.2.2.
With respect to τ . We consider again λ = α + iω ∈ C, A, τ ∈ R * + and B ∈ R. Differentiating (2.3) with respect to τ gives , Using the same decomposition as above, we get the sign of d Re λ dτ as the sign of We showed above that on the curve corresponding to ω = 0, that B =   Remark 1: in each graphs, we notice that d Re λ dτ > 0 at a point (τ, B) if and only if the implicit curve B → ω(B) is concave at this point. Verifying this requires to compute the corresponding second derivative, which is out of the scope of this paper and left for a future work. From the graphs again, we assume that there exists a value A 0 such that if A ≥ A 0 , then for all points of the plane (τ, B) there exists no more than one positive value of ω such that (ω, A, τ, B) satisfies (2.5), that is no curve crosses itself or another branch). Region of stability Remark 2: to show that the branch for ω = 0 is always below the other branches corresponding to positive values of ω, we study the consequences of B <
In the next section we perform a centre manifold analysis for the bifurcation for B > 0 (with λ = 0, a steady-state bifurcation) and for the bifurcation for B < 0 (a Hopf bifurcation).

Region of stability with transversality condition for A = 1
Region of stability 3. Centre manifold analysis for simple bifurcations. When an eigenvalue with positive real part appears through a steady-state bifurcation, the change of stability is generically of one of two types. Either one new locally stable equilibrium point appears, in which case the bifurcation is called "transcritical", or two new locally stable equilibrium points appear, in which case the bifurcation is called "pitchfork bifurcation". Similarly, when a pair of eigenvalues with positive real parts appears through a Hopf bifurcation, the change in stability is generically of one of two types. If a stable limit cycle appears when the Hopf bifurcation occurs, it is said to be a "supercritical Hopf bifurcation", and if an unstable limit cycle disappears when the Hopf bifurcation occurs, it is said to be a "subcritical Hopf bifurcation". In both instances, to determine the type of the bifurcation we can perform a centre manifold analysis for the steady-state bifurcation and for simple Hopf bifurcation.
3.1. Theoretical basis for centre manifold analysis (from [8]). We consider the general delay differential equation with two delays expressed as

Region of stability with transversality condition for A = 2
Region of stability where τ 1 , τ 2 ≥ 0, L : R 3 → R a linear operator and f : R 3 → R. The linear restriction of the equation is therefore At a point in parameter space where (3.2) possesses m eigenvalues with zero real parts, all the other eigenvalues having negative real parts, there exists an mdimensional invariant manifold C = P ⊕ Q in the state space. The long term behavior of solutions of the nonlinear equation is well approximated by the flow in this manifold. The flow on this centre manifold is given by where Φ is a basis for P, h ∈ Q and z satisfies the ordinary differential equatioṅ Here, B is the (m × m) matrix of eigenvalues with null real part. To specify b, we introduce the bilinear form associated to (3.2): Then, b is given by We now perform the centre manifold analysis for the steady-state bifurcation.
3.2. Steady-state bifurcation. In order to study the bifurcation corresponding to the branch λ = 0, we use the result of Song & Jiang on a general delay differential equation ( [26], Theorem 1). We use a rescaled version of (1.1): with A, τ > 0, which yields the characteristic equation with B = g (x * ) and x * is an equilibrium of (3.5). We introduce the variable µ defined as µ := B − A 1−e −Aτ , such that for µ = 0, λ = 0 is a root of (3.6). Therefore, we choose Φ = φ(θ) := 1 as a basis of P , and given (3.4), Ψ, Φ = 1 implies that is a basis of P * the dual space of P . The components of the decomposition of (3.5) as We then compute the linear term L 1 of the Taylor expansion of L and the second term F 2 of the Taylor expansion of F on Φx + y where x ∈ R and y ∈ Q : [26], the normal form for the steady-state bifurcation is given byẋ Furthermore, using Theorem 1 of [26] we obtain the following result: 1−e −Aτ , then λ = 0 is an eigenvalue of (3.6) and there are no other eigenvalues with zero real part. Also, the constant solution x * of (3.5) is locally asymptotically stable for B < B 0 and unstable for B > B 0 .
Furthermore, if g (x * ) = 0, then the constant solution x * with B = B 0 is stable and equation (3.5) undergoes a transcritical bifurcation at the critical value B = B 0 .
Remark: if g (x * ) = 0, the possibility of a pitchfork bifurcation can be evaluated using a normal form of higher order than the one in (3.7). However, this is out of the scope of this paper.
We now perform the centre manifold analysis for a single Hopf bifurcation.
3.3. Single Hopf. In the following, we employ the Taylor expansion of (1.1) with scaling about an equilibrium x * , given by In order to obtain the type of the Hopf bifurcation, we use the method of Campbell [8] for compute the centre manifold using the symbolic algebra package Maple. We give an overview of the different steps of computation: • we start by computing the vector b using (3.4) and the fact that for the single Hopf bifurcation we have Φ = (φ 1 , φ 2 ) = (sin(ωθ), cos(ωθ)) ; • we then introduce the function and we compute the functions h ij by solving the following partial differential equation (see [8], eq. 8.31): The arbitrary constants are determined using a second partial differential equation given in [8] (eq. 8.32) ; • the equation (3.3) becomes and the type of the Hopf bifurcation is determined (see [14], eq. 3.4.11) by the sign of (3.9) The Maple commands are available online 1 . The final expression of a involves derivatives of g of order higher than two, such that it is not possible to compute its value without choosing a function g beforehand (see next section). Because this implies a constraint on B, the value a is computed on single points rather than on the whole stability boundary as it can be seen in [1].
In the next section, we present a detailed analysis of the bifurcations occurring at parameter values where two pairs of purely imaginary eigenvalues exist. 4. Centre manifold analysis for double Hopf bifurcation. In the double Hopf case, the type of the bifurcation is indicated by four coefficients a 11 , a 12 , a 21 , a 22 which are the equivalent of a in the simple Hopf case. However, there is no explicit expression of these coefficients (an equivalent to (3.9)) in the literature. Therefore they need to be computed using the symbolic algebra package Maple. Computing a 11 , a 12 , a 21 , a 22 using Maple. The computation of a 11 , a 12 , a 21 , a 22 from a system of the forṁ
We give an overview of the different steps of computation performed in Maple : • We introduce the Taylor expansion to the second order of Φ, Ψ, h, g, that we substitute in PDE1 and PDE2. Equating coefficients of w i , w j , w i w j and w i w j we obtain expressions of the partial derivatives of Φ, Ψ as functions of partial derivatives of h, g. • We introduce the Taylor expansion to the third order of h, g that we substitute in PDE1 and PDE2. Equating coefficients of w 2 1 w 1 , w 1 w 2 w 2 , w 1 w 1 w 2 and w 2 2 w 2 , we obtain expressions of the coefficients c 11 , c 12 , c 21 , c 22 as functions of the partial derivatives of h, g. The coefficients of interest a 11 , a 12 , a 21 , a 22 are given as a ij = Re(c ij ).
• To obtain expressions of the partial derivatives of h, g as functions of the partial derivatives of f i , i = 1, 2, 3, 4, we use the Taylor expansion to the third order of f i , i = 1, 2, 3, 4 and the fact that The partial derivatives of f i , i = 1, 2, 3, 4 are computed as functions of the F l ijk , i, j, k, l = 1, 2, 3, 4 of (4.1). From there, the computation of the F l ijk , i, j, k, l = 1, 2, 3, 4 of (4.1) is performed following the adaptation of the method of Campbell [8] And the elements needed to write (3.3) are In order to obtain the type of the Hopf bifurcation, we use the method of Campbell [8] for calculating centre manifold using the symbolic algebra package Maple. The algorithm using Maple and its results are available online 4 . For sake of clarity, we decide to only give an overview of the different steps of computation to them: • we start by computing the vector b using (3.4); • we then introduce the function and we compute the functions h ii by solving the following partial differential equation (see [8], equation 8.31): The arbitrary constants are determined using a second partial differential equation (see [8], equation 8.32); • we then obtain the equation (3.3) under the form (4.1).
Combining with the expressions of a 11 , a 12 , a 21 , a 22 as functions of the coefficients of (4.1), we have an expression of a 11 , a 12 , a 21 , a 22 as a function of A, B, C, D, ω and τ . The Maple commands are available online 5 . From the values of a 11 , a 12 , a 21 , a 22 to the descriptions of the flows near the double Hopf. We introduce the continuous functions µ 1 , µ 2 , Ω 1 , Ω 2 on R 2 + describing the branches of eigenvalues associated with the purely imaginary eigenvalues ω 1 and ω 2 , i.e. such that for a point (τ, B),
Because the expressions of a ij involve C and D, we need to fix g such that the point (τ, B) corresponds to the crossing of the lower stability boundary seen on Figure 2.2 : we choose g(x) := f 0 σ n σ n +x n + β 0 with f 0 = 5 × 10 10 , β 0 = 0.01, σ = 8 × 10 9 and n = 10.166. In this case, and assuming that the curve decreases with τ corresponds to ω 1 , we have a 11 > 0, a 12 > 0, a 21 < 0, a 22 < 0, θ := a 12 a 22 < 0, δ := a 21 a 11 < 0, θδ − 1 > 0, (4.9) which corresponds to subcase VI of the "complex" case in Kuznetsov numbering scheme (see [17] p. 363). The Maple commands are available online 6 . Figure 4.1 represents the different phase portraits for the system (4.8) associated with different sectors of the (µ 1 , µ 2 ) plane. We recall that given the definitions of  For Figure 4.2, we have n = 11, τ 2 = 4.75×τ 1 , and the corresponding eigenvalues have real parts µ 1 = −0.0651, µ 2 = 0.0014, such that µ 2 < µ 1 /θ = 0.1023. This implies that it is under the line T 1 defined by Kuznetsov [17], such that this point corresponds to the lowest wedge of the µ 1 < 0, µ 2 > 0 quadrant of Figure 4.1. And as predicted by the unfolding, the solution converges to a stable limit cycle.
For Figure 4.3, we have n = 11, τ 2 = 4.24×τ 1 , and the corresponding eigenvalues have real parts µ 1 = −0.0078, µ 2 = 0.0188, such that µ 2 > µ 1 /θ = 0.0123 and 0334. This implies that it is above the line T 1 and under the line C defined by Kuznetsov [17], such that this point corresponds to the lowest interior wedge of the µ 1 < 0, µ 2 > 0 quadrant of Figure 4.1. As predicted by the unfolding, the solution converges to a torus. For Figure 4.4, we have n = 11, τ 2 = ×τ 1 , and this point corresponds to the µ 1 > 0, µ 2 > 0 quadrant of Figure 4.1. The unfolding cannot be used to predict the behaviour of the solution, and we see that a stable cycle appears.
x(t) In the next section, we use the computation presented in Section 3.3 to determine the criticality of Hopf bifurcation for a model of platelet production.

5.
Application to the production of platelets. As mentioned in the introduction, we now use the results obtained above on a model describing the production of platelets, in order to study the role of the two mechanisms for destruction (random and deterministic). This model may be seen as a simplification of the one given in

Solution in the pseudo-phase plane (y(t-τ 1 -τ 2 ),y(t))
x(t)  Numerical simulation of (1.1) for τ 2 = 4.24×τ 1 in the pseudo-phase space (x(t), x(t − τ 1 ), x(t − τ 1 − τ 2 )), corresponding to the lowest interior wedge of the µ 1 < 0, µ 2 > 0 quadrant of Figure  4.1. Once the transient dynamic is lost, a stable torus appears. [6]. Indeed, we represent the regulation by the platelet number not as a continuous action on the maturity speed of megakaryocyte progenitors, but more directly by a regulation of the rate of production of new megakaryocytes. In order to reproduce the fact that this production rate is bounded below and above, we express it as a function of the platelet count G(P ) = m σ n σ n +P n + m. We consider that τ 1 is the time between the birth of a megakaryocyte and the time its cytoplasm sheds into platelets. Furthermore we assume that platelets die randomly with a rate γ and are removed from the system by macrophages and hepatocytes after a survival time τ 2 [12,19]. Then the model is described by the following equation: with g(P ) = A 0 G(P ) = f 0 σ n σ n +P n +β 0 , A 0 the number of platelet per megakaryocyte.

Parameter estimates.
Based on observations of a delay of 5 days between stimulation by thrombopoietin and a rise in platelet count [16], we choose τ 1 = 5; and we set τ 2 = 8.4 using observed mean platelet survival time [27]. Furthermore, the following biological observations allow us to obtain f 0 , β 0 and σ: • at steady state, the amount of platelets is x * = 20 × 10 9 /kg [23], such that −γx * + f 0 σ n σ n + x * n + β 0 (1 − e −γτ2 ) = 0; x(t) • when the feedback is deactivated, which for the regulation system is equivalent to a virtually infinite amount of platelets, the amount of platelets is divided by 10 [9], such that −γx * /10 + f 0 β 0 (1 − e −γτ2 ) = 0; • the maximum increase reached with an artificial simulation of the hormone signal controling platelet production, equivalent to a population of 0, leads to an 10-fold increase [15], such that Therefore, only the parameters γ and n are not determined by the model. Until appropriate data is available, we decide to use Langlois et al. fitted value of γ = 0.05 [18]. Finally, we choose n = 1.7 such that the unique equilibrium is stable when parameters correspond to healthy patients (that is, τ 2 = 8.4 and γ = 0.05).

Stability analysis.
As g is decreasing, we are interested in the part of the (B, τ ) plane corresponding to B ≤ 0. We study the effect of increasing platelet death rate, that is, increasing γ; and of decreasing the platelet survival time, that is, decreasing τ 2 . Plots of the evolution of B = f 0 g (x * )τ 1 as either γ increases or as τ 2 decreases are given in Figure 5.1. We see that when τ 2 decreases of one day (to τ 2 = 7.2), then the system loses its stability. Furthermore, if γ is 10-fold then the system also loses its stability.  We see that when τ 2 decreases of one day (to τ 2 = 7.2), then the system loses its stability. Furthermore, if γ is multiplied more than 12 times (to γ = 0.625) then the system also loses its stability.
zone of stability, a Hopf bifurcation occurs. Furthermore, the criticality of the Hopf bifurcation is computed in both cases using the expression (3.9) of a computed in Maple (see online 7 ). In both cases, a is negative such that the bifurcation is supercritical and a locally stable periodic solution appears, as seen on the solutions obtained numerically.
Other characteristics of the dynamics change differently depending on the bifurcation parameter: the main consequence of decreasing τ 2 is an increase in amplitude, but increasing γ decreases the overall value of platelet count. However, from a clinical point of view a decrease in platelet survival time τ 2 seems to induce a more threatening chronic thrombocytopenia than a change in γ, as well as inducing a chronic thrombocytosis.

Conclusion.
In an attempt to model the quantity of platelets in the blood as affected by both a random and an age-related destruction, we analyze a delay differential equation with two delays presented by Bélair et al. in 1987 [2]. In order to complete the already existing body of work dedicated to two-delays differential equations, we studied the stability of the equilibrium points for a scaled linearized version of the equation described by three parameters, A, B and τ . As the purely imaginary eigenvalues can not be identified analytically, we used a numerical method, D-decomposition, to obtain curves in the plans (τ, B) and (A, B) corresponding to the existence of a pair of purely imaginary eigenvalues. Associated with a study of the transversality condition in A then in τ , we obtain the number of eigenvalues with positive real parts in the different stability regions. This analysis is then completed with a centre manifold analysis for the steady-state bifurcation, as well as for the single Hopf bifurcation, allowing us to test whether leaving the region of stability through a Hopf bifurcation is associated with the onset of periodic solutions. Furthermore, the existence of self-intersection in the stability boundary associated with Hopf bifurcation revealed the possibility of a double Hopf bifurcation, a phenomenon known to generate a variety of complex dynamical behaviors. Therefore we performed a centre manifold analysis near the point where a double bifurcation occurs, revealing the possibility of torus-like dynamics near the region where two eigenvalues with positive real parts exist. Finally, coming back to the model for platelet production, we use the results on the single Hopf bifurcation to explore the impact of an increase in death rate or a decrease in survival time on the onset of periodic dynamics, as we expected from the platelet-specific antibodies observed in patients with cyclic thrombocytopenia. We choose a megakaryocyte production rate which decreases when the platelet count increases, while always staying strictly positive and finite. The parameters of this feedback function are identified using the normal mean platelet count as an equilibrium. Finally, we show that although stability is gained by going from left to right in the plane (B, γ), taking in account the role of γ in the computation of B = g (x * )τ 1 enables to lose stability when increasing the death rate γ of platelets. The same counter-intuitive conclusion is obtained for a decrease of the survival time of platelets τ 2 . The extent of the change in τ 2 necessary to induce oscillations seems more reasonable than that of γ. In the cases of cyclic thrombocytopenia of the auto-immune type, platelet-specific antibodies are observed in patients blood. Therefore, our work indicates that auto-immune cyclic thrombocytopenia is more likely to be due to a over-reactive system of old platelet removal, rather than to an increase random destruction of platelets of any age. However, this relies on the likelihood of a small decrease of τ 2 being more important than the likelihood of a important increase of γ, which needs to be confirmed with clinicians.
Notice that although regularities appear in the graphs of the stability curve, such as an ordering in the curves corresponding to different intervals [π + 2nπ, 2(n + 1)π], n ∈ N or a equivalence between the transversality condition and the concavity of the function B → B(ω), we were not able to obtain an analytical proof for such results. Besides, we encountered difficulties in simulating the solutions of (1.1) corresponding the sections 1, 2 and 6 of the unfolding (Figure 4.1). These questions are still open and will be the object of a future work.
Most works on the stability of two-delays differential equations focus on a limited version of (1.3). For example, Bélair & Campbell [1] studied the case where A = 0 and B = 1, and focused on the features of the stable regions in the plane (A, r 2 ) as r 1 increases. In particular, authors obtained intersecting stability curves associated with a double Hopf bifurcation. Using centre manifold analysis, they identified two sets of parameters such that near the double Hopf bifurcation, the unfolding corresponds to type Ib of Guckenheimer & Holmes [14] for one, and VIa for the other (which corresponds to the subcase VI of the "complex" case in Kuznetsov numbering scheme (see [17] p. 363)). The authors were also able to show that when r 1 increases, the stability region becomes composed of disconnected stable regions, a feature previously thought to be impossible. Interestingly, while the choice of fixing A to 0 simplifies the computation of the stability regions, the identification of bifurcation type still relies on numerical methods. Using a normalization such that r 1 = 1, Mahaffy et al. [21] focused on the stability regions in the 3D-space (A, B, C). Indeed, once A and r 2 are fixed, the computation of the values of B and C corresponding to the existence of a pair of purely imaginary eigenvalues is straightforward. In particular, the authors identified regions in the plane (R, A) corresponding to different configurations of the stability regions in the plane (B, C), and noticed that there exists two disconnected regions for A < 0. Mahaffy & Busken [7] returned later to the study of the stability regions in the plane (B, C) while this time focusing on the differences observed between the case r 2 ∈ Q and r 2 / ∈ Q. Finally, Besse [4] analyzed the stability of (1.3) by presenting the stability regions for C < 0, r 1 and r 2 fixed. The author introduced a change of variable x = A + B and y = −A + B, such that the stability for x < −C or for y < C is straightforward. The stability for y ≥ C, x ≥ −C, however, is non-trivial, as for high r 2 it involves intertwined loops. The author were nevertheless able to obtain a region of the plane (x, y) in which the system is stable for any r 2 ≥ 0, at the condition C ≥ −π/r 1 . In our case, a change of stability occurs once along an increase of a parameter of the model. However, it is known that such a change in stability is sometimes reversed by increasing the parameter further. This is often the case when the coefficients of the characteristic equation involve the delay, as in [24,25]. In our case, the coefficients do not involve the delays, and for the application that we study, we do not have stability switches. However, it is not impossible that given different parameters or functions our model would present stability switches as the use of two delays implies an increase in complexity such that delay-dependent coefficients are not needed for stability switches. For example, disconnected stability regions have been observed for such models by Bélair and Campbell [1], which is a feature conducive to stability switches.
Although most of the aforementioned papers rely on the method of Ddecomposition, they do not use it to study the effect of parameters changes as we did in Section 5. Therefore we mention two papers dedicated to models of erythropoiesis (the production of red blood cells) applying it. In the first one, Bélair et al. [3] studies the stability of the equilibrium of a system of two differential equations with two delays, and the analysis relies on a characteristic equation given by (λ + γ)(λ + k) = −A(e −λT1 − e −γT1 e −λ(T1+T2) ), where A depends on the steady state. It is a more complicated version of (1.3) with an additional λ 2 . Authors identify numerically the stability region in the plane (γ, A), and similarly to the conclusion presented in Section 5, stability is gained when one moves from left to right in the plane (A, γ) but the role of γ in the computation of A implies that increasing γ can lead to instability. A similar phenomenon was obtained later by Mahaffy et al. [20] in a second paper on a model of erythropoiesis with a state-dependent delay.
Finally, we notice that equation (1.1) can also be obtained by adding a survival time for platelets to a model of megakaryopoiesis whose stability was recently analyzed [5]. In the case of this model, we have τ 1 > τ 2 , which according to preliminary numerical explorations (not shown) seems to imply that stability is not affected of the value of γ.