On the Backward Bifurcation of a Vaccination Model with Nonlinear Incidence

A compartmental epidemic model, introduced by Gumel and Moghadas [1], is considered. The model incorporates a nonlinear incidence rate and an imperfect preventive vaccine given to susceptible individuals. A bifurcation analysis is performed by applying the bifurcation method introduced in [2], which is based on the use of the center manifold theory. Conditions ensuring the occurrence of backward bifurcation are derived. The obtained results are numerically validated and then discussed from both the mathematical and the epidemiological perspective.


Introduction
Epidemic models have contributed greatly to get insight on transmission dynamics of infectious diseases in host populations, as well as on how an infectious disease may be managed, reduced and possibly eradicated (see, e.g., [3][4][5][6]).Two of the main aspects of modelling an infectious disease are: (i) the functional form of the force of infection, namely the function describing the mechanism of disease transmission; (ii) the description of the intervention policy to contrast the disease spread (vaccination, treatment, health campaign, etc.).
In the seventies, V. Capasso and his coworkers [7][8][9], stressed the importance to consider nonlinear incidence rates for some specific diseases.Since then, many authors have proposed several nonlinear forms for the incidence rate (see, e.g., the brief surveys contained in [10] and [11]).Many studies also deal with epidemic models with general incidence rates, see, e.g., [12][13][14].
c Vilnius University, 2011 For many diseases (pertussis, mumps, influenza, etc.) vaccination of susceptibles is an effective method for controlling the epidemic spread, and sometimes it has been a total success (smallpox, polio, measles, . . . ) [3].When it is included in mathematical compartmental models, vaccination is often represented by a linear transfer between the susceptibles and removed compartments.When vaccine is imperfect, i.e., not totally effective or does not offer permanent immunity, then a backward transfer must be considered because vaccinated individuals may return to be susceptibles or become directly infected, through the (nonlinear) transmission.When these aspects are included in the model, a rich dynamical behaviour may arise, including bistability and backward bifurcation [15][16][17].
In 2003 Gumel and Moghadas have proposed a compartmental deterministic model, which incorporates the two aspects mentioned above: a nonlinear incidence rate and an imperfect preventive vaccine given to susceptible individuals [1].Their model is interesting because the structure is simple (it is a generalization of classical SIR models) so that possible unusual dynamical behaviours may be directly lead back to the coupling of nonlinear incidence and imperfect vaccine.The authors perform a qualitative analysis and one of the main results is that the model undergoes a backward bifurcation.
The main reason to investigate the occurrence of backward bifurcations it that they play a relevant role for disease control and eradication.In fact, it is now well known that in disease transmission modelling, a classical necessary condition for disease eradication is that the basic reproductive number R 0 [18], must be less than unity.However, when a backward bifurcation occurs, an endemic equilibria may also exist for R 0 < 1.This means that the occurrence of a backward bifurcation may have important public health implications.Indeed it might not be sufficient to reduce R 0 below 1 to eliminate the disease.The basic reproductive number must be further reduced in order to avoid endemic states and guarantee the eradication.Many epidemic spread characterized by backward bifurcation may be found in the literature, for both generic and specific diseases (see, e.g., [2,15,17,19,20]).More recent contributions are given in [21][22][23].The bifurcation analysis is based on the use of the center manifold theory [2,24,25].
In [1], even if the qualitative analysis stressed the existence of backward bifurcation, the related phenomenology was not extensively discussed in terms of bifurcation theory.We think such a feature is worth to be deeply investigated.For this reason, in this paper we aim to provide a precise indication of the bifurcation thresholds and, through the bifurcation method introduced in [2], to derive conditions, expressed in terms of the parameters of the system, ensuring that either forward or backward bifurcation occurs.
The paper is organized as follows: in Section 2 we describe the model introduced in [1] and briefly discuss the qualitative analysis performed therein.In Section 3 we investigate the existence of equilibria and provide some local stability results.In Section 4 we perform the bifurcation analysis: threshold values are determined and sufficient conditions for both forward and backward bifurcation scenarios are derived.In Section 5 we present a detailed numerical verification and compare our findings with the results obtained in [1].The epidemiological implications are discussed in Section 6. Concluding remarks are given in Section 7.

