Equilibrium Analysis of a Yellow Fever Dynamical Model with Vaccination

We propose an equilibrium analysis of a dynamical model of yellow fever transmission in the presence of a vaccine. The model considers both human and vector populations. We found thresholds parameters that affect the development of the disease and the infectious status of the human population in the presence of a vaccine whose protection may wane over time. In particular, we derived a threshold vaccination rate, above which the disease would be eradicated from the human population. We show that if the mortality rate of the mosquitoes is greater than a given threshold, then the disease is naturally (without intervention) eradicated from the population. In contrast, if the mortality rate of the mosquitoes is less than that threshold, then the disease is eradicated from the populations only when the growing rate of humans is less than another threshold; otherwise, the disease is eradicated only if the reproduction number of the infection after vaccination is less than 1. When this reproduction number is greater than 1, the disease will be eradicated from the human population if the vaccination rate is greater than a given threshold; otherwise, the disease will establish itself among humans, reaching a stable endemic equilibrium. The analysis presented in this paper can be useful, both to the better understanding of the disease dynamics and also for the planning of vaccination strategies.


Introduction
Yellow fever (YF), a hemorrhagic fever caused by a Flavivirus, family Flaviviridae [1,2], is characterized by fever, chills, loss of appetite, nausea, muscle pains particularly in the back, and headaches [3]. There are more than 200,000 infections and 30,000 deaths every year [3]. About 90% of YF cases occur in Africa [4], and a billion people live in an area of the world where the disease is common [3]. It also affects tropical areas of South America, but not Asia [3,5,6]. The number of cases of yellow fever has been increasing in the last 30 years [3,7], probably due to fewer people being immune, more people living in cities, people moving frequently, and changing climate [3]. The origin of the disease is Africa, from where it spread in South America through the slave trade in the 17th century [8,9].
The yellow fever virus was the first human virus discovered [10], and its family comprises approximately 70 viruses [2], most of which are transmitted by arthropod insects (hence the name arthropod borne viruses or arboviruses).
A safe and effective vaccine against yellow fever exists and some countries require vaccinations for travelers [3]. In rare cases (less than one in 200,000 to 300,000 doses), the vaccination can cause yellow fever vaccine-associated viscerotropic disease (YEL-AVD), which is fatal in 60% of cases, probably due to the genetic morphology of the immune system. Another possible side effect is an infection of the nervous system, which occurs in one in 200,000 to 300,000 cases, causing yellow fever vaccine-associated neurotropic disease (YEL-AND), which can lead to meningoencephalitis, fatal in less than 5% of cases [6]. In some rare circumstances, however, the fatality rate of vaccine induced diseases 2 Computational and Mathematical Methods in Medicine can reach alarming proportions, as observed recently by Mascheretti et al. [11], who found 1 death per million doses applied in a Southeastern Brazilian region.
In this paper, we propose an equilibrium analysis of a dynamical model of yellow fever transmission in the presence of a vaccine. Such a kind of analysis can be useful, both to the better understanding of the disease dynamics and also for the planning of vaccination strategies.

