Analysis of unstable behavior in a mathematical model for erythropoiesis

We consider an age-structured model that describes the regulation of erythropoiesis through the negative feedback loop between erythropoietin and hemoglobin. This model is reduced to a system of two ordinary differential equations with two constant delays for which we show existence of a unique steady state. We determine all instances at which this steady state loses stability via a Hopf bifurcation through a theoretical bifurcation analysis establishing analytical expressions for the scenarios in which they arise. We show examples of supercritical Hopf bifurcations for parameter values estimated according to physiological values for humans found in the literature and present numerical simulations in agreement with the theoretical analysis. We provide a strategy for parameter estimation to match empirical measurements and predict dynamics in experimental settings, and compare existing data on hemoglobin oscillation in rabbits with predictions of our model.


Introduction
Dynamical diseases are characterized as those that occur through changes in the qualitative dynamics of physiological processes. This definition leads naturally to the description of dynamical diseases as a nonlinear system undergoing one or more bifurcations (Mackey and Glass 1977). Indeed, mathematical models have proven to be an effective way to describe several dynamic pathologies (Reimann 1951;Mackey andGlass 1977, 1989;Glass and Mackey 1988;Milton and Mackey 1989;Bélair and Campbell 1994). In this paper, we focus on characterizing irregular dynamics of one specific physiological control system, the regulation of erythropoiesis via the feedback loop between erythropoietin and hemoglobin.
Erythropoiesis, the production of erythrocytes (red blood cells), involves interaction between two organ systems: the bone marrow, where the precursors to mature red blood cells are located, and the kidney, which produces the hormone erythropoietin. When the partial pressure of oxygen in the blood drops, suggesting hypoxia, specialized epithelial cells in the kidney release erythropoietin into the bloodstream. The erythropoietin then travels to the bone marrow where it both accelerates erythropoiesis and recruits additional stem cells into the committed pathway to become mature red blood cells.
Hemoglobin is an iron-containing metalloenzyme essential to erythropoiesis. Because a high percentage of the dry volume of a red blood cell consists of hemoglobin, low levels of this protein result in suboptimal red blood cell production regardless of erythropoietin release. On the other hand, high availability of hemoglobin allows for maximal levels of erythropoiesis. The function of hemoglobin is to act as a transport protein for oxygen in the bloodstream. For this reason, the increase in hemoglobin content (or red blood cell number) caused by elevated erythropoietin raises the oxygen carrying capacity of the blood and corrects hypoxia, thus removing the original pressure for the release of erythropoietin. This is an example of a negative feedback system, in which an excess of product suppresses its production. Many biological control systems involve negative feedback, and the disruption of such a loop can result in severe pathology (Grodins et al. 1954;Foley and Mackey 2009).
A mathematical model of hematopoiesis was introduced by Mackey (1978), who proposed a model of two nonlinear differential equations with one constant delay representing the average cell cycle duration. Mackey's model is called an age-structured model, where populations are classified and tracked according to their age. Modifications of this model have been analyzed to be applied to a variety of hematological diseases by different authors (Bernard et al. 2003;Pujo-Menjouet and Mackey 2004;Adimy et al. 2010a,b). Application of similar models to erythropoiesis can be found in Bélair et al. (1995), Mahaffy et al. (1998) and Adimy et al. (2010b). For an extended review of past work, we refer to Foley and Mackey (2009).
Many studies use bifurcation analysis to provide an understanding of the underlying system dynamics leading to pathology (Mackey 1979;Guevara and Glass 1982). Specifically, periodic diseases can be simulated by systems which undergo a Hopf bifurcation, a local bifurcation in which a limit cycle arises from an equilibrium point, resulting in oscillatory behavior (Kuznetsov 2004). Such a bifurcation in a model of erythropoiesis would result in oscillatory behavior in both erythropoietin and hemoglobin levels (Bélair and Campbell 1994;Bélair et al. 1995;Adimy et al. 2010a). In particular, this behavior has been observed in patients with renal failure treated with erythropoiesisstimulating agents, or ESAs (Fishbane and Berns 2007).
In this paper we consider a modification of the age structured model presented by Bélair et al. (1995) to describe red blood cell production. Since the total population of red blood cells actually consists of two distinct subpopulations, committed precursor cells and mature erythrocytes, the model analyzed in this paper implements two discrete delays corresponding to the lifespan of each of these subpopulations.
We perform a complete theoretical bifurcation analysis to identify all instances for which a supercritical Hopf bifurcation can occur in the system. Specifically, since the hemoglobin-driven production of erythropoietin is modeled with a Hill function, we show that the dynamical model is stable if the Hill exponent is one and determine all the values of this parameter at which Hopf bifurcations arise. The derivation of analytical expressions which describe the conditions for which Hopf bifurcations occur provide a good understanding of how such models are relevant in reproducing oscillatory dynamics associated with dynamical diseases. Furthermore, we provide a procedure to estimate model parameters that can be used to fit our model to experimental results, and demonstrate its ability to predict pathological oscillatory behavior in an experiment involving induced hemolytic anemia in rabbits.
The paper is organized as follows. We introduce the mathematical model for regulation of erythropoiesis in Sect. 2. In Sect. 3, we perform a linear stability analysis after reducing the model to a system of ordinary differential equations with two discrete delays. We determine the unique equilibrium point of the system and derive the characteristic equation for the linearized system around this point. In Sect. 4 we perform a bifurcation analysis providing a complete parameter study and determining the conditions for Hopf bifurcations. In this analysis, we demarcate three subregions of the parameter space under different conditions on the parameters, in which we demonstrate the existence of Hopf bifurcations. In this manner, we identify all possible Hopf bifurcations and establish analytical expressions for the scenarios in which they arise. We present simulations which provide examples of oscillatory dynamics resulting from Hopf bifurcations for physiologically relevant parameters in Sect. 5. In Sect. 6, we provide a simple strategy to estimate model parameters to fit experimental data, and compare the predictions of our model to existing data of dynamical hematopoietic disease in rabbits. A concluding discussion is given in Sect. 7.