The model
The compartmental deterministic model proposed in [1] is given by: where S, V and I denote the size of compartments of susceptible, vaccinated and infective individuals, respectively.The incidence rate is a saturating Michaelis-Menten functional.All the parameters are positive constants, with the following interpretation: π is the recruitment rate of susceptibles; β 1 and β 2 are the transmission probabilities of susceptibles and vaccinated individuals, respectively; c is the average number of contact partners; ξ is the vaccination rate of susceptibles; α is the therapeutic treatment rate of infected individuals; µ is the natural death.The following assumptions are also considered: (i) β 2 < β 1 , due to the fact that vaccination can reduce or eliminate the incidence of infection; (ii) the prevalent disease does not kill infected individuals; (iii) treatment does not offer permanent immunity.
In [1] the existence of endemic equilibria is performed by using a vaccination function ξ = ξ(I), given by: which can be obtained by the algebraic system of equilibria coming from (1).It can be shown that the existence and the number of endemic equilibria depends on the monotonicity of (2) in the neighbourhood of I = 0, i.e., by the sign of ξ (0): Precisely, a threshold value ξ c exists such that the following result may be stated [1]: if ξ (0) > 0, then an unique endemic equilibrium exists for ξ < ξ c , whereas two endemic equilibria may coexist if ξ > ξ c .If ξ (0) ≤ 0, then only one endemic equilibrium may exist.In other words, the model undergoes a backward bifurcation when ξ (0) > 0 and a forward bifurcation when ξ (0) ≤ 0.
However, the stability properties in the backward bifurcation scenario are only sketched.Hence, in the next section we will discuss such phenomenology in terms of bifurcation theory.

Equilibria: existence and local stability
In this section we provide some preliminary results concerning with threshold values for the existence of endemic equilbria.We begin by observing that system (1) admits the disease-free equilibrium , 0 .
Endemic equilibria E = (S * , V * , I * ) are such that , where I * is given by the real positive solutions of the equation, where Note that the coefficient A is always negative.Hence if C > 0 then ∆ = B 2 − 4AC > 0.
By applying the Descartes' rule of signs, one positive and one negative real equilibria exist, whatever is the sign of B. It follows that the coefficient C is required to be negative in the range R 0 < 1 in order to be in line with the backward bifurcation scenario.If C < 0, then by adding the conditions B > 0 and ∆ > 0, we get two positive real equilibria and the condition ∆ = 0 provides the critical value for the saddle node bifurcation, which is related with the appearence/disappearence of the two positive equilibria in the backward bifurcation framework.To get an insight on this aspect, we inspect more closely the condition ∆ > 0, observing that where ( Nonlinear Anal.Model.Control, Vol. 16, No. 1, 30-46, 2011 Under the hypothesis β 2 < β 1 , it is easy to verify that the quantity where Taking into account of (3), it follows that 1 is a negative quantity.Hence, we may conclude that in this case: (i) the inequality c > c * 2 implies the existence of two real positive equilibria; (ii) when c < c * 2 there are no feasible equilibria; (iii) c = c * 2 is the critical value of the parameter c, for the saddle-node bifurcation.
Now, we focus on the disease-free equilibrium E 0 and investigate the occurrence of the transcritical bifurcation at R 0 = 1.
The Jacobian matrix of (1) evaluated at the disease-free equilibrium E 0 is given by so that the eigenvalues λ are real and given by: Introduce now the basic reproductive number R 0 , so that: (a) if R 0 < 1, then the eigenvalues are all negatives and E 0 is locally stable; (b) if R 0 > 1, then two eigenvalues are negative and one is positive, so that E 0 is unstable (saddle point).Observing that it follows that the disease-free equilibrium E 0 is locally stable when c < c * , whereas it looses its stability when c > c * .As a consequence, the critical value c = c * is a bifurcation value.
Next step is to investigate the nature of the bifurcation involving the disease-free equilibrium E 0 at c = c * (or equivalently at R 0 = 1).This will be performed in the next section.

Bifurcation analysis
As mentioned in Section 1, when forward bifurcation occurs, the condition R 0 < 1 is usually a necessary and sufficient condition for disease eradication, whereas it is no longer sufficient when a backward bifurcation occurs.In fact, the backward bifurcation scenario involves the existence of a subcritical transcritical bifurcation at R 0 = 1 and of a saddle-node bifurcation at R 0 = R sn 0 < 1.The backward bifurcation scenario may be qualitatively described as follows.In the neighborhood of 1, for R 0 < 1, a stable disease-free equilibrium coexists with two endemic equilibria: a smaller equilibrium (i.e., with a smaller number of infective individuals) which is unstable and a larger one (i.e., with a larger number of infective individuals) which is stable.These two endemic equilibria disappear by saddle-node bifurcation when the basic reproductive number R 0 is decreased below the critical value R sn 0 < 1.For R 0 > 1, there are only two equilibria: the disease free-equilibrium, which is unstable, and the larger endemic equilibrium, which is stable.The qualitative bifurcation diagrams describing the two types of bifurcation at R 0 = 1 are depicted in Fig. 1.As a consequence, in the backward bifurcation scenario, if R 0 is nearly below unity, then the disease control strongly depends on the initial sizes of the various sub-populations of the models.On the contrary, reducing R 0 below the saddle-node bifurcation value R sn 0 , may result in disease eradication, which is guaranteed provided that the disease free equilibrium is globally asymptotically stable.Hence, determining the sub-threshold R sn 0 may have a crucial importance in terms of disease control.
In the following we will make use of Theorem A, summarized in the appendix, which has been obtained in [2] and is based on the use of the center manifold theory [25].Theorem A prescribes the role of the coefficients a and b of the normal form representing the system dynamics on the central manifold, in deciding the direction of the transcritical bifurcation occurring at φ = 0 (see appendix and the notation defined therein).More precisely, if a < 0 and b > 0, then the bifurcation is forward; if a > 0 and b > 0 then the bifurcation is backward.
We apply Theorem A to show that system (1) may exhibit a backward bifurcation when c = c * = (µ+α)µ(ξ+µ) π(β1µ+ξβ2) .First of all, observe that the eigenvalues of the matrix, are given by: Thus λ 3 = 0 is a simple zero eigenvalue and the other eigenvalues are real and negative.Hence, when c = c * (or equivalently when R 0 = 1), the disease-free equilibrium E 0 is a nonhyperbolic equilibrium: the assumption (A1) of Theorem A is then verified.Now denote by w = (w 1 , w 2 , w 3 ) T a right eigenvector associated with the zero eigenvalue λ 3 = 0.It follows: Furthermore, the left eigenvector v = (v 1 , v 2 , v 3 ) satisfying v • w = 1 is given by: The left eigenvector v is thus: The coefficient a and b defined in Theorem A, may be now explicitly computed.Taking into account of system (1) and considering only the nonzero components of the left eigenvector v, it follows that: where In view of ( 5) and ( 6), we then get: where The coefficient b is always positive so that, according to Theorem A, it is the sign of the coefficient a -and hence the sign of the quantity a 0 -which decides the local dynamics around the disease-free equilibrium for c = c * .
In the following, we investigate the role specifically played by vaccination (ξ), treatment (α) and transmission (β 1 , β 2 ) parameters in the occurrence of backward or forward bifurcations.For our convenience, we introduce the parameter β = β 1 − β 2 .In this way, where Both A and C are positive coefficients, whereas B is always positive provided that α < µ+2π.On the contrary, if α > µ+2π, then B > 0 if and only if α−3µ−4π .Now, by applying Theorem A, we may conclude that if a 0 = Aβ 2 + Bβ + C < 0, a backward bifurcation occurs, whereas if a 0 > 0, system exhibits a forward bifurcation.We discuss these two cases separately.Backward bifurcation.It is easy to verify that a 0 < 0 if ∆ > 0 and B < 0. Hence, the following conditions allow the existence of a backward bifurcation at c = c * : where β a and β b are the real positive roots of the equation Forward bifurcation.In this case it is easy to check that a 0 is always positive if (i) α < 3µ + 4π or (ii) α > 3µ + 4π and ξ < ξ = 4µ(µ+π) α−3µ−4π .Furthermore, if (iii) ∆ > 0 and B < 0, then a 0 > 0 if and only if β < β a or β > β b .More explicitly, each of the following conditions allows the existence of a forward bifurcation at c = c * : or or

Numerical verification
In this section we aim to provide a numerical verification of the above results and to show their agreement with the ones obtained by Gumel and Moghadas, regarding the existence of the endemic equilibria and their stability properties.We also put into evidence that the bifurcation analysis may offer a unifying perspective of their results.We consider the same parameter values used in [1]: π = 700, β 2 = 0.000003, µ = 0.03.Then, in order to satisfy the inequalities (8) we take α > 3µ + 4π = 2800.09,i.e., α = 3150; these parameter values imply ξ = 0.24, so that we choose ξ = 0.  Setting β = 0.000097, system (1) exhibits a backward bifurcation at c = c * = 17718.9187;the saddle-node bifurcation value is c * 2 = 17716.8243,Fig. 3(a).By choosing instead β < β a or β > β b , a forward bifurcation occurs as stated by condition (11).Fig. 3(b) illustrates the phenomenology related to the case β = 0.00002 < β a whereas Fig. 3(c) shows the case β = 0.00012 > β b .
We can now compare our findings with the results obtained in [1], where two sets of parameter values are considered, leading to different dynamical situations.
Our bifurcation analysis allows to state that both these parameter sets fall into the same bifurcation framework.In fact, in both these cases, condition (9) is verified (i.e., α = 0.7 < 3µ + 4π = 2800.09),allowing a forward bifurcation to occur.Moreover the qualitative phenomenology related to this bifurcation framework can give direct indication on the stability properties of the equilibria involved.For example, for the parameter set II, the forward bifurcation value is c = c * = 2.6472.Hence, being c = 4 > c * , it is possible to conclude that, for this value of c, the disease free-equilibrium is unstable and the endemic equilibrium E * is the only attractor for the system.Fig. 5 depicts this situation.We are thus in perfect accordance with the existence and stability results presented in [1].The above examples also suggest how the bifurcation analysis may be considered, in this context, a useful and rich qualitative tool, since different cases related to specific parameter sets (as parameter set I and parameter set II), may be read in the unified perspective of the same bifurcation scenario.