Model Formulation
The mathematical model described below addresses the transmission dynamics of an infectious agent in a homogeneous population in the presence of an imperfect vaccine. We consider a nonlinear system of ordinary differential equations involving the human and the vector-mosquitoes and their eggs-populations. The term "eggs" also includes the intermediate stages, such as larvae and pupae. It is also worth highlighting that the model proposed here is based on previous papers [12,13], and we updated the model originally developed by Amaku et al., 2013 [14, 15], for the purpose of investigating the impact of vaccination on population. By including vaccination, in particular vaccines which may have serious adverse effects, the model may help the designing of realistic (from the cost and logistic point of view) vaccination strategies.
All variables and parameters in the human system will carry the subscript , while those in the vector system will carry one of the subscripts (mosquitoes) or (eggs). In our model the total human population, denoted by , is split into four subclasses which are susceptible humans ( ), vaccinated humans ( ), infected humans ( ), and recovered (and immune) humans ( ), so that = + + + . The total vector population, which is formed by both total mosquitoes population, denoted by , and the total eggs population, denoted by , are split into susceptible mosquitoes ( ), infected and latent mosquitoes ( ), infected and infectious mosquitoes ( ), and noninfected eggs ( ), so that = + + and = . A flow diagram of the model is depicted in Figure 1, and the associated variables and parameters are described in Tables 1 and 2, respectively (values from references [16][17][18][19][20][21][22]).
The model supposes a homogeneous mixing of human and mosquito population based on the idea that the mosquito has a human biting habit, so that each mosquito bite has an equal probability of transmitting the virus to the susceptible human in the population or acquiring infection from an infected human. The equations are derived based on the fact that, in presence of the yellow fever in the population, both mosquitoes and humans can infect each other upon contact. While an infected mosquito remains infected until death, it is assumed that infected humans can recover from the disease (see [23]). We define a logistic recruitment rate of humans, mosquitoes, and eggs, and all new born humans and newly emerged mosquitoes are susceptible (no vertical transmission; see [23]). Susceptible humans become infected through the bite by an infected mosquito and the susceptible mosquitoes become latent infected as result of  biting infectious humans. Upon acquiring infection, the susceptible individuals move into the infected compartment. The incidence of new infections is given by the standard incidence (see [24, pp. 602]). Deaths can occur amongst the human population, mosquitoes, and eggs, naturally. In contrast, in the presence of the yellow fever, the human population can either die due to the additional effects of the disease or recover. It is also assumed that recovered human individuals acquire immunity against reinfection, so that they do not acquire yellow fever for a second time.
Although there is a vaccine for yellow fever, it is expected that it is imperfect; that is, it does not offer 100% protection against infection in all population. Thus, it is instructive to assess the potential impact of an imperfect yellow fever vaccine. After the duration of protection wanes down, the vaccinated individuals, , move to the susceptible class, , and they may then acquire a new infection. Hence, the vaccinated population is decreased by the waning of vaccine-induced immunity and by natural death.
Combining the above formulation and assumption, it follows that the model for the transmission dynamics of the yellow fever disease in the presence of an imperfect vaccine is given by the following system of nonlinear ordinary differential equations: with the conditions = + + + , = + + , and = and the initial conditions ≥ 0, and (0) ≥ 0. Note that the transmission from mosquitoes to humans is given by ( / ) (also called "force of infection") and from humans to mosquitoes is given by ( / ). This means that the transmission from mosquito to humans depends on the number of infective mosquitoes, but the transmission from human to mosquitoes depends on the density of infective humans.
Since system (1) models human, mosquito, and eggs populations, it is assumed that all variables in the system are nonnegative. This assumption yields the epidemiologically feasible domain: Since the right-hand sides of equations of system (1) and their partial derivatives are continuous in Ω, we will use the techniques described in [21] that there exist a unique solution It can also be verified that the given initial conditions make sure that ≥ 0. Thus, the total population remains positive and bounded for all finite time > 0. Similar arguments can be applied to both mosquito and eggs equations with corresponding expressions. Therefore, all solutions of the model with initial conditions in Ω remain in Ω for all ≥ 0, the region Ω is positively invariant with respect to model (1), and its solutions are considered epidemiologically and mathematically well posed in Ω.