Model for regulation of erythropoiesis
In the process of erythropoiesis, new red blood cells are created from precursor cells in the bone marrow at a rate proportional to the amount of the hormone erythropoietin present in the body. These cells age over a period of months through abrasion in the capillaries, eventually losing all ability to transport oxygen, at which point they are destroyed by phage cells. Mature red blood cells carry hemoglobin, the concentration of which is fairly constant among the population of red blood cells and therefore a good measure of red blood cell levels. Erythropoietin and hemoglobin are involved in a negative feedback loop (Graber and Krantz 1978). The mathematical model we consider for modeling these interactions is based on the age-structured system of equations introduced in Bélair et al. (1995).
Let t represent time and ν represent red blood cell age. Let m(t, ν) be the density of red blood cells in the blood with respect to time t and red blood cell age ν. We assume that the cells age at a constant rate α, since the age of a red blood cell depends mainly on the number of times it passes through the capillary system. The level of hemoglobin M(t) is determined as a proportion of the total density of red blood cells and influences the level of the hormone erythropoietin A(t). Consistent with the idea of negative feedback, the concentration of erythropoietin will decrease with increasing numbers of mature erythrocytes.
Because we assume that all red blood cells age at a constant rate and die only when they reach the end of a prescribed lifespan, we can consider the population of mature cells as a conserved quantity. We therefore describe the evolution of the red blood cell population by the transport equation where ν max is the maximum age of a mature cell. We impose the following time-dependent Dirichlet boundary condition where T p is the time required for a stem cell committed to erythropoiesis to reach maturity, and E is a monotone increasing function that represents the rate of erythropoietin-driven production of red blood cells. We consider E to be a Hill function, as is common to enzyme kinetic models (Cornish-Bowden 2004). We express it as where E max and E 50 are the maximum concentration and pre-image of 50 % saturation, respectively. This indicates that as the erythropoietin concentration increases, the amount of incoming red blood cells at time t increases as well, until it asymptotically reaches a saturation point. Note that there is a delay of time T p in Eq. (2) due to the time required for committed stem cells to mature fully. This delay signifies that the input at every time t is dependent on A(t − T p ). We consider A (t) to be known for t ∈ [−T p , 0]. Since hemoglobin concentration is directly related to the total population of red blood cells, we have where R is a scaling constant representing the mean corpuscular hemoglobin-i.e. the amount of hemoglobin in a typical red blood cell. We consider the negative feedback loop between the levels of hemoglobin M(t) and erythropoietin A(t). The time evolution of A(t) is given by the differential equation where k > 0 represents the elimination rate of erythropoietin, and we consider the rate of hemoglobin-driven production of erythropoietin F to be given by a decreasing Hill function where F max and F 50 are defined analogously to E max and E 50 in Eq. (3). The parameter r indicates the speed of response of the rate of hemoglobin-driven production to the level of hemoglobin. Consistent with the concept of negative feedback, an increase in hemoglobin concentration leads to a fall in the rate of erythropoietin production, as seen by the decrease of F(M(t)), and consequently of the left hand side of Eq. (5), as M(t) rises. We review the parameters we have introduced and their associated units in Table 1. The existence of Hopf bifurcations in similar models to the one we have presented have been suggested in different approaches (Mackey 1979;Bélair and Campbell 1994;Bélair et al. 1995). In the following sections we determine the conditions for existence of Hopf bifurcations in the dynamics of the erythropoietin-hemoglobin regulation loop for the complete parameter space. In order to perform a bifurcation analysis we first present a linear stability analysis to determine the conditions under which the unique equilibrium point of the system can lose stability.

Linear stability analysis
The system introduced above consists of three equations describing red blood cell population (Eq. 1), hemoglobin levels (Eq. 4) and erythropoietin levels (Eq. 5). In this section we express the system as a pair of ordinary differential equations with two delays. Such a reduction is natural as it allows us to concentrate on the feedback loop between rates of erythropoietin and hemoglobin production. We then show the existence of a unique, physiologically relevant equilibrium point, and analyze the stability of the system linearized around this point.
3.1 Reduction to a system of ordinary differential equations with delay We first use the method of characteristics to solve the transport equation (1). The directional derivative of m along the vector (1, α) is zero, implying that m is constant in the direction of this vector, and thus Combining Eq. (7) with the boundary condition (2), we can express the solution of the transport equation for large time t and all age ν as We now express the levels of hemoglobin M and erythropoietin A in terms of a system of equations that does not depend explicitly on the partial differential equation (1). Integrating Eq. (1) over age ν and then using Eq. (8) we obtain a differential equation with two delays describing the time evolution of the level of hemoglobin: is the lifespan of erythrocytes. Thus, together with Eq. (5), we describe the time evolution of the levels of hemoglobin and erythropoietin by the following non-linear system of ordinary differential equations:

Existence and uniqueness of steady state
We show existence of an equilibrium point of system (10) and verify its uniqueness. Using the solution of the transport equation as expressed in Eq. (8) and applying a change of variables τ = ν α , we write Eq. (4) as Let us define the values of M(t) and A(t) at equilibrium as M ∞ and A ∞ , respectively. To find an expression for the equilibrium values in terms of the reduced system, we let M (t) = M ∞ and A (t) = A ∞ , also implying d A(t) dt = 0. From Eqs. (11) and (10b) we have Solving Eq. (12b) for A ∞ and plugging the result back into Eq. (12a) we find Every solution to this equation with respect to M ∞ gives an equilibrium point of the system. We see that the right hand side of Eq. (13) takes on a positive value when M ∞ = 0, is negative for large values of M ∞ , and is a continuous function of M ∞ . By the intermediate value theorem, Eq. (13) must have at least one positive solution. Furthermore, the right hand side of Eq. (13) is a monotone decreasing function, and so there exists a unique positive M ∞ such that Eq. (13) holds. Since A ∞ is a function of M ∞ , we conclude that there is a unique equilibrium point (M ∞ , A ∞ ) with both M ∞ and A ∞ positive.

Conditions for stability of the steady state
We study the behavior of the system (10) around the equilibrium point (M ∞ , A ∞ ) through the roots of the characteristic equation, which we derive from the linearized system of equations around the equilibrium point. We linearize the system by calculating the Taylor expansions in the neighborhood of the equilibrium point for the functions E and F As is standard in local stability analysis for nonlinear systems, linearization of system (10) around the equilibrium point is obtained by considering the first order terms of the previous Taylor expansions. Thus,system (10) becomes where we used Eq. (12b) to get Eq. (14b). This system locally describes the time evolution of (M(t), A(t)) around the equilibrium point (M ∞ , A ∞ ). We introduce perturbations from the stationary values Plugging these into system (14) we obtain We are interested in finding the characteristic values of this linearized system. We search for solutions of the formM where λ is the characteristic value of the system. Plugging this form ofM andÃ into system (15), we arrive at the following homogeneous system of linear equations in To find a non-trivial solution the determinant of the system has to be zero. This condition provides the characteristic equation associated with system (14), We note that, as a consequence of the two delays in the model, the characteristic equation is transcendental. The steady state of the linearized system is asymptotically stable if all solutions to the characteristic equation have negative real parts, and is unstable if there exists a solution with positive real part (Hale and Verduyn Lunel 1993). We are interested particularly in the loss of stability through Hopf bifurcations.
Hopf bifurcations occur in the system when a pair of complex eigenvalues crosses the imaginary axis. We therefore search for purely imaginary eigenvalues, i.e. those of the form λ = iω, ω ∈ R. Enforcing this form into Eq. (16) and using trigonometric identities we obtain We present three conditions to ensure that the complex numbers on either side of Eq. (17) are equal. Two conditions are to satisfy equal arguments of the complex numbers: first, we need that the ratios of the imaginary part over the real part of both sides are equal, and second, the real parts on both sides of (17) need to have the same sign. Finally, the third condition is the requirement that the absolute values of both sides need to be equal. The three conditions are represented respectively by the following relations: We name Eqs. (18a), (18b) and (18c) the argument equation, the sign condition and the modulus equation, respectively.
In the following section, we determine the regions of the parameter space for which all conditions in (18) are met, implying that there may exist a Hopf bifurcation in that region.