Epidemiological aspects of the results
The bifurcation analysis performed in the previous section deserves a deep discussion under the epidemiological point of view.In fact, conditions ( 8)- (11), which are summarized in Table 1, allow to answer two important questions.First, which are the mechanisms responsible for the occurrence of a certain bifurcation scenario, i.e., forward or backward?Second, inside each scenario, how can we act in way to eradicate the disease or eventually to prevent disease outbreaks?Table 1.Synoptic table of conditions ( 8)- (11).
We begin by answering at the former.Our results show that the occurrence of a certain kind of bifurcation critically depends on the interplay among three biological mechanisms explicitly included in the model: vaccination, treatment and transmission.
Vaccines, as disease-control tools, aim to reduce the degree of susceptibility of a healthy susceptible individual against a particular infection.Vaccination is mainly characterized in terms of (i) vaccination rate, i.e., the rate at which susceptible individuals are immunized (ii) efficacy, i.e., the percentage of susceptibles left unprotected even if vaccinated.Both these features are included in model (1): the vaccination rate is represented by the parameter ξ whereas the efficacy of vaccine is hidden in the preliminary assumption β 2 < β 1 , so that β 2 = σβ 1 where 0 < σ < 1 is the efficacy parameter.
Treatment of infected is represented through the parameter α and, according to the model assumptions, treatment of infected individuals does not offer here permanent immunity.Transmission mechanisms are related to the parameters β 1 and β 2 .
Our analysis shows that -inside a fixed bifurcation framework -disease control (eradication or prevention) may be obtained by suitably acting on the contacts c with infected individuals.In model (1) the dependence of the basic reproductive number R 0 on the parameter c may be obtained by observing that so that R 0 is an increasing function of c: reducing the average number of contacts c, causes the decreasing of R 0 .Hence, in order to eradicate the disease or to prevent disease outbreaks, one has to identify the critical contact value c * necessary to reduce R 0 below 1.However, in the presence of multiple endemic equilibria (as it may happen in the backward bifurcation scenario), this might not be sufficient to eliminate the disease.For this reason, it is partic-ularly important to elucidate those mechanisms that can control or avoid such backward situations.
Our analysis suggests that a key role in distinguishing between forward or backward scenario, is played in primis by the therapeutic treatment α of infected individuals.If the treatment α is kept below a certain threshold α * = α * (µ, π) = 3µ + 4π, a forward bifurcation occurs at c = c * , or equivalently at R 0 = 1.
If the therapeutic treatment α is above such threshold α * , then it is the vaccination rate level ξ to prescribe the occurrence of a backward or a forward scenario.If ξ is below the critical value ξ, then a forward bifurcation occur at R 0 = 1.
Finally, if the treatment α is above α * and the vaccination rate ξ is higher than ξ, then the transmission parameters decide between the two bifurcation scenario.More precisely, if β = β 1 − β 2 is sufficiently small or sufficiently large (i.e., β < β a or β > β b ), then a forward bifurcation occurs.On the contrary, if β = β 1 − β 2 remains limited in a certain range (i.e., β a < β < β b ), then a backward bifurcation will occur.
In all the cases of forward bifurcation scenario, disease control may be performed by reducing c below the transcritical bifurcation value c = c * , so as to lower R 0 below 1.
In the case of backward bifurcation scenario, disease control is not so smooth.Two threshold values have in fact to be kept in mind: the critical value c = c * (equivalent to R 0 = 1) and the critical value c = c * 2 (equivalent to R 0 = R sn 0 ) related to the saddle-node bifurcation which causes the disappearing of the two (respectively stable and unstable) endemic equilibria in the range R 0 < 1. Disease eradication may be obtained by lowering c below c * 2 so that R 0 < R sn 0 : in this way, the disease free equilibrium will be the only attractor for the system.
In conclusion, the above results indicate that a backward bifurcation scenario may be encouraged by the standing of all these features: (a) sufficiently high therapeutic treatment, i.e., higher than a certain threshold, only depending on natural death and recruitment; (b) sufficiently high vaccination rate; (c) difference of transmission terms, β 1 − β 2 bounded in a certain range.
Under the point of view of the disease control campaigns, public policy makers might thus move on two fronts.First, try to avoid or to prevent the dangerous backward scenario, by which disease endemicity could persist also under the classic threshold R 0 < 1.At this aim, they may act by keeping suitably low the therapeutic treatments of infected individuals or the vaccination of the susceptible individuals.In this case the disease control may straightly follow the classic road of reducing R 0 below 1.
Second, if the backward scenario cannot be avoided, public policy makers have to be particularly careful since simply keeping R 0 < 1 may still result in endemicity of the disease: in fact, in this case a reduction of the basic reproduction number below the subthreshold R sn 0 is required in order to avoid the dangerous range [R sn 0 , 1].This goal may be obtained by reducing c under the threshold c * 2 , through aimed public education programs which, affecting individual behaviors, can act to properly reduce the contacts and hence the disease transmission.
(iii) a > 0, b < 0. When φ < 0, with |φ| 1, x = 0 is unstable and there exists a locally asymptotically stable negative equilibrium; when 0 < φ 1, x = 0 is stable and a positive unstable equilibrium appears; (iv) a < 0, b > 0. When φ changes from negative to positive, x = 0 changes its stability from stable to unstable.Correspondently, a negative unstable equilibrium becomes positive and locally asymptotically stable.
Remark.Taking into account of Remark 1 in [2], we observe that if the equilibrium of interest in Theorem A is a non negative equilibrium x 0 , then the requirement that w is non negative is not necessary.When some components in w are negative, one can still apply the theorem provided that w(j) > 0 whenever x 0 (j) = 0; instead, if x 0 (j) > 0, then w(j) need not to be positive.Here w(j) and x 0 (j) denote the jth component of w and x 0 respectively.

Fig. 1 .
Fig. 1.Qualitative bifurcation diagrams for the forward (a) and backward (b) bifurcations respectively.The bifurcation parameter is the basic reproductive number R0.The solid lines (-) denotes stability; the dashed line (--) denotes instability.
6.According to (8), we have a backward bifurcation at c = c * , for β a = 0.000037 ≤ β ≤ β b = 0.000106.This backward scenario implies a range in the bifurcation parameter for which both the disease free equilibrium and the endemic equilibrium with the larger I value are stable.In terms of the bifurcation parameter c, bistability occurs for c ∈ B, where B = [c * 2 , c * ].It is interesting to observe how the size of such bistability region B depends on the value of β.As shown in Fig. 2, for the chosen parameter values, the maximum of the function f (β) = c * (β) − c * 2 (β) occurs for β ≈ 0.00006, which in turn corresponds to largest range B. In order to make B as small as possible, one should take, for example, β ≈ β b .