Methods for deriving necessary and sufficient conditions for backward bifurcation

Backward bifurcation has significant implications for disease control. Deriving necessary and sufficient conditions for backward bifurcation is of paramount importance to understand the reasons for its occurrence and devise effective control strategies. In this paper, we review the methods that lead to necessary and sufficient conditions for backward bifurcation in infectious disease models. We review separately the methods that apply to ODEs and methods that apply to PDEs. We further propose a new method, applicable to both ODEs and PDEs. We illustrate the methods on three examples: a novel ODE model of cholera with vaccination, a PDE version of the cholera model with vaccination, and on an eight equation model of dengue, taken from the literature.


Introduction
From the time of the first definition of R 0 and the first mathematical models of infectious diseases, it has been believed that R 0 is a threshold, such that if R 0 < 1, the disease will die out, while if R 0 > 1 the disease will persist. In the 1990s, bifurcation analysis near the critical value of R 0 revealed that the bifurcation of the endemic equilibria might not be supercritical but subcritical in R 0 , thus allowing for endemic equilibria when R 0 < 1 [10]. Hadeler and Castillo-Chavez, finding backward bifurcation in a core group model with vaccination, realized and discussed the public health implications of this phenomenon [8]. Since then backward bifurcation has been a very studied phenomenon.
In bifurcation theory, a forward bifurcation is a particular type of local bifurcation where the system transitions from one fixed point, which is typically stable, to two nonnegative fixed points, one stable, one unstable, when a bifurcation parameter, usually taken in mathematical epidemiology as the basic reproduction number R 0 crosses a critical value, say R 0 = 1 (see Figure 1, left). Backward bifurcation is a particular type of local bifurcation where the system transitions from two nonnegative fixed points, one stable, one unstable, to one fixed point, typically unstable, when a bifurcation parameter, usually taken in epidemiology as the basic reproduction number R 0 , crosses a critical value, R 0 = 1 (see Figure 1, right). One important caveat accompanying the definition is that R 0 is a compound parameter, and it is not necessary to use R 0 as a bifurcation parameter. One can  use any of the model parameters that make up R 0 as a bifurcation parameter. One has to be careful in this case, however. Parameters p which increase R 0 , that is R 0 (p) is an increasing function, produce forward and backward bifurcations as illustrated in Figure 1. Parameters p which decrease R 0 , that is R 0 (p) is a decreasing function, produce inverted forward and backward bifurcations as illustrated in Figure 2. The critical bifurcation value of the parameter p, denoted by p * is the value such that R 0 (p * ) = 1. To simplify matters and avoid confusion in this article we will assume that p is a parameter so that R 0 increases as p increases.
A quick search in MathSciNet revealed that there are more than 200 articles including the term 'backward bifurcation' in the context of natural sciences. Many articles examine potential causes of backward bifurcation. Those causes include vaccination, exogenous reinfection (e.g TB), multi-group models, treatment, public health education campaigns, superinfection, disease-induced mortality in vector-borne models with standard incidence and many others [7]. In general, it has been noticed that multiple groups of individuals with different susceptibility to the disease lead to backward bifurcation [1,4,11]. This is in particular the underlying mechanism in all of the above causes, except the backward bifurcation in vector-host models generated by standard incidence.
Presumably the presence of backward bifurcation has several implications for the disease and its control. First, in the presence of backward bifurcation, if R 0 ≈ 1, the number of infectives jumps significantly from its initial values and equilibrates at high equilibrium numbers, rather than increase gradually as R 0 increases above one [1]. Secondly, even if R 0 < 1, the disease may persist if the initial number of infected is sufficiently large. For this reason, when backward bifurcation is present for the range of R 0 given by R * < R 0 < 1, it is usually recommended that additional control measures should be applied to lower R 0 below its critical value R * . Rather than decreasing R 0 below R * , in some cases, as [6] suggests, it is possible and perhaps more efficient to eradicate the disease whose R 0 < 1 with temporary control measures which may be applied for some period of time.
A key component of understanding how to avoid backward bifurcation in the first place is to obtain a necessary and sufficient condition for that phenomenon to occur and investigate the reason for its occurrence as well as potential control strategies that may change parameters in a beneficial way with the goal to avoid the backward bifurcation. Necessary and sufficient conditions for backward bifurcation also allow us to find examples where such a phenomenon occurs for plausible parameters [13]. This is important as experimental evidence for the occurrence of backward bifurcation is scant and it is not clear how to test or observe it in infectious disease systems.
Deriving necessary and sufficient conditions for backward bifurcation is not new. The first method we list here seems to date back to 1990s [9] which can also be applied to PDEs [15]. One more recently developed method, based on the central manifold theorem, has acquired a significant popularity [2] due to its easy applicability to very complex models.
In this paper, our main objectives are (1) to review existing methods for deriving necessary and sufficient conditions for backward bifurcation in ODE and age-structured PDE models and (2) to propose a new method, called sensivities method, that may be used both for ODEs and PDEs. The new sensitivities method is applicable for simple and complex models and allows for computation of some steps using computer algebra systems (such as Mathematica). We illustrate the methods on a new ODE model of cholera with vaccination, on an age-since-infection structured version of the cholera model, and on a dengue model with vaccination taken from [5]. In the next section, we introduce the methods for ODEs, in section 3 we introduce the methods for PDEs, and in section 4, we apply the new method for ODEs and PDEs to the existing model in [5].