Bifurcation analysis
In this section, we study the conditions under which the model presented in Sect. 2 can reproduce a dynamical disease. We determine the instances where the equilibrium point loses stability through a Hopf bifurcation. A Hopf bifurcation occurs when a pair of complex conjugate eigenvalues of the system crosses the imaginary axis through the varying of one or more parameters (Strogatz 2000;Kuznetsov 2004). In order to find Hopf bifurcations, we first find all purely imaginary solutions of the characteristic equation, i.e. when Eqs. (18a) and (18c) have a solution simultaneously and condition (18b) is satisfied. We then verify that we indeed have Hopf bifurcations at these points.
If λ = iω is a root of the characteristic equation, ω must satisfy Eqs. (18a), (18b), and (18c). The possible values of ω that satisfy the argument equation (18a) are characterized in the following proposition.
such that ω is a solution of the argument equation (b) There are infinitely many solutions ω > 0 to the argument equation such that the "sign condition" (18b) is satisfied.
Determining when ω satisfies the modulus equation is more involved. Let us rewrite modulus equation (18c) as and define the function g as the right hand side of this equation The left hand side of Eq. (21) is found in the characteristic equation (16) as the coefficient of a combination of exponential functions, the absolute value of which is bounded by 2 when λ is purely imaginary. We note that for imaginary λ, the term |λ 2 + λk| in the characteristic equation can be arbitrarily large. As we wish to consider arbitrary values of |λ|, the named coefficient must be able to take on arbitrarily large values to fulfill the characteristic equation. Thus, we focus on the exponent r , as it is directly responsible for the speed of increase of this term. This exponent will act as the key parameter in the following analysis.
We express the left hand side of Eq. (21) as the function b depending on r .
As we have a large number of parameters in our model, we facilitate our analysis by introducing two strictly positive dimensionless parameters, P and Q: where D and G are defined as the upper bounds of M ∞ and A ∞ , respectively: We will use the dimensionless parameters P and Q to perform the bifurcation analysis. These allow a simple and precise description of the imaginary roots of the characteristic equation.
The next proposition determines the dependence of the equilibrium point on r .
Proof Using Eq. (13) and the definition of P and Q, we obtain We eliminate the constant F 50 by introducing a dimensionless positive function ψ(r ): We substitute (31) into (30) and obtain implying that S(ψ) = 0. As S(0) = −4P Q < 0 and we see that there must exist at least one ψ(r ) such that S(ψ) = 0 with 0 < ψ(r ) < 2Q, and so there exists a solution to Eq. (29). From we see that for ψ ≥ 0, S (ψ) > 0, and so the function S(ψ) is strictly increasing for all r ≥ 1. Therefore, we can conclude that S(ψ) = 0 has a unique positive solution.
The following proposition shows the dependence of the left hand side of the modulus Eq. (21) on r . (21) can be expressed as

Proposition 3 The left term of modulus equation
Proof We consider the derivatives of E (see Eq. 3) and F (see Eq. 6): and substitute into (23) We define f as and using the following relations that follow from (24), (25), (28a) and (28b): where ψ(r ) is the unique root of (29) as proved in Proposition 2, we have that Therefore (34) holds. Now that we have introduced some notation which facilitates our analysis of the imaginary roots of the characteristic equation, we separate the following analysis into two cases based on the value of r : r = 1 and r > 1. For the case r = 1, we show that there exists no possible solution for the modulus equation (18c). For r > 1, we partition the parameter space based on relations between the dimensionless parameters P and Q, and determine the conditions under which Hopf bifurcations occur within each subregion.

Bifurcation analysis for r = 1
Theorem 1 For r = 1, there exists no possible solution for the modulus equation (18c), implying that no Hopf bifurcations exist in this case.
Proof Let us first look at the right hand side of the modulus equation (21). At ω = 0, the function g(ω) (see (22)) has a removable discontinuity: lim ω→0 g(ω) = k T rbc , and the infimum of g(ω) occurs in the limit as ω → 0: Thus, g(ω ) ≥ k T rbc for each = 1, 2, . . .. Now let us look at the left hand side of the modulus equation (21). From Proposition 3 we have that since P, Q and ψ (1) are all positive.
Thus we see that the left hand side of (21) is always less than k T rbc , while the right hand side of (21) is always at least k T rbc , showing that the modulus equation (21) is never satisfied for r = 1.

Bifurcation analysis for r > 1
We now consider the case r > 1. In the following, we divide the parameter space into three regions according to the dependence of the equilibrium point (M ∞ , A ∞ ) on r . We define the nullcline where the derivative M ∞ (r ) is identically zero, and the suband supercline where this derivative is negative and positive, respectively. Within each region, we determine the conditions for which there exists pure imaginary roots to the characteristic equation (16). Finally, we show that all such imaginary roots correspond to a loss of stability of the equilibrium point via a Hopf bifurcation.
It is immediately seen that ψ(r ) = 1 is a solution of this equation. Then from Eq. (28a) we see that M ∞ (r ) = D 2Q , i.e. M ∞ (r ) is a constant. Thus, M ∞ (r ) = 0. Now assume that M ∞ (r ) = 0, i.e. we are on the nullcline. From Proposition 4 we know that ψ(r ) = 1, for all r > 1 in the nullcline, and from Proposition 2 we know that ψ(r ) is the unique root of S(ψ) = 0. Substituting ψ(r ) = 1 in S(ψ) = 0 we obtain Isolating Q, we obtain Q = 1+P 2P .
In the following, we identify the family of values of r where bifurcations can occur in the nullcline. Recall that there exists a sequence ω , = 1, 2, . . . which satisfies the argument equation stated in Proposition 1.
This implies that Hopf bifurcations can only occur at these values of r .
Proof We assume that P and Q satisfy the nullcline condition, Q = 1+P 2P , and consequently ψ(r ) = 1, for all r ≥ 1 from Proposition 4. From the expression of b(r ) in Proposition 3 we obtain that Recall that b(r ) is the left hand side of the modulus equation (21). Evaluating the modulus equation (21) at the roots ω of the argument equation (18a), we obtain k T rbc r null 2(1 + P) = g(ω ), and so r null satisfies (45).

