Global analysis of an epidemic model with vaccination.

In this paper, an epidemic dynamical model with vaccination is proposed. Vaccination of both newborn and susceptible is included in the present model. The impact of the vaccination strategy with the vaccine efficacy is explored. In particular, the model exhibits backward bifurcations under the vaccination level, and bistability occurrence can be observed. Mathematically, a bifurcation analysis is performed, and the conditions ensuring that the system exhibits backward bifurcation are provided. The global dynamics of the equilibrium in the model are also investigated. Numerical simulations are also conducted to confirm and extend the analytic results. © Korean Society for Computational and Applied Mathematics 2017.

Introduction

Mathematical models have become important tools in analyzing the spread and control of infectious diseases [2]. Based on the theory of Kermack and Mckendrick [19], the spread of infectious diseases usually can be described mathematically by compartmental models such as models (where S represents the class of the susceptible population, E is the exposed class in the latent period, I is infectious class, R is the removed class, which has recovered with temporary or permanent immunity). In recent years, a variety of compartmental models have been formulated, and the mathematical analysis of epidemiology models has advanced rapidly, and the analyzed results are applied to infectious diseases [2, 18, 32]. Vaccination campaigns have been critical in attacking the spread of infectious diseases, e.g., pertussis, measles, and influenza. The eradication of smallpox has been considered as the most spectacular success of vaccination [44]. Although vaccination has been an effective strategy against infectious diseases, current preventive vaccine consisting of inactivated viruses do not protect all vaccine recipients equally. The vaccine-based protection is dependent on the immune status of the recipient [2, 32]. For example, influenza vaccines protect 70–90% of the recipients among healthy young adults and as low as 30–40% of the elderly and others with weakened immune systems (such as HIV-infected or immuno-suppressed transplant patients) (see, [14, 30, 44]).

Since vaccination is the process of administering weakened or dead pathogens to a healthy person or animal, with the intent of conferring immunity against a targeted form of a related disease agent, the individuals having the vaccine-induced immunity can be distinguished from the recovered individuals by natural immunity. Thus, vaccination can also be considered by adding some compartment naturally into the basic epidemic models. Over the past few decades, a large number of simple compartmental mathematical models with vaccinated population have been used in the literature to assess the impact or potential impact of imperfect vaccines for combatting the transmission diseases [1, 3, 11, 16, 20, 21, 23, 31, 43, 45]. In some of these studies (e.g., papers [16, 31, 43]), authors have shown that the dynamics of the model are determined by the disease’s basic reproduction number . If the disease can be eliminated from the community; whereas an endemic occurs if . Therefore, if an efficient vaccination campaign acts to reduce the disease’s basic reproduction number below the critical level of 1, then the disease can be eradicated. While in other studies, such as Alexander et al. [1] and Arino et al. [3], they have shown that the criterion for is not always sufficient to control the spread of a disease. A phenomenon known as a backward bifurcation is observed. Mathematically speaking, when a backward bifurcation occurs, there are at least three equilibria for in the model: the stable disease-free equilibrium, a large stable endemic equilibrium, and a small unstable endemic equilibrium which acts as a boundary between the basins of attraction for the two stable equilibria. In some cases, a backward bifurcation leading to bistability can occur. Thus, it is possible for the disease itself to become endemic in a population, given a sufficiently large initial outbreak. These phenomena have important epidemiological consequences for disease management. In recent years, backward bifurcation, which leads to multiple and subthreshold equilibria, has been attracting much attention (see, [1, 3, 4, 6, 11, 16, 17, 20, 21, 23, 24, 33, 34, 37, 40]). Several mechanisms with vaccination have been identified to cause the occurrence of backward bifurcation in paper [33].

In this paper, we shall investigate the effects of a vaccination campaign with an imperfect vaccine upon the spread of a non-fatal disease, such as hepatitis A, hepatitis B, tuberculosis and influenza, which features both exposed and infective stages. In particular, we focus on the vaccination parameters how to change the qualitative behavior of the model, which may lead to subthreshold endemic states via backward bifurcation. Global stability results for equilibria are obtained. The model constructed in this paper is an extension of the model in paper [31], including a new compartment for the latent class (an important feature for the infectious diseases eg. hepatitis A, hepatitis B, tuberculosis and influenza) and the disease cycle. It is one of the aims of this paper to strengthen the disease cycle to cause multiple endemic equilibria.

The paper is organized as follows. An epidemic model with vaccination of an imperfect vaccine is formulated in Sect. 2, and the basic reproduction number, and the existence of backward bifurcation and forward bifurcation are analyzed in Sect. 3. The global stability of the endemic equilibrium is established in Sect. 4. The paper is concluded with a discussion.

The model and the basic reproduction number

In order to derive the equations of the mathematical model, we divide the total population N in a community into five compartments: susceptible, exposed (not yet infectious), infective, recovered, and vaccinated; the numbers in these states are denoted by S(t), E(t), I(t), R(t), and V(t),  respectively. Let . The flow diagram of the disease spread is depicted in Fig. 1.

Fig. 1

Flowchart diagram for model (2.1)

All newborns are assumed to be susceptible. Of these newborns, a fraction of individuals are vaccinated, where . Susceptible individuals are vaccinated at rate constant . The parameter is the rate constant at which the exposed individuals become infectious, and is the rate constant that the infectious individuals become recovered and acquire temporary immunity. Finally, since the immunity acquired by infection wanes with time, the recovered individuals have the possibility of becoming susceptible again. is the transmission coefficient (rate of effective contacts between susceptible and infective individuals per unit time; this coefficient includes rate of contacts and effectiveness of transmission). Since the vaccine does not confer immunity to all vaccine recipients, vaccinated individuals may become infected but at a lower rate than unvaccinated (those in class S). Thus in this case, the effective contact rate is multiplied by a scaling factor (, where describes the vaccine efficacy, represents vaccine that offers 100% protection against infection, while models a vaccine that offers no protection at all). It is assumed that the natural death rate and birth rate are and the disease-induced death rate is ignored. Thus the total population N is constant. Since the model consider the dynamics of the human populations, it is assumed that all the model parameters are nonnegative.

Thus, the following model of differential equations is formulated based on the above assumptions and Fig. 1,with nonnegative initial conditions and . System (2.1) is well posed: solutions remain nonnegative for nonnegative initial conditions. We illustrate here that there are limiting cases in system (2.1): if , the vaccine is perfectly effective, and , there is no vaccination, system (2.1) will be reduced to the standard SEIRS model in [28]; if and the limit , system (2.1) will be equivalent to an SVIR model in [31]. If we let and , system (2.1) can be reduced to an SVEIR epidemic model in [16], where authors aim to assess the potential impact of a SARS vaccine via mathematical modelling. To explore the effect of the vaccination period and the latent period on disease dynamics, an SVEIR epidemic model with ages of vaccination and latency are formulated in paper [10]. In papers [10, 16, 28, 31], authors have shown that the dynamics of the model are determined by the disease’s basic reproduction number . That is, the disease free equilibrium is globally asymptotically stable for ; and there is a unique endemic equilibrium which is globally asymptotically stable if . If and limit , system (2.1) will be reduced into an SIV epidemic model in [36], where authors investigate the effect of imperfect vaccines on the disease’s transmission dynamics. In [36], it is shown that reducing the basic reproduction number to values less than one no longer guarantees disease eradication. In this paper, we show that if a vaccination campaign with an imperfect vaccine and the disease cycle is considered, a more complicated dynamic behavior is observed in system (2.1). For example, the backward bifurcation occurs in system (2.1). In the following, first, it is easy to obtain that the total population N in system (2.1) is constant. To Simplify our notation, we define the occupation variable of compartments S, E, I, V,  and R as the respective fractions of a constant population N that belong to each of the corresponding compartments. We still write the occupation variable of compartments as S, E, I, V and R,  respectively. Thus, it is easy to verify thatis positively invariant and globally attracting in . It suffices to study the dynamics of (2.1) on . Thus, system (2.1) can be rewritten as the following system:In the case , system (2.3) reduces to an SEIRS model without vaccination [28], where , is considered as the basic reproduction number of the model. The classical basic reproduction number is defined as the number of secondary infections produced by a single infectious individual during his or her entire infectious period. Mathematically, the reproduction number is defined as a spectral radius (which is a threshold quantity for disease control) that defines the number of new infectious generated by a single infected individual in a fully susceptible population [39]. In the following, we shall use this approach to determine the reproduction number of system (2.3). It is easy to see that system (2.3) has always a disease-free equilibrium,Let System (2.3) can be rewritten aswhereThe Jacobian matrices of and at the disease-free equilibrium are, respectively,where, is the next generation matrix of system (2.3). It follows that the spectral radius of matrix isAccording to Theorem 2 in [39], the basic reproduction number of system (2.3) isThe basic reproduction number can be interpreted as follows: A proportion of of exposed individuals progress to the infective stage before dying; represents the number of the secondary infection generated by an infective individual when he or she is in the infectious stage. Those newborns vaccinated individuals have generated the number of the secondary infection. Average vaccinated individuals with vaccination rate have generated the fraction of the secondary infection.

Equilibria and bifurcations

Now we investigate the conditions for the existence of endemic equilibria of system (2.3). Any equilibrium (S, V, E, I, R) of system (2.3) satisfies the following equations:From the second and third equation of (3.1), we have Since , this equation can be true only for ; hence, there exists no endemic equilibrium for . For , the existence of endemic equilibria is determined by the presence in (0, 1] of positive real solutions of the quadratic equationwhere,From (3.2) and (3.3), we can see that the number of endemic equilibria of system (2.3) is zero, one, or two, depending on parameter values. For (the vaccine is totally effective), it is obviously that there is at most one endemic equilibrium () in system. From now on we make the realistic assumption that the vaccine is not totally effective, and thus .

We notice that if , then we haveSince all the model parameters are positive, it follows from (3.3) that . Furthermore, if , then . SinceThus, is a continuous decreasing function of for , and if , then and . Therefore, it follows that P(I) of Eq. (3.2) has a unique positive root for .

Now we consider the case for In this case, , and . From (3.3), it is easy to see that is an increasing function of . Thus, if , then for . Thus, P(I) has no positive real root which implies system have no endemic equilibrium in this case. Thus, let us consider the case . In this case, let . It is obvious that if , then . Notice that is an linear increasing function of . Thus, there is a unique such that and thus Since is a quadratic function of with positive coefficient for , has a unique root in Thus, for we have , , and for Therefore, P(I) has two possible roots and system (2.3) has two endemic equilibria , ) for From the above discussion, we have for , and for . Therefore, it follows that system (2.3) has no endemic equilibria for .

If , we have . In this case, system has a unique endemic equilibrium for and no endemic equilibrium for

Summarizing the discussion above, we have the following Theorem:

Theorem 3.1

If , system (2.3) has a unique endemic equilibrium ; If there exists , system (2.3) has two endemic equilibria , for and has no endemic equilibria for ; If , system (2.3) has a unique endemic equilibrium for and no endemic equilibrium for

According to Theorem 2 of van den Driesche and Watmough [39], we have the following result.

Theorem 3.2

The disease-free equilibrium is locally asymptotically stable when and unstable when .

In the following, we first give a global result of the disease-free equilibrium of system (2.3) under some conditions.

Theorem 3.3

If , is globally asymptotically stable in the feasible positively invariant region.

Proof

Consider the following Lyapunov functionalBy directly calculating the derivative of L along system (2.3) and notice that , thus, we haveIt is easy to verify that the maximal compact invariant set in is when . The global stability of follows from the LaSalle invariance principle [22].

From the above discussion, we know that system (2.3) may undergo a bifurcation at the disease-free equilibrium when . Now we establish the conditions on the parameter values that cause a forward or backward bifurcation to occur. To do so, we shall use the following theorem whose proof is found in Castillo-Chavez and Song [5], which based on the use of the center manifold theory [15].

For the following general system with a parameter .Without loss of generality, it is assumed that is an equilibrium for system (3.4) for all values of the parameters , that is, for all .

Theorem 3.4

Assume that:Let be the kth component of f andThen the local dynamics of system (3.4) around are totally determined by a and b.

is the linearization matrix of system (3.4) around the equilibrium with evaluated at 0. Zero is simple eigenvalue of A and all other eigenvalue of A have negative real parts;

Matrix A has a (non-negative ) right eigenvector and a left eigenvector v corresponding to the zero eigenvalue.

When with , is locally asymptotically stable and there exists a positive unstable equilibrium; when is unstable and there exists a negative and locally asymptotically equilibrium;

When , with , is unstable; when is locally asymptotically stable and there exists a negative unstable equilibrium;

When , with , is unstable and there exists a locally asymptotically stable negative equilibrium; when is stable and a positive unstable equilibrium appears;

When changes from negative to positive, changes its stability from stable to unstable. Correspondingly, a negative unstable equilibrium becomes positive and locally asymptotically stable.

Now by applying Theorem 3.4, we shall show system (2.3) may exhibit a forward or a backward bifurcation when . Consider the disease-free equilibrium and choose as a bifurcation parameter. Solving givesLet denote the Jacobian of the system (2.3) evaluated at the DFE with . By directly computing, we haveLetIt is easy to obtain that has eigenvalues given byThus, is a simple zero eigenvalue of the matrix and the other eigenvalues are real and negative. Hence, when , the disease free equilibrium is a non-hyperbolic equilibrium. Thus, assumptions (A1) of Theorem 3.4 is verified. Now, we denote with , a right eigenvector associated with the zero eigenvalue .

Thus,Thus, we haveThe left eigenvector satisfying is given byFrom the above, we obtain thatLet a and b be the coefficients defined as in Theorem 3.4.

Computation of . For system (2.3), the associated non-zero partial derivatives of f (evaluated at the DFE , ) are given bySince the coefficient b is always positive, according to Theorem 3.4, it is the sign of the coefficient a, which decides the local dynamics around the disease-free equilibrium for . If the coefficient a is positive, the direction of the bifurcation of system (2.3) at is backward; otherwise, it is forward.

Thus, we formulate a condition, which is denoted by Thus, if holds, we have , otherwise, .

Summarizing the above results, we have the following theorem.

Theorem 3.5

If holds, system (2.3) exhibits a backward bifurcation at (or equivalently ). Otherwise, system (2.3) exhibits a forward bifurcation at (or equivalently when ).

Remark 1

From Theorem 3.5, it can follows that the occurrence of either a backward or forward bifurcation may be expected. In fact, in system (2.3), let It is easy to verify that the condition is satisfied. By applying Xpp plot software and choosing the above parameters, we can describe the backward bifurcation diagram of system (2.3) (see, Fig. 2). Let all parameter values be same as in Fig. 2, except is changed as 0.2. The condition is not satisfied and we have forward bifurcation diagram Fig. 3 at . So it is clear that there is one threshold value of say such that backward bifurcation occurs of and forward bifurcation occurs if . Both of these bifurcation diagrams are obtained by considering as bifurcation parameter and later it is plotted with respect to .

Fig. 2

The backward bifurcation diagram for model (2.3)

Fig. 3

The forward bifurcation diagram for model (2.3)

Global stability of the endemic equilibrium

In this section, we shall investigate the global stability of the unique endemic equilibrium for Here we shall apply the geometric approach [25, 27, 38] to establish the global stability of the unique endemic equilibrium. In recent years, many authors [3, 16, 26, 28, 29] have applied this method to show global stability of the positive equilibria in system. Here, we follow the techniques and approaches in paper [3, 16] to investigate global stability of the endemic equilibrium in system (2.3). Here, we omit the introduction of the general mathematical framework of these theorems and only focus on their applications.

In the previous section, we have showed that if , system (2.3) has a unique endemic equilibrium in . Furthermore, implies that the disease-free equilibrium is unstable (Theorem 3.2). The instability of , together with , implies the uniform persistence of the state variables. This result can be also showed by using the same arguments from Proposition 4.2 in [27] and Proposition 2.2 in [29]. Hence, there exists a constant such that any solution of of system (2.3) with the initial conditions satisfiesThus, we first give the following result:

Proposition 4.1

System (2.3) is uniformly persist in for .

To prove our conclusion, we set the following differential equationwhere , is open set and simply connected and .

Letwhere, P(x) be a nonsingular matrix-valued function, which is on and is the derivative of P(x) in the direction of the vector field f(x). is also matrix, the second additive compound of the Jacobian matrix is the measure with respect to a vector norm . The following result comes from Corollary 2.6 in paper [25].

Theorem 4.1

Suppose that is simply connected and is a compact set which is absorbing with respect to system (4.1); For some matrix P, there exists a positive constant such that for all . Then the unique equilibrium in system (4.1) is global asymptotically stable.

From Proposition 4.1, it is easy to verify that the condition in Theorem 4.1 holds. Therefore, to prove our conclusion, we only verify that in Theorem 4.1 holds. According to paper [35], the Lozinski measure in Theorem 4.1 can be evaluated as follow:where is the right-hand derivative.

Now we state our main result in this section.

Theorem 4.2

Suppose that the parameters in system (2.3) satisfy the following inequalitiesThen the unique equilibrium in system (2.3) is globally asymptotically stable for .

Let , where and . Then, the Jacobian matrix of system (2.3) can be written asThe second additive compound [25](see, “Appendix”) of Jacobian matrix is the matrix given byLetSet , where is the derivative of P in the direction of the vector field f. Thus, we haveFrom (2.3), we haveThus, we obtain thatwhereAs in [3, 16], we define the following norm on :where , with components andand letNow we demonstrate the existence of some such thatBy linearity, if this inequality is true for some z, then it is also true for . Similar to analyzing methods in paper [3, 16], our proof is subdivided into eight separate cases, based on the different octants and the definition of the norm (4.4). To facilitate our analysis, we use the following inequalities:for all

Case 1. Let , and

Then we have and . Taking the right derivative of , we haveSince , and , thus, we obtain Case 2. Similarly, it is easy to verify that Eq. (4.6) also holds for and when

Thus, if we require that holds, then the inequality (4.5) holds for case 1 and case 2.

Case 3. Let , and Thus, we have and So,we haveUsing the inequalities , from the above, we obtain that Case 4. By linearity, Eq. (4.7) also holds for and when

Thus, if we require that holds, then the inequality (4.5) holds for case 3 and case 4.

Case 5. Let , and Thus, we have and . By directly calculating, we obtain thatUsing the inequalities , we have Case 6. By linearity, Eq. (4.8) also holds for and when

Thus, if we require that holds, then the inequality (4.5) holds for case 5 and case 6.

Case 7. Let , and Thus, we have and . Thus, we haveUsing the inequalities and , we have Case 8. By linearity, Eq. (4.9) also holds for and when

Thus, if we require that holds, then the inequality (4.5) holds for case 7 and case 8.

Therefore, from the discussion above, we know that if inequalities (4.3)hold, then there exists such that for all and all nonnegative S, V, E and I. All conditions in Theorem 4.1 can be satisfied when inequalities (4.3) hold. Therefore, by Theorem 4.1, we can determine that if inequalities (4.3) hold, then the unique endemic equilibrium of system (2.3) is globally stable in for .

Remark 2

In Sect. 3, we have shown that system (2.3) exhibit a backward bifurcation for . As stressed in [3], for cases in which the model exhibits bistability, the compact absorbing set required in Theorem 4.1 does not exist. By applying similar methods in [3], a sequence of surfaces that exists for time and minimizes the functional measuring surface area may be obtained. Therefore, the global dynamics of system (2.3) in the bistability region can be further investigated as it has been done in paper [3].

Discussion

In this paper, an epidemic model with vaccination has been investigated. By analysis, it is showed that the proposed model exhibits a more complicated dynamic behavior. Backward bifurcation under the vaccination level conditions, and bistability phenomena can be observed. The global stability of the unique endemic equilibrium in the model is demonstrated for . Note that the model (2.3) can be solved in an efficient way by means of the multistage Adomian decomposition method (MADM) as a relatively new method [8, 9, 12, 13]. The MADM has some superiority over the conventional solvers such as the R-K family. To illustrate the various theoretical results contained in this paper, the effect of some important parameter values on the dynamical behavior of system (2.3) is investigated in the following.

Now we consider first the role of the disease cycle on the backward bifurcation. If , [i.e., the disease cycle-free in model (2.3)], then the expression for the bifurcation coefficient, a, given in Eq. (3.5) reduces toThus, the backward bifurcation phenomenon of system (2.3) will not occur if This is in line with results in papers [16, 31], where the disease cycle-free model (2.3) has a globally asymptotically stable disease-free equilibrium if the basic reproduction number is less than one.

Differentiating a, given in Eq. (3.5), with respect to givesHence, the bifurcation coefficient, a is an increasing functions of . Thus, the feasibility of backward bifurcation occurring increases with disease cycle.

Now we consider the role of vaccination on the backward bifurcation. Let , then the expression for the bifurcation coefficient, a, given in Eq. (3.5), is reduces toThus, the backward bifurcation phenomenon of system (2.3) will not occur if (i.e., the model (2.3) will not undergo backward bifurcation in the absence of vaccination). This is also in line with results in paper [26], where the vaccination-free model (2.3) has a globally asymptotically stable equilibrium if the basic reproduction number is less than one. Furthermore, the impact of the vaccine-related parameters () on the backward bifurcation is assessed by carrying out an analysis on the bifurcation coefficient as follows. Differentiating a, given in Eq. (3.5), partially with respect to givesThus, the backward bifurcation coefficient, a is a decreasing function of the vaccination rate . Hence, the possibility of backward occurring decreases with increasing vaccination rate ( i.e., vaccinating more susceptible individuals decrease the likelihood of the occurrence of backward bifurcation).

Differentiating the bifurcation coefficient a, given in Eq. (3.5), partially with respect to giveswithThus, the bifurcation coefficient, a is a decreasing function with respect to . That is, the likelihood of backward bifurcation occurring decreases with increasing vaccine efficacy. Let By direct calculating, it is easy to verify that is negative and also condition is satisfied. Figure 4 depicts the backward bifurcation occurring phenomena with lower vaccine efficacy with ; Fig. 5 depicts the likelihood of backward bifurcation occurring with higher vaccine efficacy .

Fig. 4

The backward bifurcation diagram with respect to vaccine efficacy

Fig. 5

The backward bifurcation diagram with respect to vaccine efficacy

In addition, it is obvious that our expression for the basic reproduction number in system (2.3), i.e.,is independent of the loss rate of immunity . From the above analysis, we have found that the dynamics of the model are not determined by the basic reproduction number, and the phenomena of the backward bifurcation in system may occur. Moreover, it is found that the occurrence feasibility is increasing with the loss rate of immunity .

From the following expression,it is easy to see that the policy of vaccinations with imperfect vaccines can decrease the the basic reproduction number . Thus, the imperfect vaccine may be beneficial to the community. This is also a positive point, sice it is know that the use of some imperfect vaccine can sometime result in detrimental consequences to the community [3, 20].

At last, we must point out that although the system (2.3) with (2.2) is well posed mathematically, we acknowledge the biological reality that the fraction of the constant total population which occupies a compartment can only be within the subset Q of rational values within , and furthermore only within a sub-subset of values within Q belonging to n / N where n belongs to the integers . In addition, we also point out that the analysis of the model (2.1) may become somewhat different if disease fatalities and more complex vital dynamics are included, in particular, if the population size is no longer constant. In the future, we may investigate many various modeling possibilities to simulate a real world biological process based on model (2.1). On the other hand, we note that the population in our model (2.1) is assumed to be homogeneously mixed. In fact, different individual may have different number of contacts. Thus, a complex network-based approach on diseases transmission may be closer to a realistic situation [7, 41, 42]. In the future, we shall investigate dynamics of the proposed model based on a complex network.