Methods for deriving necessary and sufficient conditions for backward bifurcation in ODE models
The early papers on backward bifurcation used techniques from bifurcation theory to detect the phenomenon. These approaches often did not produce clear conditions on the parameters for the backward bifurcation to occur. Having necessary and sufficient conditions on the parameters for the backward bifurcation to arise facilitates the analysis and the design a specific example, but more importantly, it delineates the specific conditions under which the backward bifurcation exists and may suggest what public health efforts can be used to eliminate it or reduce its depth. In this subsection, we list the methods that lead to necessary and sufficient conditions for backward bifurcation and are applicable to ODE models. We illustrate the methods on the following example of cholera dynamics with vaccination: Example 2.1: To introduce the model let S be the number of susceptibles, V the number of vaccinated individuals, I is the number of infected individuals, R is the number of recovered and B is the amount of bacteria in the environment. While cholera can transmit directly, the main route of transmission is environmental and we incorporate only environmental transmission in the model. The model with vaccination takes the form Vaccination is an existing control measure of cholera, but few cholera models with vaccination have been considered (see [3]). The authors in [3] determine that there is a unique endemic equilibrium which is globally stable. Although full analysis of model (1) is not of interest to us, we show that this model exhibits backward bifurcation. The system (1) has a disease-free equilibrium E 0 = (S 0 , V 0 , 0, 0, 0) where The control reproduction number is given by The parameters of the model are listed in Table 1. The following result is not hard to establish: Proposition 2.1: If R(ψ) < 1 then the disease-free equilibrium is locally stable. If R(ψ) > 1 then the disease-free equilibrium is unstable.

General description of method
This method consists in solving the system of equilibria and expressing all dependent variables in terms of one whose value is zero at the disease-free equilibrium (usually the infective individuals, I). That produces an equation for the equilibrium value I: F(I) = 0. The equation F(I) = 0 should not have zero as a solution (that is, we must cancel I while simplifying). The equation F(I) = 0 depends on the parameters of the model, including the basic or controlled reproduction number R that provides a criterion for stability of the disease-free equilibrium. The idea is to treat I as dependent on a parameter p of the model. We can then rewrite the equation for I in the form: The parameter p is typically either R or one of the parameters in R. We recall that for simplicity we assume that R increases with the parameter p. At the critical value of R = 1, which gives a critical value of p = p * , the bifurcation is backward in p (and R) iff the derivative of I with respect to p is negative in the case when R is increasing with p. That is, Using implicit differentiation, we obtain Hence, It appears that this method was first applied to ODEs in [9].