Corollary 1
In the nullcline, the steady state of the model does not depend on the exponent r . Furthermore, the following explicit formulas hold in this region: Proof From Proposition 4, we have that ψ(r ) = 1 for all r > 1 in the nullcline. Then, from Eq. (31), we arrive at Eq. (46). From the definition of G (Eq. 27), we use Eq. (28b) and arrive at Eq. (47). From the nullcline relation Q = 1+P 2P (from Theorem 2), we isolate E max using definitions of P (Eq. 24) and Q (Eq. 25) to arrive at Eq. (48).
We have determined conditions for the existence of imaginary roots to the characteristic equation in the nullcline. In the following sections we perform a similar analysis for the subcline and the supercline.
Conversely, suppose M ∞ (r ) < 0. Then from the definition of M ∞ (r ) in Proposition 2 we have that dψ dr < 0. Equation (44) then implies that ψ(r ) > 1, which then means that S(1) < 0, which is equivalent to Q > 1+P 2P . Let us define the set of values r sub > 1 such that both the modulus equation and the argument equation are satisfied in the subcline. Recall again that there exists a sequence ω , = 1, 2, . . . which satisfies the argument equation stated in Proposition 1.

Theorem 4
In the subcline, there exists a sequence of exponents r sub for which the modulus equation (18c) and the argument equation (18a) are satisfied simultaneously. That is, there exists a unique r sub > 1 for every = 1, 2, . . . such that b(r sub ) = g(ω ).
Proof Since b(r ) > 0 and from Proposition 2 we know that ψ(r ) > 0, it follows from Eq. (34) in Proposition 3 that This provides a sharper estimate for the upper bound of ψ: From Proposition 5 we have that ψ (r ) < 0. Therefore from (50) we have that b (r ) > 0 for all r ≥ 1 and thus b(r ) will be strictly increasing for r > 1. Since ψ(r ) is strictly decreasing, we also have that the function 1 − 1+2P 4P Q ψ(r ) is strictly increasing. From here we have that where we used that ψ satisfies S(ψ) = 0 (see (29)). Therefore for all r > 1. Multiplying the inequality (51) by r and k T rbc we obtain From this we see that lim r →∞ b(r ) = +∞, implying that b(r ) may take on any value larger than b(1). Since b(r ) is strictly increasing, there exists a unique r sub > 1 for every = 1, 2, . . . such that Eq. (49) is satisfied.

Supercline
Recall that we define the supercline as the region of the parameter space for which M ∞ (r ) > 0. Proof This proof is analogous to the proof of Theorem 2 and Proposition 5. If Q < 1+P 2P then S(1) > 0 and so 1 > ψ(r ). Since ψ(r ) < 1, Eq. (44) implies that dψ dr > 0 and from the definition of M ∞ (r ) in Proposition 2 we conclude that M ∞ (r ) > 0.
In the following proposition, we divide the supercline into two subregions with different behaviors.

Proposition 7 Let us assume we are in the supercline region, that is, Q
Proof (a) From Proposition 3 we have that and from Proposition 6 we have that ψ(r ) < 1. This immediately implies that so for β < 1, this implies that lim r →∞ b (r ) = +∞.
We have identified all possible pure imaginary roots λ = iω to the characteristic equation. To ensure that the equilibrium loses stability through a Hopf bifurcation at these points, the roots must cross the imaginary axis. In the following section, we show that all imaginary roots in our system define Hopf points.

Existence of Hopf points
Recap: In Sect. 3, we defined three conditions for the existence of purely imaginary roots of the characteristic equation (16). We have shown that there is an increasing sequence of frequencies ω > 0, = 1, 2, . . . for which the argument equation (18a) is satisfied and the sign condition (18b) holds. These frequencies depend only on T p > 0, T rbc = ν max α > 0, and k > 0. Fulfilling the third condition, the modulus equation (18c), is more involved. To facilitate the analysis we introduced three dimensionless parameters: P, Q and r ≥ 1.
For r = 1, we showed that the modulus equation is never satisfied and therefore there exist no purely imaginary roots to the characteristic equation. For r > 1, we define three regions with distinct properties on the first quadrant of the (P, Q) plane, namely the nullcline, subcline and supercline.
The nullcline consists of the set of points (P, Q) belonging to the positive branch of Q = 1+P 2P with asymptotes Q = 1 2 and P = 0. In the nullcline, we establish explicit formulas for the exponents r > 1 for which the modulus equation is satisfied and therefore λ = iω are purely imaginary roots of the characteristic equation.
The subcline consists of the set of points (P, Q) such that Q > 1+P 2P . We show the existence of a unique sequence of values r > 1 such that λ = iω is a root of the characteristic equation.
The supercline consists of the set of points (P, Q) such that Q < 1+P 2P . We consider two subregions in this region. For 1+P 4P ≤ Q < 1+P 2P , we obtain a similar result as the one for the subcline for the existence of a unique sequence of values r > 1. For Q < 1+P 4P , we prove that there exist at most a finite number of solutions of the characteristic equation.
In the following, we show that all purely imaginary roots λ = iω define Hopf bifurcations. This implies that the equilibrium point (M ∞ , A ∞ ) changes stability from being spirally stable to being repellent.
Theorem 7 A Hopf bifurcation occurs for each of the pure imaginary roots of the characteristic equation, that is, there exists a small neighborhood of a critical value r = r > 1 such that when |r − r | < , < 1, the real part of the root of the characteristic equation changes sign, i.e., the root crosses the imaginary axis.
Proof We consider here only the case of the nullcline (where we have explicit expressions for r > 0). This analysis can be extended easily to the other regions.
For η = 0 and ω = ω the following equation is satisfied is the sign condition (Eq. 18b). We consider the right-hand side of Eq. (57) for w = w and |η| ≤ as a function of η: where φ = ω (T p + T rbc ).
Since L(r ) = k T rbc r 2(1+P) we can write r as a function of η We show that the behavior of r with respect to η for |η| ≤ depends on the sign of h (0). Differentiating Eq. (59) and computing h (0) we can see that if h (0) < 0 (and h (0) > 0, respectively) there exists an 0 < such that h (η) < 0 (and h (η) > 0 respectively) for |η| ≤ 0 since h(η) is continuous in a neighborhood of 0. For singular values h (0) might vanish depending on the values of ω , T p and T rbc .
Thus r (η) is strictly increasing for |η| ≤ 0 if h (0) > 0 and strictly decreasing if h (0) < 0. Therefore η changes sign when perturbing r around r = r * and so that in both cases the bifurcation at η = 0 results in a change of stability of the equilibrium.

Examples of Hopf bifurcations and simulations
In this section we provide illustrative examples of Hopf bifurcations. We use parameters (shown in Table 2) which reproduce results within physiological ranges of hemoglobin and erythropoietin (Orr et al. 1968;Miller et al. 1990;Erslev 1990Erslev , 1991. As the dynamics that results from supercritical Hopf bifurcations is oscillatory in nature, the goal of this approach is to study how the model can reproduce behavior characteristic to that of a dynamical hematological disease. As hematological pathology is often associated with a change in the turnover rate of red blood cells, we consider examples varying α in the interval [1, 5], a range larger by a factor of two than would be expected for healthy individuals. An example could be the oscillatory dynamics in levels of erythropoietin and hemoglobin often seen in patients with renal failure who are treated with erythropoiesis-stimulating agents (ESAs), which act on the longevity of red blood cells (Fishbane and Berns 2007). We then choose three different values of α from which we can calculate T rbc using expression (9). Together with the values of T p and k we use a bisection method to solve Eq. (20) for = 1 and obtain ω 1 .
Each α places us in a different region of the parameter space. In order to determine the region represented by each α, we compare P and Q, where a value of Q equal to, greater than, or less than 1+P 2P places us in the nullcline, subcline, or supercline, respectively. We compute the values of these parameters from expressions (24) and (25) using the definitions of D (see Eq. 26) and G (see Eq. 27). We calculate critical values of r where a Hopf bifurcation arises in each of the regions. We compute these critical values of r from expressions (45), (49), or (53) from Theorems 3, 4, or 5, respectively, depending on whether we are in the nullcline, subcline, or supercline.
In Table 3, we show results obtained for values of α corresponding to the nullcline, subcline, and supercline, respectively. We observe that for the given parameter values in Table 2 and any value of α, P = 0.25 remains constant and so 1+P 2P = 2.5 is constant as well. On the other hand, as Q depends on α, it varies in the interval [1. 82, 9.13] for the given interval of α, therefore placing us in different regions of the parameter space. Each pair (ω 1 , r 1 ) satisfies the argument and modulus equations (18a) and (18c), respectively, for purely imaginary roots of the characteristic equation.
As stated in Theorem 7, we expect a change of behavior in the evolution of the system's dynamics when varying the value of parameter r around the critical value r 1 . We   Table 3. We compute the solutions for a range of r values around the corresponding critical values r 1 in each region in Table 3. Figure 1 displays the values of η, the real part of the rightmost root of the characteristic equation, as a function of r . The pictures, from left to right, show how η varies for the nullcline case for r ∈ [3, 6.5] (r 1 = 4.60264), the subcline case for r ∈ [3, 6] (r 1 = 4.42738), and the supercline case for r ∈ [4, 8] (r 1 = 5.92321).
From Fig. 1, we observe in each case that for values of r < r 1 the characteristic equation is only satisfied for negative η values, implying that the equilibrium point is spirally stable. Since the roots of the characteristic equation never cross the imaginary axis for r = 1, this implies that all eigenvalues of system (10) will have negative real part for any choice of parameter values. Thus, the steady state (12) is asymptotically stable for any choice of parameters when r = 1 (Hale and Verduyn Lunel 1993). For r > r 1 , the characteristic equation is only satisfied for positive η values, implying that the equilibrium point loses stability, repelling all trajectories to a stable limit cycle. Since the slope of η as a function of r in each case is positive at the critical r 1 , we have that the Hopf bifurcations are supercritical (Diekmann et al. 1995).
Next, we illustrate this behavior by simulating the long term dynamics of the model for the presented values. We examine the evolution of the levels of hemoglobin M(t) and erythropoietin A(t) in the three regions until time t = 1,000 days for the same initial values of M(t) and A(t).
We display three figures, one for each region of the parameter space. In each, we present the long term dynamics of the model for some r < r 1 and r > r 1 . The top of  Fig. 3 Evolution of hemoglobin and erythropoietin until time t = 1,000 for the subcline case. Top Convergence to the equilibrium point for r = 3.8 < 4.42738 = r 1 . Bottom Convergence to a limit cycle for r = 4.6 > 4.42738 = r 1 from a Hopf bifurcation in the nullcline for r = 4.8. Figure 3 displays similar behavior in the subcline. The top picture shows a stable, attracting equilibrium point for r = 3.8 while the bottom picture shows trajectories which converge to a stable limit cycle with positive amplitude after the emergence of a Hopf bifurcation for r = 4.6. Similarly, in the supercline case in Fig. 4, we observe the stable equilibrium point for r = 5.2 attracting all trajectories in the top picture and the formation of a stable limit cycle resulting from a Hopf bifurcation for r = 6.2. We have given examples showing that for = 1, ω 1 > 0 gives a root of the argument equation satisfying the sign condition, and we can compute a unique exponent r 1 > 1 at which a supercritical Hopf bifurcation occurs. The real part of the root of the characteristic equation is positive for r > r 1 , and therefore the equilibrium point becomes unstable and a stable limit cycle emerges, its amplitude increasing with r . All trajectories starting near the equilibrium point are attracted to the limit cycle. In all three cases, for r < r 1 the equilibrium point is spirally stable and attracts all trajectories.

Predicting oscillatory behavior from experimental data
In the previous section, we provided examples of oscillatory dynamics for a given set of parameters within reasonable physiological ranges for humans. In this section, we provide a simple strategy for parameter estimation to fit the model to any dataset given some empirical measurements and compare predictions from our model to results from a specific experiment.
6.1 Procedure for fitting parameters of the model from experimental data Let us provide a simple strategy to estimate the parameters of our model for a given dataset. In order to recreate experimental results, we assume we are provided with values or reasonable ranges for parameters which correspond to easily measured physiological quantities. These include blood levels of erythropoietin A and hemoglobin M, the destruction rate of erythropoietin k (or its half-life, from which we can calculate k as an exponential decay constant (Bélair et al. 1995), red blood cell lifespan (corresponding to ν max and/or T rbc from which we can compute the aging rate of red blood cells α using Eq. 9), mean corpuscular hemoglobin R or red blood cell count (allowing for direct calculation of R, which is equivalent to M divided by the red blood cell count), and time spent in the precursor state T p .
The values for the parameters E max , E 50 , F max , F 50 and the exponent r which correspond to our Hill functions E and F (Eqs. 3 and 6, respectively) remain to be calculated. For our parameter estimation, we assume that the equilibrium values of hemoglobin and erythropoietin do not depend on a bifurcation parameter. This condition is fulfilled only in the nullcline-as stated in Corollary 1, the values of the steady state in this region are independent of the exponent r . In the nullcline we can take advantage of having explicit formulas  for the computation of F max , F 50 , and E max for some chosen E 50 . We can then define a set of parameters for which the values of A ∞ and M ∞ computed by our model are equivalent to the given steady states of erythropoietin and hemoglobin (or are chosen to be contained within their given ranges). Note that variations in the exponent r affect only the stability of the equilibrium values and do not change their absolute values for solutions in the nullcline.

Simulating hemolytic anemia in rabbits
We utilize the procedure outlined above to recreate pathological hemoglobin oscillations in rabbits as shown in experiments by Orr et al. (1968). In these experiments hemolytic anemia in rabbits was induced every 2-3 days via injection of an incompatible red blood cell isoantibody, creating long term oscillations in blood hemoglobin levels.
We use experimental data taken from references (Burwell et al. 1953;Hewitt et al. 1989;Bélair et al. 1995). We consider the maturation time of red blood cells in rabbits to be approximately 3 days, so we set T p = 3 days. The lifespan of red blood cells in healthy rabbits, ν max , is given to be in a range of 45-50 days. The half life of erythropoietin is 2.5 h, and from this we can easily compute k = 6.654 day −1 as an exponential decay constant. Physiological ranges of hemoglobin and erythropoietin are M ∈ [94, 174] g/l and A ∈ [10, 18] U/l, respectively. From the value of M and total red blood cell count, 3.8−7.9 × 10 12 cells/l, we compute R ≈ 4 × 10 −11 g.
In healthy animals, we assume that red blood cells survive the duration of the maximum life span, so T norm rbc = ν max and therefore α norm = 1. From the physiological ranges given, we choose M ∞ = 134 g/l and A ∞ = 14 U/l to be equilibrium values of hemoglobin and erythropoietin and setting E 50 = 80 U/l, we compute values for Hill function parameters using Eqs. (46)-(48). The complete set of model parameters for a healthy rabbit is displayed in the first column of Table 4.
We now consider the case of destabilization of the dynamics regulating red blood cell production. In this case pathological values of T rbc are calculated as the average distance between peaks in hemoglobin in one experimental rabbit (Figure 4 in Orr et al. 1968). From this, we have that T path rbc = 12.4 days and, fixing ν max ∈ [45, 50] days, we have α path ≈ 3.8. The second column of values in Table 4 shows new values of model parameters corresponding to T path rbc .  We simulate long term dynamics of the model for estimated healthy and pathological parameter values for different values of the exponent r in Hill function F. For a healthy rabbit, we show convergence to a stable equilibrium point, and in a pathological rabbit, we present an example for which the equilibrium values of A and M are stable and an example for which a supercritical Hopf bifurcation occurs.
Calculation of the critical value of the exponent r at which a supercritical Hopf bifurcation occurs is explained in detail in Sect. 4. We recall that we first need to compute ω 1 > 0 from argument equation (20) in Proposition 1 using given values of k, T p and T rbc and then use formula (45) from Theorem 3 to obtain critical exponent r 1 .
The critical value of the bifurcation parameter r for the healthy case using T norm rbc and α norm as given in the first column of Table 4 is calculated to be r norm 1 = 17.942. Figure 5 shows stable behavior for a healthy rabbit with r = 5.7 < r norm 1 = 17.942. The critical value of the bifurcation parameter r for the pathological case computed from values in the second column of Table 4 is r path 1 = 5.6237 (corresponding to ω 1 = 0.33599). Figure 6 shows 500 days of evolution for our model with parameters corresponding to rabbits with induced hemolytic anemia from initial values M 0 = 160 g/l and A 0 = 11 U/l. The top picture represents stable behavior for hemoglobin and erythropoietin for r = 5 < r While oscillatory dynamics can technically be simulated in our model for healthy rabbits, we note that r path 1 is significantly lower than r norm 1 , suggesting that it is far more likely for sick rabbits to exhibit unstable hemoglobin levels. Furthermore, though explicit measurements of r cannot easily be made, this exponent can be thought of roughly as the sensitivity of the rate of red blood cell production to erythropoietin levels. Therefore, it must be contained within some reasonably tight range, making it even more unlikely that a value higher than r norm 1 would be reached in vivo. From Eq. (45), we see explicitly that r 1 decreases with T rbc , and so for "sufficient" pathology (i.e. T rbc low enough), r 1 may enter the range of physiologically reasonable values of r , allowing for possible spontaneous oscillations in hemoglobin levels due to small natural perturbations in other parameters. This is reflected in our simulations, as r path 1 < r norm 1 . Oscillatory dynamics are induced for r = 5.7 for an anemic rabbit (the bottom plot of Fig. 6) while only stable behavior is observed for an equivalent value in healthy animals (Fig. 5).
A rule of thumb for supercritical Hopf bifurcations is that if r 1 is the value at which a supercritical Hopf bifurcation occurs, then the amplitude of oscillations is proportional to √ r − r 1 for any r close to r 1 and r > r 1 , see Strogatz (2000). In r > r 1 with an error on the order of O(r − r 1 ). The period can therefore be written as T = 2π ω 1 + O(r − r 1 ) (see Strogatz 2000). In one experimental rabbit (Figure 4 in Orr et al. 1968), hemoglobin oscillations were observed with a period of 17 days, a mean hemoglobin level of approximately 75 % of pre-injection values, and a range of approximately ±10 % of the mean. Since we consider an r = 5.7 close to r path 1 = 5.6237 in the bottom plot of Fig. 6, we can estimate the period as above by T ≈ 2π ω path 1 = 18.7 days, where ω path 1 = 0.33599 corresponds to the critical value r path 1 . While a decrease in the mean of M is not observed, the range of M values around the mean is reasonable. Furthermore, fluctuations of erythropoietin levels A are also physiologically reasonable, ranging from about ±25 % around the mean. While it is to be expected that there are some aspects of experimental results that our model cannot exactly reproduce (for instance, the drop in mean at onset of anemia), the results predicted are qualitatively sound as well as within quantitative physiologically sound ranges.

Conclusions
Periodic diseases can be simulated by systems in which the equilibrium point loses stability through a Hopf bifurcation. In a model of erythropoiesis, this would result in oscillations in erythropoietin and hemoglobin levels. Clinically, such behavior is observed in patients who use erythropoiesis-stimulating agents (ESAs) after renal failure.
The goal of this paper was to provide a clearer understanding of the dynamics of an age-structured model of erythropoiesis through a theoretical bifurcation analysis. We have reduced the model presented in Sect. 2 to a system of nonlinear ordinary differential equations with two constant delays and determined the conditions under which the unique equilibrium point can lose stability. We have performed a complete bifurcation analysis and derived analytical expressions which allow the identification of all Hopf bifurcations that occur over the parameter space. We have presented examples and shown through a numerical analysis of the roots of the characteristic equation that the Hopf bifurcations in our system are supercritical, implying long-term oscillations in hemoglobin and erythropoietin levels. Numerical simulations are in agreement with our analysis, showing the existence of oscillatory dynamics in both erythropoietin and hemoglobin levels for long time computations. Furthermore, we have also provided a simple procedure to estimate model parameters to fit experimental data and predict pathological oscillations in real systems. We have presented an example for which our model replicates oscillations in hemoglobin levels shown in rabbits after induced hemolytic anemia.