Large Amplitude and Multiple Stable Periodic Oscillations in Treatment–donation–stockpile Dynamics

A transmission–treatment–donation–stockpile model was originally formulated for the 2014–2015 West Africa Ebola outbreak in order to inform policy complication of large scale use and collection of convalescent blood as an empiric treatment. Here we reduce this model to a three dimensional system with a single delay counting for the duration between two consecutive donations. The blood unit reproduction number R 0 is calculated and interpreted biologically. Using the LaSalle's invariance principle we show that the zero blood bank storage equilibrium is globally asymptot-ically stable if R 0 < 1. When R 0 > 1, the system has a non-zero equilibrium with potential occurrence of Hopf bifurcations. The geometric approach previously developed is applied to guide the location of critical bifurcation points. Numerical analysis shows that variations of the single delay parameter can trigger bi-stable large amplitude periodic solutions. We therefore suggest that this time lag must be carefully chosen and maintained to attain stable treatment availability during outbreaks.


Introduction
Convalescent blood transfusion (also called passive immune therapy) concerns the treatment that transfuses blood from convalescent survivors to the ills.This treatment method has strong track records and wide applications not only in the pre-antibiotic era [19], but also in many modern treatment of diseases such as tetanus, hepatitis A/B, rabies, measles, complications of smallpox vaccination [16], H1N1, Spanish Influenza, severe acute respiratory syndrome (SARS), and Chikungunya [10,17,22,23].During the emergence of a new virus/bacteria disease, treatment methods are usually limited and vaccines are rarely available, hence convalescent blood transfusion is often considered as a potential treatment option, such as for H5N1, middle east respiratory syndrome (MERS), and Ebola [22,31,32].
During the 2014-2015 Ebola outbreak in West Africa, the World Health Organization (WHO) initiated an official guidance about the use of convalescent blood transfusion as the empiric therapy for Ebola virus disease (EVD) [32].This guidance provided detailed guidelines on donor selection, screening, donation and handling of blood products.In the event of large scale use of convalescent therapy, keeping the blood storage on a desired level is crucial for maintaining a positive feedback cycle: a high level of blood bank storage provides adequate treatments to increase the survival rate, those survivors become potential donors, and their donations in turn boost the blood storage.
Two modelling studies about Ebola blood transfusion were published during the outbreak [14,15], with a major difference in modelling the patient inflow rate for blood transfusion.In [14], it is assumed that any available blood will be utilized upon request, so the inflow rate of patients for blood transfusion is set to be the minimization of the available blood units and the number of untreated patients.In [15], the inflow rate function is considered as a multiplication of two factors: the rate of newly identified cases (which corresponds to disease transmission) and the fraction of patients who can receive blood transfusion (which corresponds to deployment policy).
In this paper, we detach the treatment-donation-stockpile sub-dynamics from the model in [15], with the aim of investigating the dynamics of convalescent therapy being used as a long-term empirical treatment for an endemic disease.To do this, we assume that the inflow rate of patients for blood transfusion is a function f that depends on the blood bank level at the same time.The flow diagram of the detached treatment-donation-stockpile dynamics is given in Figure 1.1, from which we have a 5-dimensional ordinary differential equation system with a single delay given below: (1.1) We note that system (1.1) can be formulated as a 3-dimensional system.By denoting D(t) := D 1 (t) + D m (t), we have In this way, equations of T, D, B are independent from D 0 , hence can be isolated from (1.1).We thus have the following three dimensional system: ) make donations at rate α and then enter the recovery class (D 0 ) and stay at least τ days until being considered again as potential donors for their next donation.We assume a loss rate ξ of all donors (D 1 , D 0 , and D m ) due to loss of immunity/contact.Donated blood units are recruited to the blood bank (B) at rate ω, and have an expiry rate δ, the level of blood bank storage in turn affects the rate of new patients admitted for transfusion treatment f (B).
For simplicity, we assume ω = α, since in Ebola convalescent blood transfusion, one effective blood donation can be used to treat one patient [32].The admission rate of patients for blood transfusion, f (B), can be decomposed and understood as a composition of two parts.Namely, which means that the inflow rate is affected by (i) a rate ηB proportional to the blood bank level, and (ii) a decreasing probability function g(B) since the efficiency of treatment administration could be reduced due to the burden of screening and matching with large stockpile amount.
In Section 2, we start with preliminary analysis of the system (1.2).In section 3, we compute the blood unit reproduction number R 0 , and prove that the zero equilibrium is globally asymptotically stable when R 0 < 1.In Section 4, we discuss the occurrence of possible bifurcations near the positive equilibrium.In Section 5, we provide a numerical example to illustrate rich dynamics including co-existence of multiple stable large amplitude oscillations.
Future challenges lie in locating the critical Hopf bifurcation values, validating the transversality conditions, and applying the global Hopf bifurcation theory [11,34] to explain the appearance of bifurcation branches, such as the analysis accomplished in [27].Although such a similar analysis is not available now for our 3-dimensional system, we show that the geometric approach developed in [4] enables one to numerically find the critical parameter values for Hopf bifurcation.

Preliminaries
We adopt the aforementioned assumption of α = ω in (1.2) and rewrite the system for all where τ > 0 is a delay and all the 6 parameters are positive, T(t) represents the population under blood transfusion treatment at time t, D(t) is the number of potential donors at time t, and B(t) denotes the number of blood bank stockpiles at time t.
For each τ > 0 we denote by where we require Based on our discussion in the introduction, we further formalize our assumptions on the patient admission rate function f as where C 1 (R + , R + ) refers to the Banach space of continuously differentiable functions, and BC(R + , R + ) refers to the space of bounded continuous functions.f ∞ denotes the upper bound for f ∈ BC(R + , R + ).
From the first equation in (2.1) and (A.2) we have Finally from the third equation of (2.1) we have Hence which completes the proof.
Furthermore, we can prove that the solution is ultimately bounded.
Proof.In this proof, we continue using the notation of m 1 introduced in (2.3).Assume for contradiction that D(t) is unbounded, we can find a t * > 0 such that Then integrating both sides of (2.

Equilibria and stability
In the following, we discuss the existence of positive equilibria of (2.1).An equilibrium (T * , D * , B * ) of (2.1) must satisfy the following where (0, 0, 0) is obviously the zero-equilibrium under any condition of parameters.Before moving into the discussion about the positive equilibria, we define the basic reproduction number R 0 where r 1 = εγ γ+µ represents the ratio of patients who survive and become potential donors after blood transfusion therapy, r 2 = α α+ξ−αe −ξτ is the expected number of donations made by each donor, and r 3 = f (0) δ+ f (0) is the probability of making use of one blood into treatment before it expires.We now state the following theorem for the existence of the equilibria.Theorem 3.1.Assume (A.2) and (A.3) hold, where Then, the following results hold.
(1) Assume R 0 ≤ 1, let (T * , D * , B * ) be a nonnegative equilibrium of (2.1).From (3.1), it is easy to see that where B * satisfies the following equation Then for any B * > 0, by (A.3) we have which means that equation (3.3) can never be satisfied, so it has no positive root when r 1 r 2 > 1 either.Hence, (0, 0, 0) is the only equilibrium of (2.1) given R 0 ≤ 1.
(2) By the proof of ( 1), it follows that if (2.1) has at least one positive equilibrium, there must be R 0 > 1.We only need to prove the sufficiency, assume R 0 > 1, since f (0) > 0 (hence r 3 < 1), it follows that r 1 r 2 > 1 and and the existence of positive equilibrium (T * , D * , B * ) follows. ( ) is onto and one to one, there exists a unique B * ∈ (0, ∞) such that (3.4) holds.It follows that the obtained (T * , D * , B * ) is the unique positive equilibrium of (2.1).
Next we show that under the assumptions (A.2) and (A.3), R 0 = 1 is the threshold for the local stability of the zero equilibrium.Theorem 3.2.Given (A.2) and (A.3), the zero equilibrium is locally asymptotically stable if R 0 ≤ 1 and is unstable if R 0 > 1.
Proof.Firstly, we prove the stability of zero equilibrium when R 0 ≤ 1.We denote z (t) := (T (t) , D (t) , B (t)) T and the system can be linearized at the equilibrium E 0 as where Compute the characteristic equation associated with the linearized system, we get where a = µ + γ, b = α + ξ, and c = δ + f (0).We rewrite (3.5) as follows, Taking the norm on both sides, we get Now we assume for the purpose of contradiction that there is an eigenvalue with positive real part Re λ > 0, then we have which contradicts to our assumption.
Moreover, we can prove that the zero equilibrium is globally asymptotically stable when R 0 ≤ 1. Theorem 3.3.Assume both (A.2) and (A.3) hold, if R 0 ≤ 1, then the zero equilibrium is globally asymptotically stable.
Proof.The preliminary results in section 2 show that R + × C + × R + is positively invariant with respect to the system (2.1) and with ϕ satisfying (A.1) being the initial condition, the solution (T(t), D t , B(t Next we define a Liapunov functional V : Now we are able to calculate the derivative of V: When Thus we are able to identify the set We next determine the largest positively invariant set in S with respect to (2.1).For a positively invariant set M ⊆ S, and for any (ψ 1 , ψ 2 , ψ 3 ) ∈ M ⊆ S being the initial condition, let (T(t), D(t), B(t)) be the unique nonnegative solution of (2.1).We have (T(t), D t , B(t)) ∈ M ⊆ S for all t > 0, that is, B(t) = 0 for all t > 0. From the third equation of (2.1) we have 0 = B (t) = αD(t), so we have D(t) = 0 for all t ≥ 0. So from the second equation in (2.1) we get D(t) = D(0) + εγ t 0 T(s)ds + αe −ξτ t−τ −τ D(s)ds = 0 for all t ≥ 0, so we have T(t) = 0 for t ≥ 0, and D(t) = 0 for t ∈ [−τ, +∞].Therefore, we have ψ 1 = 0, and ψ 2 (s) = 0 for s ∈ [−τ, 0].Then the largest positively invariant subset of S should contain all such elements, we have hence the omega limit set is not empty and is contained in M by the Liapunov-LaSalle principle [28, Theorem 5.17], so the zero equilibrium is globally asymptotically stable.

Hopf bifurcations from the positive equilibrium
In what follows, we treat the donor recovery period, τ, as the delay dependent bifurcation parameter, and consider the values of τ such that R 0 > 1.By Theorem 3.1, there exists at least one positive equilibrium, we denote it as E * := (T * , D * , B * ), and notice that the steady state E * depends on the value of τ.Next we determine critical values of τ at which Hopf bifurcations could occur near E * .
Linearizing the system (2.1) at E * , we rewrite the characteristic equation (3.5) as where P and Q are polynomials of λ of degrees 3 and 2, respectively, and )

General method revisited
We follow the method developed in [4] to determine critical parameter values of τ for Hopf bifurcation, and we denote the critical values as τ * .In other words, we look for τ * such that (4.1) has pure imaginary solutions.For convenience, we separate P and Q in terms of their real and imaginary parts: where P r , P m , Q r and Q m are the real and imaginary parts of P and Q, respectively.The overall assumptions in [4] are: (H 4 ) let I := {τ ∈ R + : there exists ω ∈ R + such that F(ω, τ) = 0}, for each single-valued continuous branch ω : From (H 4 ), for each τ ∈ I and each continuous branch ω(τ), a function θ : I → [0, 2π) is well-defined to satisfy the following condition: We are thus able to define the functions S n : I → R for each n ∈ Z + := {0, 1, 2, 3, . . .}: We present the partial result in [4] that guarantees the existence of τ * .
Finding solutions of (4.6) is equivalent to looking for the positive solutions of We can determine that A 2 = a 2 + (α + ξ) 2 + c 2 − α 2 e −2ξτ > 0, but A 0 , A 1 could have sign switches under the variation of τ.Furthermore, the derivative of ω(τ) is hard to compute explicitly as both c and f (B * ) are functions of τ.
Because of the complexity of the system that mentioned above, we will only focus on numerically locating all critical Hopf bifurcation values.To get the expression in (4.4), we compute and hence the function θ : I → R is defined as: Therefore, in the numerical simulation, we can plot the functions S n , n ∈ Z + through (4.6) and (4.8).Our numerical simulation in the next section shows that there exist two positive solutions of F(ω, τ) = 0, which leads to two distinct single-valued continuous branches of ω(τ), hence results in distinct S n functions for each n ∈ Z + .So we denote different branches of ω(τ) and S n (τ) functions as follows.

Critical bifurcation points
The global bifurcation diagram is shown in Figure 5.4: each branch connects a pair of Hopf bifurcation points.The blue part corresponds to stable branches, and the red part shows the unstable branches.We are able to identify 5 stability switches on the two branches, and label them by A-E, respectively.Point A projects to τ A = 21.1920 on τ-axis, where a simple positive Floquet multiplier moves into the unit circle through 1 as τ increases, as shown in Figure 5.5.A fold bifurcation happens at τ A .For τ ∈ (τ A , τ * 1 ), there are two limit cycles with different stability.The stable limit cycle has a large amplitude of 1400 (cases), and the system loses its stability sharply on the left side of the Hopf bifurcation critical value τ * 1 .Moreover, bi-stability occurs when τ ∈ (τ A , τ * 1 ), since the positive equilibrium is also locally stable for such τ values.As shown in Figure 5.6 for τ = 21.5, we can get solutions which converge to either the positive equilibrium or the large amplitude periodic solutions with different initial values.
Point B corresponds to τ B = 57.5233,where a pair of conjugate complex Floquet multipliers leave the unit circle as τ increases, as shown in Figure 5.7.A torus bifurcation occurs at τ B .We conjecture that there is at least one more bifurcation appearing as we move along the  first branch from point B to point (τ * 3 , 0): the quasi-periodic solution bifurcated from point B will finally disappear near the critical Hopf bifurcation point (τ * 3 , 0), since the limit cycle originated from the Hopf bifurcation is not accompanied with a quasi-periodic solution, and the point where quasi-periodic solution disappears could correspond to another critical value of torus bifurcation.However, we are unable to identify any critical points since the Floquet multipliers along the periodic solution (the Hopf branch bifurcated from (τ * 3 , 0) undergoes sharp changes with respect to τ).
Point E with τ E = 36.0604, is a flip (period doubling) bifurcation point, where we observe a simple Floquet multiplier moves into the unit circle through −1 as τ increases.Numerically, we only observe that two stable periodic solutions coexist on the interval (τ E , τ B ), neither of which is related to the flip bifurcation at point E. We set τ = 50 in Figure 5.8 and find two periodic solutions: one with period 125 days and a larger amplitude, and another with period 50 days and a smaller amplitude.
Both point C (τ C = 91.2307)and D (τ D = 98.5984) are flip bifurcation points, where we observe that a simple Floquet multiplier moves outside (at C) and inside (at D) of the unit circle through λ = −1, as shown in Figure 5.9.The Hopf bifurcation branch switches from being stable to unstable at C and becomes stable again at D. Moreover, we find that the period doubling limit cycle, the one bifurcated from point C is stable up to τ = 97.We verify that this limit cycle does have a period which is approximately twice of the period of the stable limit cycle for τ ∈ (τ E , τ C ).We show a double period solution for τ = 93 in Figure 5.10 (B), where the periodic solution has a period 181 days, and is almost twice of the period of the solution shown in Figure 5.10 (A) for τ = 90.Furthermore, a quadruple period solution is observed when τ = 97, and some chaotic behaviors are observed when τ ≈ 98, as shown it in Figure 5.11 (A) and (B), respectively.).We fix the same initial conditions for D(t) = 100, t ∈ [−τ, 0], and take different initial conditions on (A).T(t) = 500, B(t) = 500, (B).T(t) = 500, B(t) = 20, (C).T(t) = 50, B(t) = 500.

Implication on public health
We summarize our findings on long term behaviors of the system (5.2) by varying the length of donor recovery duration τ.The parameter τ can be likely adapted in an emergency situation, with a careful consideration of all public health and patient safety factors.1. τ ∈ [0, τ A ): the blood bank storage will stabilize at a high level.When the donor recovery duration is close to τ A days, there might be slight oscillations at the beginning of the dynamic (Figure 5.12).
2. τ ∈ (τ A , τ * 1 ): lengthening the donor recovery duration will result in possible large amplitude fluctuations in the blood bank level (Figure 5.13).It is most likely to have a stable amount of blood storage if the initial storage level is high enough (compare (A) and (B) in Figure 5.6), it corresponds to the fact that during outbreaks, sufficient stockpile can help building a positive feedback on a stable blood donation until the outbreak is over.Moreover, sometimes initial blood storage is not enough for a stable effect, provided there are no enough infections at the beginning (compare (A) and (C) in Figure 5.6).For diseases that persist for long time periods, sufficient initial blood storage should be accompanied with enough patients (who are future donors), so as to stabilize the blood bank level.

τ ∈ (τ *
1 , τ E ): delaying the donor recovery duration to such extent will result in fluctuations of blood bank level in the future.4. τ ∈ (τ E , τ B ): Figure 5.8 shows that the blood bank level will rely on initial data, but will definitely oscillate with two possible periodic intervals.The number of transfusion patients and blood stockpile either changes quickly with smaller amplitudes or changes slowly with larger amplitudes, depending on the initial stockpile level.
5. τ ∈ (τ B , τ C ): when donor recovery duration increases in this interval only periodic solutions with smaller amplitude can be observed (Figure 5.10 (A)).
6. τ ∈ (τ C , τ D ): extending the donor recovery duration will lead to oscillations of blood storage level with periods that are multiple magnitudes of the usual period with shorter donor recovery duration (Figure 5.10 (B), 5.11).

Discussions
Convalescent blood transfusion is believed to play an important role for treatment of emerging and re-emerging infectious diseases.Here we developed a novel treatment-donationstockpile model as an initial step of understanding the dynamic interactions of three components involved, and our model provides a framework to examine the feasibility and dynamic outcomes of large scale use and collection of convalescent blood in future outbreaks.We assumed in the numerical example that the rate of new infections decreases exponentially with respect to the present blood bank storage, and this assumption can be altered based on further implementation strategies.We also assumed that recovered patients can donate for multiple times with a donor recovery duration τ as the bifurcation parameter, and we showed that this recovery duration is a crucial factor in the maintenance of blood bank level.
The reproduction number R 0 for this system represents the average number of blood units reproduced by an initial blood unit, and it is a decreasing function of the delay parameter.It is not hard to understand that longer donor recovery period would result in smaller reproduction number, since the longer a donor waits for the next donation, the higher possibility for the donor to lost his/her immunity or contact.As shown in a numerical example, the system undergoes at least 5 stability switches including hopf, fold, torus, and flip bifurcations.Consequently, we observe coexistence of positive equilibrium and periodic solutions, multiple periodic solutions with different amplitudes and periods, double and quadruple periodic solutions, and even chaotic behavior.Our results provide important information for decision making: (i) under the optimal setting involving treatment survival ratio, donation collection frequency, and blood sample storage environment, the donor recovery duration should be carefully chosen and maintained in a range; (ii) it is crucial to optimize the stockpile storage before the outbreak and keep good track of the potential donors.
Our analysis of this proposed treatment-donation-stockpile system illustrates that a relatively simple model motivated by a significant public heath issue can exhibit complicated dynamics, and hence decision makers and the public in the epidemic region during a large scale outbreak should be informed the possibility of fluctuation of the stockpile level at different stages when use of blood transfusion service is made available to a large population.The model can also be modified to incorporate other realities during an outbreak such as nonmedical intervention (quarantine and isolation), and potential availability of other treatment and vaccine.We do hope that the model and the model-based analysis encourage the development of more effective use of convalescent blood, such as the production of human plasma derived immunoglobulin, and the identification and use of the serum from asymptomatic patients.
Finally, with the practical applications in the future emerging outbreaks of infectious diseases, the bifurcation behaviors of the treatment-donation-stockpile system is worthwhile to be examined with different possible settings of the patient inflow rate function f .In particular, the global continuation of the local Hopf branches in the example of this paper could be proved by using global Hopf bifurcation theory [11,34], with further analytical results on the non-existence of periodic solutions with certain period lengths.The difficulties for such proofs are elaborated in the appendix section.

A Appendix
In order to apply the global Hopf bifurcation theory [11,34], one thing left to prove is that all periodic solutions on the bifurcation branches have bounded finite periods.It would be helpful if one can prove that solutions with certain periods do not exist, such as in [27], the authors were able to show that the system do not obtain periodic solutions with periods of τ, 2τ, or 4τ.However, we show by numerical simulation that, our system (1.2) does have a periodic solution of period τ.
We look for a periodic solution with period τ when R 0 > 1, that is, D(t) = D(t − τ), ∀t > 0. So we try to determine if there are periodic solutions with period τ/n, n ∈ Z + of the ODE system: We keep all the parameters and setting of function f being the same with the numerical study in section 5.  Furthermore, the difficulty of determining the upper bound for periods of all periodic solutions lies in the fact that our system (1.2) with parameters in section 4 obtains large period solutions that are multiple times of τ, as shown in Figure A.2.We leave this as an open question arising from the novel treatment-donation-stockpile system for future projects.

Figure 1 . 1 :
Figure 1.1:The treatment-donation-stockpile flow diagram.Patients under blood transfusion (T) experience a death rate µ and a recovery rate γ.A percentage ε of recovered patients become potential first-time donors (D 1 ), all potential donors (D 1 and D m) make donations at rate α and then enter the recovery class (D 0 ) and stay at least τ days until being considered again as potential donors for their next donation.We assume a loss rate ξ of all donors (D 1 , D 0 , and D m ) due to loss of immunity/contact.Donated blood units are recruited to the blood bank (B) at rate ω, and have an expiry rate δ, the level of blood bank storage in turn affects the rate of new patients admitted for transfusion treatment f (B).
Next we use MATLAB package DDE-BIFTOOL v.3.1 to draw the bifurcation diagram and the global branch of each bifurcation.Figure 5.3 shows the critical values of τ for Hopf bifurcations, which coincide with the values found by the method stated in the above section.Numerically, we can see the local stability switch of the positive equilibrium from Figure 5.3: E * (τ) is locally stable for τ ∈ [0, τ * 1 ), unstable for τ ∈ (τ * 1 , τ * 4 ) and resumes its local stability for τ ∈ (τ * 4 , τ * e ).

Figure 5 . 3 :
Figure 5.3: Real parts of the eigenvalues of characteristic equation (4.1)with τ as the varying parameter.

Figure 5 . 4 :
Figure 5.4: Global bifurcation diagram with τ as the bifurcation parameter.We color the stable parts of the branches in blue, and color the unstable parts in red.

Figure 5 . 5 :Figure 5 . 6 :
Figure 5.5: Floquet multiplier of the upper bifurcation branch at τ A = 21.1920.A fold bifurcation happens at point A as a Floquet multiplier moves into the unit circle through λ = 1.

Figure 5 . 7 :Figure 5 . 8 :Figure 5 . 9 :
Figure 5.7: Floquet multipliers of the upper bifurcation branch at τ B .A torus bifurcation happens at point B as a pair of complex Floquet multipliers move out of the unit circle at τ B = 57.5233.
Figure A.1(a) shows the bifurcation diagram about the positive equilibrium E * of the system (A.1), and we can see that there does exist a periodic solution of (A.1) with period τ from the intersection of the two curves in Figure A.1 (b).