Application of the sign of the derivative method
In what follows we apply the sign of derivative method to the model of cholera with vaccination (1) (Example 1). The system of equilibria takes the form: The equilibrium values of the dependent variables (S, V, I, R, B) do not depend on time any more, but depend on the parameters of the model.
In this application of the method, the role of I in the general description of the method will be played by the bacteria in the environment B, while the role of the bifurcation parameter p will be played by β.
First, we solve the system for the equilibria to obtain an equation in B, F(B) = 0, that gives the endemic equilibria. With Substituting V in the equation for I, we obtain: In principle Equation (7) can be rewritten as a quadratic equation and conditions for two positive equilibria can be posed. We discuss this approach in the next subsection. The approach described below is more general and can be applied independently of the complexity of the equation F(B) = 0. Notice that zero is not a solution of this equation. We differentiate implicitly Equation (7) with respect to β Next, we evaluate at the disease-free equilibrium. We obtain where Q = (1 + σ ψ/μ). The condition ∂B/∂β < 0 becomes, after canceling β/D: where  (1) with respect to parameter β. Parameters are as in Table 1 except D = 10 7 .
The necessary and sufficient condition (10) immediately suggests that there will be no backward bifurcation if q = 0, in other words if w = 0 or γ = 0, that is, if there is no recovery or if immunity does not wane. Somewhat more elaborate computations show further that if ψ = 0 backward bifurcation also does not occur. Thus, the backward bifurcation in this model is a result of the combined effect of recovery, waning immunity and vaccination. Figure 3 shows that backward bifurcation occurs in model (1) for plausible parameter values. Furthermore, the figure shows that backward bifurcation does not occur for all vaccination levels. Since vaccination causes the backward bifurcation, conceivably at low enough vaccination levels, the phenomenon may not occur. However, the figure also shows that backward bifurcation would not occur at high enough vaccination levels. Thus it is a phenomenon associated with intermediate levels of vaccination.

General description of the method
In this method, the equations for the endemic equilibria are reduced to a quadratic equation in the key-dependent infected variable, say, I. For R = 1, this equation should have a non-zero constant term. Then one applies conditions for multiple positive equilibria for quadratic equations. In particular, if we have the quadratic equation:  We add that the quadratic equation has exactly one positive root if and only if a 0 < 0. This approach can be applied for some of the simpler models. A simpler approach that is applicable to more complex models that lead to a quadratic equation is to observe that if we assume that the reproduction number is equal to one, then typically a 0 = 0 and the quadratic equation has a zero root. If the other root is positive, then the quadratic equation has two positive roots when R ≈ 1. . This needs to be verified separately. (2) The conditions on the parameters for the bifurcation to be backward may not the same as the condition(s) for two endemic equilibria to exist. This may occur in the case when the dependent variable (I in our case) is constrained from above, and the upper branch does not return completely to the bifurcation point (see Figure 4). In this case, Theorem 2.1 will not work, regardless of the fact that backward bifurcation occurs and there is a parameter range with 2 positive equilibria. Theorem 2.1 appears to have been first proposed by Brauer [1].

Application of the method to Example 1
We apply Theorem 2.1 to model (1). We again use as key-dependent variable the free bacteria in the environment B. Using computer algebra system, we rewrite Equation (7) as a quadratic equation: Since a 2 < 0, Theorem 2.1 implies that the condition for backward bifurcation is a * 1 > 0. Keeping in mind that we simplify a * 1 > 0 to the following inequality: Further rewriting leads to the condition: where and Q = (1 + σ ψ/μ). Of course, one could have stopped at (12). We continued rewriting so we can compare the results from the different methods. This method only works if the equation for the key infectious quantity can be reduced to a quadratic equation. This is obviously the case with model (1).

Sensitivities method
We introduce here a new method, called sensitivities method. This new method is an extension of the derivative method. This extension is again based on imposing a condition on the derivative with respect to one of the parameters; however, the approach for computing that derivative is distinct, presumably easier and better adaptable for use of computer algebra systems.

General description of the method
We split the variables into infected states y and remaining variables x. Thus the system can be written as where x ∈ R n and y ∈ R m . We start again from the system for the equilibria where x ∈ R n and y ∈ R m . The method can be described in general terms with the following steps: (1) We look at the equations for the equilibria, thinking that all dependent variables depend on a given parameter p. We determine one of the infectious variables as the 'key', say y 1 and we divide the equilibrial Equations (16) for the infectious variables by y 1 . Hence, we have the following augmented system for the equilibria: where z ∈ R m−1 is given by z k = y k /y 1 for k = 2, . . . , m.