The Existence of Equilibria
Our next result concerns the existence of equilibrium points of system (1) that are biologically feasible. Thus, we will find the equilibrium points of system (1) in the region Ω by setting right hand side of all equations in it as equal to zero. First of all, we will seek the conditions for the existence of the equilibria of system (1) which are biologically feasible.
From the second and fourth equations of (1) with the right-hand side equal to zero, it can be seen that the equilibrium points must satisfy, respectively, the following relations: where Therefore, * > 0 is always satisfied. Substituting (3) and (4) into the third equation of system (1), we obtain * = ( Hence, * > 0 is also always satisfied. From eighth and ninth equations of system (1), we obtain * = * , From tenth equation of system (1), we get either with thres = /( + ). From expression (11) it follows that * > 0 ⇐⇒ < thres .
Later we will see that for * = 0 only the trivial equilibrium can exist, while for * > 0, given by (11), both the trivial and the nontrivial equilibrium may exist for system (1). Substituting (9) into the seventh equation of system (1), From (7) and (13) we obtain either From expression (14), it should be noted that, for < thres , * > 0 whenever On the other hand, from fifth equation of system (1), we also get * = [ and * > 0 whenever From expressions (16) and (18), one can note, however, that, 2 is the maximum value of * , so it follows that 1 < 2 and * > 0 if and only if condition (18) holds. Furthermore, for * > 0 the system (1) reaches an endemic equilibrium point.
The positive endemic equilibrium of model (1) is obtained by solving for * from the quadratic (20) and substituting the results (positive values of * ) into the expressions that give the coordinates of equilibrium point.
First, it is straightforward to note that thres given by (23) and given by (22) are positive. Moreover, under our assumption Ψ < 0 and 1 is also positive. Clearly, the coefficient 2 of the quadratic equation (20) is always positive.
Also, note that the coefficient 0 is positive either if > 0 and < thres or < 0 and > thres . If < 0 then > thres and the positive endemic equilibrium does not exist (see (12)) for system (1). Later we will see that, for > thres , the only equilibrium biologically feasible and stable is the equilibrium given by 0 .
Therefore, 0 > 0 if and only if > 0 (i.e., < thres ) and < thres . In this case, if 1 > 0, the quadratic equation (20) does not have a positive solution. It follows then that system (1) does not have positive endemic equilibria whenever < thres . In this case, note that the only equilibrium biologically feasible and stable for system (1) is the trivial equilibrium given by 0 (see (19)).
Quite apart from this, the existence of the positive equilibrium points for the system (1) can be also summarized as follows.  (1) has the trivial equilibrium 0 and the endemic equilibrium 1 . For < ℎ , the unique equilibrium that exists is given by the trivial equilibrium 0 . Otherwise, for > ℎ , both the trivial equilibrium, 0 , and the endemic equilibrium, 1 , exist.
Having found the scenarios in which there exist the equilibria for the system (1), it is instructive to analyse whether or not these equilibria are stable under any of these scenarios. Furthermore, together with the threshold vaccination rate, ] , and the reproduction number, vac , we will see that each scenario can be used as a check for the existence, the uniqueness, and the stability of all equilibria. This is explored below for > .
To establish the stability of both trivial equilibrium, the Jacobian of the system (1) is computed and evaluated at both 0 and 0 . We will discuss the properties of both trivial equilibrium points making an elementary rowtransformation for the Jacobian matrix.
Evaluating the system's Jacobian at 0 , the local stability of 0 is straightforward determined by the six eigenvalues given by 1 It is easy to verify that both the traces of the matrices tr( 0 ) and tr( 1 ) are always negative. Moreover, the determinant of the matrix, det( 0 ), is always positive, but det( 1 ) is positive if and only if > thres . In other words, it means that the four eigenvalues of matrix 0 are either negative or have negative real part whenever > thres . Therefore, all the eigenvalues of the characteristic equation associated with the system (1) at 0 have negative real parts if and only if > thres and > . We state then the following result.
The stability of the disease-free equilibrium 0 is now examined by linearizing the system (1) around 0 . The characteristic equation of the Jacobian matrix of the system (1) at 0 is given by where with 0 , 0 , and 0 given by (25) and where 0 = 0 (see (27)) and It can be seen after some calculations that the polynomial (33) is equivalent to with vac < 1 given by (5); 0 is the threshold quantity or the basic reproductive number of the diseases defined by (23) and vac < 0 . It is worth remembering that 0 given by (25) exists whenever > and < thres , so 0 > 0 and vac > 0. By using the Routh-Hurwitz criteria for the polynomial (33), it follows that 2 > 0, 1 > 0, and 0 > 0 if and only if vac < 1 and 1 2 − 0 > 0. Therefore, the polynomial (33) has negative (or has negative real part) roots if vac < 1.
Hence, for > , all the eigenvalues of the characteristic equation associated with the system (1) at 0 are negative or have negative real parts if and only if < thres and vac < 1. Hence, the disease-free equilibrium 0 is locally asymptotically stable when vac < 1 and < thres . It is worth remembering that if < thres and < thres there are no positive solutions of the quadratic equation (20) and thus there is no endemic equilibrium of system (1). However, by (19), it follows that the disease-free equilibrium 0 is the only equilibrium point that exists when < thres . Otherwise, for > thres , there is a unique positive real solution of the quadratic equation (20) which indicates the possibility of a unique positive endemic equilibrium given by 1 . On the other hand, from (25) it can be also noted that the existence of the disease-free equilibrium 0 does not depend on thres , and thus 0 could also exist in this case. Therefore, the disease-free equilibrium 0 coexists with an endemic equilibrium 1 whenever > thres . It follows from above analyses that the disease-free equilibrium 0 is a unique equilibrium which is locally asymptotically stable whenever < thres and < thres . In such a scenario, disease elimination would depend upon the birth rate of humans, , and the mosquitoes natural mortality rate . Otherwise, if < thres and > thres , then 0 is locally asymptotically stable when vac < 1. The epidemiological implication of this is that the requirements of < thres and > thres are, although necessary, no longer sufficient for disease elimination. The disease elimination would depend upon the vaccination coverage ( ) and thus will eliminate the disease from the population if the usage of an imperfect vaccine results in making (and keeping) vac < 1.
For the case when < thres , for < thres the quadratic equation has no positive solution. Hence, system (1) has no positive solution (endemic equilibrium) for < thres . As the only existing equilibrium for < thres is the trivial one, 0 , then we can conjecture that, for < thres , 0 is the only possible equilibrium for system (1), being, therefore, stable. On the other hand, for > thres the quadratic equation has a positive solution and hence, system (1) shows an endemic equilibrium, 1 . As the existence of 0 , however, does not depend on thres (see (25)), we can state that 0 also exists for > thres . In addition, if vac < 1, then 0 is stable. Note that, in all cases, the mosquito's mortality rate is the critical parameter, in the absence of vaccination, determining the stability of the disease.
How about the case when vac > 1, 0 becomes unstable? In this case, the critical value of vaccination is given by (38); that is, if > thres , then 0 is stable. In contrast, if < thres , then 0 is unstable and 1 becomes stable. We establish then the following stability result for the system (1).

Lemma 3. For
> and < ℎ the diseasefree equilibrium of the model (1), 0 , exists and it is locally asymptotically stable if (a) V < 1 and > ℎ , The quantity vac given by (36) is called the vaccinated reproduction number, since it represents the expected average number of new infections produced by a single infective when introduced into a human community where a fraction of the susceptible population has been vaccinated. For a disease in which the susceptible population is vaccinated it has been demonstrated that vac , which is the basic reproduction number 0 [26] modified by vaccination, must be reduced below one in order to ensure that the disease dies out [27]. If there is no vaccination, then vac = 0 . Therefore, the aim of the vaccination must be to reduce vac below one and to provide prolonged protection against the infection. Now, expression (36) for vac can be written in terms of 0 , 0 , and 0 such as Setting vac = 1, and solving (37) for , the threshold vaccination rate is found to be thres = ( + ) ( 0 − 1) .
In contrast, if 0 > 1, thres is positive but we also have either vac > 1 or vac < 1. In this situation, the disease could be eliminated whenever > thres and 0 is globally asymptotically stable. Otherwise, when < thres , the use of an imperfect vaccine will fail to eliminate the disease from the community, 0 is unstable, the disease will persist in the community, and 1 becomes globally asymptotically stable.
Finally, concerning the analysis of the polynomial (20), Lemmas 2 and 3, we can establish the following conjecture for the system (1).

Conjecture 4. For
> the disease-free equilibrium 0 is globally asymptotically stable if > ℎ . If < ℎ , 0 becomes unstable and for < ℎ , 0 is globally asymptotically stable equilibrium point. If > ℎ , then 0 is globally asymptotically stable if > ℎ . Otherwise, if < ℎ , then 0 becomes unstable and the globally asymptotically equilibrium point stable is then given by 1 .
The stability of the equilibrium points can then be summarized in Table 3.

Numerical Analysis
The numerical analysis of the stability of the equilibria was done with the parameters of the model fixed at the baseline values indicated in Table 2. We explore the implications of variable vaccination coverage ( ), birth rate of humans ( ), and mosquitoes natural mortality rate ( ) which are chosen for simulations purposes only, so we can illustrate our theoretical results. For the baseline parameters values in Table 2, we have thres = 0.000035085, thres = 4.751131, and thres = 8.0117. Figures 2 and 3 show the graph of (20), with ( * ) plotted versus . Figure 2 shows the graph for < thres and < thres , with increasing values of . One can see two real negative roots of (20), which are not biologically viable for the system (1). Therefore, the unique equilibrium point that is biologically viable and locally asymptotically stable is given by 0 (see (19)). In this case, the disease can be eradicated from the population. Figure 3 shows the graph for < thres and > thres , with increasing values of . One can see one positive and one negative real root of (20). Hence, the system (1) has a unique endemic equilibrium 1 , and the disease persists at an endemic level. Parameter values used are as given in Table 2 (baseline values), except for . Figure 4 shows the graph of (37) for < thres and > thres , with increasing values of . Parameter values used are as given in Table 2 (baseline values), except for . For > thres , vac < 1; the disease is then eradicated from population ( 0 is locally asymptotically stable). For < thres , then vac > 1; the disease persists into the population ( 1 is the locally asymptotically stable).  Figure 5 shows the prevalence of infectious individuals as a function of . All parameters values used are given in Table 2 (baseline values), except for . For > thres , * = 0, and * = 0; thus, 0 is globally asymptotically stable. For < thres , then * ̸ = 0 and * ̸ = 0 and 1 is globally asymptotically stable. Figures 6 and 7 show the profiles of both infectious populations * and * . In Figure 6, for < thres and < thres there is no positive real solution of ( * ) and the DFE 0 is globally asymptotically stable. In Figure 7, for < thres and > thres the quadratic equation ( * ) has one positive real solution. For < thres , the endemic equilibrium 1 is, therefore, a unique equilibrium globally asymptotically stable. Parameter values used are as given in Table 2 (baseline values), except for , , and .

Sensitivity Analysis
In this section, we present the sensitivity analysis of the model to find out the degree to which the parameters influence the outputs of the model. Using the equation described in [15,28,29], we investigate only two of the most significant epidemiological concepts that affect the disease dynamics: the force of infection and the prevalence of infection.
The force of infection is defined as [28,29] = * * , where and are given in Table 2. Thus, substituting (7) into expression (40) gives Figure 5: Prevalence of (a) infectious individuals ( * ) and (b) infectious mosquitoes ( * ) as a function of . 0 is a unique equilibrium globally asymptotically stable for > thres and 1 is globally asymptotically stable for < thres .    Figure 6: For < thres and < thres , 0 is a unique equilibrium globally asymptotically stable for any vaccination coverage, . (a) Profile of population of both infectious humans ( * ) and mosquitoes ( * ). (b) * is plotted versus * . The system approaches * = 0 and * = 0, that is, the system approaches DFE, 0 .
Finally, applying (45) to estimate the sensitivity of force of infection and the prevalence to the vaccination effort ( vac ) we can show that for every 1% of variation in vac results in a variation of approximated 0.6% and 844% in and prev, respectively. Therefore, the prevalence is 1400 more sensitive to vaccination than the force of infection. All parameters values used are given in Table 2 (baseline values).

Summary and Conclusions
In this paper, we have formulated a model for yellow fever disease with both human and vector populations as variables. We found four threshold parameters that control the development of the disease and the infectious status of the human population in the presence of a preventive vaccine whose protection may wane over time.
Our analysis is based on the assumption that the growth rate of the human population is positive, that is, − > 0, which is the case of the yellow fever affected populations. We, therefore, can conclude the following: (a) if the mortality rate of the mosquitoes is greater than the threshold, > thres , then the disease is naturally (without intervention) eradicated from the population; (b) if, in contrast, the mortality rate of the mosquitoes is less than the threshold, < thres , then the disease is eradicated from the populations only when  Figure 7: For < thres and > thres , 1 is a unique equilibrium locally asymptotically stable for < thres . (a) Profile of populations of both infectious humans ( * ) and mosquitoes ( * ). (b) * is plotted versus * . The system approaches * ̸ = 0 and * ̸ = 0, that is, the system approaches endemic equilibrium, 1 . the growing rate of humans is less than a threshold, < thres . Otherwise, > thres ; then the disease is eradicated only if vac < 1; (c) in case vac > 1, then the disease will be eradicated from the human population if the vaccination rate is greater than a threshold, > thres . Otherwise, < thres ; then the disease will establish itself among humans, reaching a stable endemic equilibrium. This conclusion derives from a rearrangement of (37); that is, ) or in words, the vaccination effort should be such that 1 minus it should be greater than the proportion of susceptible at equilibrium; (d) the prevalence is 1400 more sensitive to vaccination than the force of infection.
The model and analyses presented in this paper are intended to serve as a framework for testing alternative vaccination schedules taking into account the vaccine and disease induced mortality rates. This will be the subject of a future publication