Remark 2.1:
There are two solutions bifurcating from the critical point (p * , 0) in the (p, y 1 ) plane: the disease-free equilibrium y 1 ≡ 0 and the endemic equilibrium y * 1 > 0. Thus, dy 1 /dp has two values: one is zero and the other one is non-zero. The normalization by y 1 is necessary to avoid computing the tangent to the disease-free equilibrium.
(2) Compute the solution of system (17) corresponding to the disease-free equilibrium, that is with y 1 = 0. This will produce the usual values of the disease-free equilibrium, as well as values, typically nonzero, of the quotients z k = y k /y 1 .

Application of the method to Example 1
We illustrate the method on system (1). We first decide on a parameter. We take p to be β in the system. The role of our key dependent variable y 1 will be played by B. We consider the system for equilibria and we augment it by the system obtained from dividing the equations for the infectious variables I and B by the key variable B. We obtain: where I/B is treated as a monolitic variable and can be denoted by z. The augmented disease-free equilibtium is given by (S 0 , V 0 , 0, 0, 0, δ/η).
From the last equation, we can see that and replace it in the previous equation. However, this cannot be always done, so we continue by differentiating the full system. We will denote by primes the derivatives with respect to β. At the same time we are evaluating at the disease-free equilibrium and β * : We solve the equations for the proportions to obtain one equation for the key variable y 1 or in our case B. We obtain from the last two equations: Expressing I in terms of B from Equation 5 in system (23), and canceling B , we see that equation three in system (23) is trivial, because R(ψ) = 1 (since β = β * ). Solving the first five equations for S and V in terms of B we obtain: Hence, Substituting in Equation (24) we have Solving for B , we obtain B = A n /A d , where We see that the numerator is positive. The necessary and sufficient condition for backward bifurcation then becomes A d to be negative, that is, After some rewriting, we can see that this condition is equivalent to the one we had before, namely, and Q = (1 + σ ψ μ ).

Center manifold method
This method was first formulated as a Theorem in [2], although the result has its roots in prior research, most notably [19]. In application of the main Theorem, use of computer algebra systems in some of the computations is also possible. The method allows for easy application and guarantees the necessary and sufficient condition for backward bifurcation. This form of the Theorem only applies to ODEs.

General description of the method
We recall the theorem here [2]: Theorem 2.2 (Castillo-Chavez and Song): Consider the following general system of ODEs with a parameter φ.
where 0 is an equilibrium point of the system, that is f (0, φ) ≡ 0 for all φ. Assume The local dynamics of the system around 0 is completely determined by the sign of a and b: When φ < 0 with |φ| 1, 0 is locally asymptotically stable and there exists a positive unstable equilibrium; when 0 < φ 1, 0 is unstable and there exists a negative and locally asymptotically stable equilibrium.
In practice the following two observations are important [14]: (1) The equilibrium 0 is the disease-free equilibrium, φ is one of the parameters in the reproduction number, and the critical value of φ is the value of that parameter that makes the reproduction number equal to one. (2) Since the disease-free equilibrium has positive entries, the right eigenvector w need not be nonnegative. Components of the right eigenvector w that correspond to positive entries in the disease-free equilibrium could be negative. However, components that correspond to zero entries in the disease-free equilibrium should be nonnegative. For further details, please see Remark 1 in [2].

Application of method to Example 1
We apply the Center manifold method to system (1). We use β as the bifurcation parameter φ. The Jacobian computed at the disease-free equilibrium is If β * is the value of β for which R(ψ) = 1, then it is not hard to check that the above Jacobian has a simple eigenvalue of zero and all other eigenvalues are negative. Computing left v and right w eigenvectors of the above matrix we have: Above q = wq 0 as before. The left eigenvector has all zero components except v 5 = 1 and v 3 = η/(μ + γ ). Thus in the sums for a and b the only non-zero terms are the ones that correspond to k = 3,5. From system (1), we see that From (33) we see that the second partial derivatives of f 5 are zero. Hence, the only contributing terms are those corresponding to f 3 , the right hand side of the equation for I in system (1). Using a computer algebra system, we get: The quantity b consists of one term and is clearly positive. The quantity a takes the form: That leads to the following condition for backward bifurcation: which is equivalent to the ones previously obtained.

Methods for deriving necessary and sufficient conditions for backward bifurcation in age-structured models
Virtually all of the methods listed in the above section can be applied to age-structured PDE models, except Theorem 2.2 which is specifically designed for ODEs.

Example 3.1:
To illustrate the application of the methods to PDEs, we recast model (1) into an age-since-infection structured PDE. In particular, if τ is the age-since-infection of infectious individuals, and θ is the age of the bacteria in the environment, we denote by i(τ , t) the density of the infectives and by u(θ , t) the density of the bacteria in the environment. With these notations in addition to the ones for the (1), the age-structured model becomes: The parameters have the same meanings as in Table 1; however, β and δ are functions of θ, while γ and η are structured by τ . B in this model denotes the total concentration of bacteria in the environment This system has a disease-free equilibrium E 0 = (S 0 , V 0 , 0, 0, 0) where We define the probability of survival of infected individual in the infectious class and the probability of survival of the bacteria in the environment: The control reproduction number is given by (derivation can be found in the appendix): As before the following result is not hard to establish: Proposition 3.1: If R(ψ) < 1 then the disease-free equilibrium is locally stable. If R(ψ) > 1 then the disease-free equilibrium is unstable.
To apply the methods we have to select bifurcation parameter and key infectious quantity. This model has only a few parameter constants. We could use ψ as a bifurcation parameter. However, to stay consistent with the previous discussion, we will write β(θ) = ββ 0 (θ ) where sup β 0 (θ ) = 1. We will useβ as a bifurcation parameter. As a key infectious quantity we will use u(0). Although the sensitivities method can be applied directly to the system for the equilibria, it seems easier if methods are applied to the system for the equilibria where the ODEs in (37) have been solved: and replaced in the remaining equations. Thus we obtain the system: where the following notations have been used:

Sign of derivative method
To apply the sign of derivative method, we have to select bifurcation parameter and key infectious quantity. We will useβ as a bifurcation parameter. As a key infectious quantity we will use u(0). To apply the sign of the derivative method, we need to obtain an equation in u(0). From the equations in u(0) and i(0), we have where S and V can be expressed in terms of u(0) from the remaining equations: and q = w (μ + w) .
Equation (41) becomes: Differentiating with respectβ and evaluating at the disease-free equilibrium and the critical value ofβ, we get: Hence, u (0) < 0 if and only if (after cancelingβ/D): The sign of the derivative method was first applied to age-since-infection structured models in [15] where it was used to establish the presence of backward bifurcation in a model with two exposed classes. In age-structured models backward bifurcation was first found in [11] where the quadratic equation method was used. This method can be applied to chronological-age structured models as well as diffusion models.

Sensitivities method
The new sensitivities method can also be applied to the age-structured models. Again the "key" quantity is u(0) so we augment system (38) with the equations for the proportions of the infectious variables, i(0), u(0) in this case, with the 'key' quantity u(0).
where we treat i(0)/u(0) as a monolitic variable, which we may denote with z. The augmented disease-free equilibrium is (S 0 , V 0 , 0, 0, 0, 1/ ). We differentiate with respect toβ and evaluate at the disease-free equilibrium and the critical value ofβ. The system for the sensitivities takes the form: From the last two equations, we have It remains to determine S and V in terms u (0). Solving the equations for R , S and V we get: where q = w (μ + w) . Thus, Substituting in Equation (46) and solving for u (0), we obtain u (0) = A N /A D where Hence, the condition u (0) < 0 becomes:

Applications of the methods to an existing model
Mathematical models of infectious diseases are becoming more and more complex. Methods to detect backward bifurcation need to be capable to address that complexity. In this section, we test the methods on an existing example taken from [5]. The paper introduces several models, studies backward bifurcation and the reasons for it. We test the methods for deriving necessary and sufficient conditions for backward bifurcation on model (21) in the paper. The model describes dengue transmission with vaccination.
Example 4.1: Model (21) in [5] consists of eight equations. The dependent variables are susceptible humans and vectors S H , S V , exposed humans and vectors E H , E V , infectious where For the description of parameters, we refer to Table 2 and [5]. The disease-free equilibrium is given by E 0 = (S * H , P * H , 0, 0, 0, S * V , 0, 0) where

Centre manifold method
The authors of [5] apply the Centre manifold method [2] and obtain a rather hard to interpret condition for backward bifurcation. The condition does guide them to obtain an example of backward bifurcation. We refer to [5] for the details of the application.

Sensitivities method
We apply here the Sensitivities method to complement the results in [5]. We begin with the equations for the equilibria. For models with standard incidence and disease-induced death rate, it is often convenient to first normalize the system of equations for the equilibria. We do that before applying the sensitivities method. We denote the normalized susceptible humans and vectors by s H = S H /N H , s V = S V /N V , exposed humans and vestors by e H = E H /N H , e V = E V /N V , infectious humans and vectors by i H = I H /N H , i V = I V /N V , vaccinated and recovered humans by p H and r H . The forces of infection take the form: where we recall that The system of equilibrium equations for the proportions becomes: We can apply the sesitivities method directly to this system but to simplify matters and reduce the number of equations, we first eliminate e H and e V by expressing them in terms if i H and i V .
The system for equilibria becomes: Figure 5. The function F(ξ ) plotted against ξ gives the dependence of condition (60) on vaccination rate ξ . For small vaccination rates (ξ ≈ 10) backward bifurcation becomes more likely. No amount of vaccination can eliminate the backward bifurcation with these parameters (see [5], Figure 4 for parameter values). However, the condition is very sensitive to δ H . Reducing δ H may eliminate the backward bifurcation, at least for some vaccination levels.

Discussion
In this paper, we consider four methods for deriving necessary and sufficient conditions for backward bifurcation: The sign of derivative method, the quadratic equation method, the sensitivities method and the centre manifold method. Each of these methods can be applied for ODEs and all but the Center Manifold Method can be applied also for age-structured PDEs.
The sensitivities method is a new method that we propose here for the first time. Unfortunately, at this time no rigorous proof exists for the method, but applications to multiple examples suggest that the method works well, if applied correctly. The intuition behind the sensitivities method is that it combines the sign of derivative method with sensitivities-type of approach for computing the derivative with respect to a specified parameter p.
For each specific model, the four methods will produce equivalent necessary and sufficient conditions for backward bifurcation. This is illustrated in the paper with a novel model of cholera that serves as an application example for the methods. The simplest method of the four is the quadratic equation method. If it works for a specific model, it should be preferred; however, this is the only method that may not be applicable to all models exhibiting backward bifurcation. In terms of ODEs, the most often preferred method is the Castillo-Chavez and Song method. We demonstrate on a specific example, taken from the literature, that the novel sensitivities method may produce simpler and more easily interpretable conditions for backward bifurcation than the Castillo-Chavez and Song method (also termed Center manifold method). Unfortunately at this point, the author does not believe there is a way to tell which method will produce the result faster just by looking at the model. One has to try each method to determine which one is easiest in each specific scenario, and of course the chosen method will certainly depend on the personal preferences of the people applying it.
The main purpose of this paper is to collect, describe and illustrate methods for deriving conditions for backward bifurcation both for ODEs and age-structured PDEs. As DDEs are a special case of the age-structured PDEs, we surmise that the methods for age-structured PDEs will be applicable for DDEs; however, we leave the examples in this direction for future articles. Because the present paper is a methods paper, we offer little biological interpretation, and the one included in the paper is only related to the derived necessary and sufficient conditions for backward bifurcation. Full analysis and biological interpretation of the novel cholera model is not included because we believe it is beyond the scope of this paper.
Substituting the equation for w(0) from the second equation in the first equation, and canceling y(0), we obtain the following characteristic equation η(τ )e λτ π(τ ) dτ .

R(ψ)
is defined as the value of the righthand side for λ = 0: