Modeling the Dynamics of Arboviral Diseases with Vaccination Perspective

In this paper, we propose a model of transmission of arboviruses, which takes into account a future vaccination strategy in human population. A qualitative analysis based on stability and bifurcation theory reveals that the phenomenon of backward bifurcation may occur; the stable disease-free equilibrium of the model coexists with a stable endemic equilibrium when the associated reproduction number, R0, is less than unity. We show that the backward bifurcation phenomenon is caused by the arbovirus induced mortality. Using the direct Lyapunov method, we prove the global stability of the trivial equilibrium. Through a global sensitivity analysis, we determine the relative importance of model parameters for disease transmission. Simulation results using a nonstandard qualitatively stable numerical scheme are provided to illustrate the impact of vaccination strategy in human communities. Keywords-Mathematical model; Arboviral disease; Vaccination; Stability; Backward bifurcation; Sensitivity analysis; Nonstandard numerical scheme.


I. INTRODUCTION
Arboviral diseases are affections transmitted by hematophagous arthropods.There are currently 534 viruses registered in the International Catalogue of Arboviruses and 25% of them have caused documented illness in humans [20], [49], [42].Examples of these kinds of diseases are dengue, yellow fever, Saint Louis fever, encephalitis, West Nile Fever and chikungunya.A wide range of arbovirus diseases are transmitted by mosquito bites and constitute a public health emergency of international concern.According to WHO, dengue, caused by any of four closely-related virus serotypes (DEN-1-4) of the genus Flavivirus, causes 50-100 million infections worldwide every year, and the majority of patients worldwide are children aged 9 to 16 years [72], [84], [86].The dynamics of arboviral diseases like dengue or chikungunya are influenced by many factors such as humans, the mosquito vector, the virus itself, as well as the environment which affects all the present mechanisms of control directly For all mentioned diseases, only yellow fever has a licensed vaccine.However, some works are underway for development of a vaccine for dengue [10], [11], [33], [50], [73], [85], Japanese encephalitis [73], and Chikungunya [53], [54], [55], [46].Moreover, the existence of different strains of dengue virus, for example, makes the developpement of an efficient vaccine a challenge for scientists.However, according to the French newspaper Le Figaro, the SANOFI laboratory hopes to market in the second half of 2015, the first vaccine against dengue fever, with an overall efficacy of 61% vaccine among young people aged 9 to 16 years and the rate of protection against severe dengue 95.5% [39].Therefore, it is important to assess the potential impact of such vaccines in human communities.
In this paper, we formulate a model, described by differential equations, in which we include two aspects: vaccination in the human population and the aquatic stage in the vectors population.We perform a qualitative analysis of the model, based on stability and bifurcation theory.In particular, we use an approach based on the center manifold theory [19], [31], [43] to evaluate the occurrence of a transcritical backward bifurcation and, as a consequence, the presence of multiple endemic equilibria when the basic reproduction number R 0 is less than unity.Under the point of view of disease control, the occurrence of backward bifurcation has relevant implications for disease control because the classical threshold condition R 0 < 1, is no longer sufficient to ensure the elimination of the disease from the population.
The global stability of the trivial equilibrium and the disease-free equilibrium (the equilibrium without disease in both populations), whenever the associated thresholds (the net reproductive number N and the basic reproduction number R 0 ) are less than unity, is derived through the use of Lyapunov stability theory and LaSalle's invariant set theorem, and the approach of Kamgang and Sallet [48], respectively.
Through global sensitivity analysis, we determine the relative importance of model parameters for disease transmission.The analysis of the model is completed by the construction of a nonstandard numerical scheme which is qualitatively stable.
The rest of this paper is organized as follows.In Section II, we develop the mathematical model and give the invariant set in which the model is defined.In Section III, we compute two thresholds: the net reproductive number N and the basic reproduction number R 0 .Depending of the values of these thresholds, we identify two disease-free equilibria: the trivial equilibrium which corresponds to the extinction of vectors, when N ≤ 1, and the disease-free equilibrium (DFE) when N > 1 and R 0 < 1.The results concerning the local and global stability of these two equilibria are also given.The section IV is devoted to the existence of endemic equilibria and bifurcation analysis.Vaccine impact is evaluated in Section V. Uncertainty and sensitivity analysis and the construction of a stable numerical scheme, are made in sections VI and VII respectively.A conclusion completes the paper.
In this section we describe the mathematical model that we shall study in this paper.The formulation is mostly inspired, with some exceptions, by the models proposed in [30], [40], [68], [80].We assume that the human and vector populations Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 are divided into compartments described by timedependent state variables.This said, the compartments in which the populations are divided are the following ones: -For humans, we consider susceptible (denoted by S h ), vaccinated (V h ), exposed (E h ), infectious (I h ) and resistant or immune (R h ).Humans susceptible population is recruited at a constant rates Λ h .Each human subpopulation comes out from the dynamics at natural mortality rates µ h .The human susceptible population decreased following infection, which can be acquired via effective contact with an exposed or infectious vector at a rate λ h (the incidence rate of infection for humans), given by where m denote the alternatively sources of blood [1], [80], a is the biting rate per susceptible vector, βhv denotes the probability of transmission of infection from an infectious vector (E v or I v ) to a susceptible human (S h or V h ).We obtain the expression of λ h as follows: the probability that a vector chooses a particular human or other source of blood to bite can be assumed as . Thus, a human receives in average a N v N h + m bites per unit of times.Then, the infection rate per susceptible human is given a βhv In expression (1), the modification parameter 0 < η v < 1 accounts for the assumed reduction in transmissibility of exposed mosquitoes relative to infectious mosquitoes.It is worth emphasizing that, unlike many of the published modelling studies on dengue transmission dynamics, we assume in this study that exposed vectors can transmit dengue disease to humans.This is in line with some studies (see, for instance [34], [40], [87], [90]).However, it is significant to note that, in the case of Chikungunya for example, the exposed vectors do not play any role in the infectious process, in this case η v = 0.The vaccinated population is generated by vaccination of susceptible humans at a rate ξ.Further, it is expected that any future dengue vaccine would be imperfect [40], [68].Thus, vaccinated individuals acquire infection at a rate (1 − )λ h where is the vaccine efficacy.Exposed humans develop clinical symptoms of disease, and move to the infectious class at rate γ h .Infectious humans may die as consequence of the infection, at a disease-induced death rate δ, or recover at a rate σ.It is assumed that individuals successfully acquire lifelong immunity against the virus.
-Vectors population is classified into four compartments: susceptible (S v ), exposed (E v ), infectious (I v ) and aquatic (A v ).The aquatic state includes the eggs, larvae, and pupae.The vector population does not have an immune class, since it is assumed that their infectious period ends with their death.The population of vectors consists essentially of females which reached adulthood.A susceptible vector is generated by the adulthood females at rate θ.The susceptible vector population decreased following infection, which can be acquired via effective contact with an exposed or infectious human at a rate λ v (the incidence rate of infection for vectors), given by (2) where βvh is the probability of transmission of infection from an infectious human (E h or I h ) to a susceptible vector (S v ), where the modification parameter 0 ≤ η h < 1 accounts for the relative infectiousness of exposed humans in relation to infectious humans.Here too, it is assumed that susceptible mosquitoes can acquire infection from exposed humans [23], [40].Exposed vectors move to the infectious class with the rate γ v .As in the case of the outbreak of Chikungunya on Réunion Island, it has been shown that lifespan of the infected mosquitoes is almost halved.This particular feature of infected mosquitoes therefore influences the dynamics of the disease [32], [30].Thus, following Dumont and coworkers [29], [30], we assume in this work that the lifespan of a vector depends on its epidemiological status.So the average lifespan for susceptible mosquitoes is 1/µ v days, the average lifespan of the exposed mosquitoes is 1/µ 1 days and the average adult lifespan for infected vector is 1/µ 2 .Thus, we have 1/µ 2 ≤ 1/µ 1 ≤ 1/µ v (with equality for other arboviral diseases).We do not consider vertical transmission in this work, so only susceptible humans are recruited, and vectors are assumed to be born susceptible.
We are now in position to write the model (the meaning of the state variables and parameters are summarized in Table I and Table II: In model (3) the upper dot denotes the time derivative.K denote the carrying capacity of breeding sites.The parameter K is highly dependent on some factors such as (the location, temperature, season).In this work, we follow Dumont and Chiroleu [30], and consider, in our sensitive analysis, that K depend of the location, thus K = χN h , where χ is a positive integer which represent the location and N h the human population size.For example, in the year 2005 at Saint-Denis and Saint-Pierre in Réunion island, χ = 2) [30].µ b represent the number of eggs at each deposit per capita and µ A is the natural mortality of larvae.We set Let N h the total human population and N v the total of adult vectors, i.e.
and the other entries of A(X) are equal to zero; and F = (Λ h , 0, 0, 0, 0, 0, 0, 0, 0) T .Note that A(X) is a Metzler matrix, i.e. a matrix such that off diagonal terms are non negative [8], [47], for all X ∈ R 9 + .Thus, using the fact that F ≥ 0, system (4) is positively invariant in R 9  + , which means that any trajectory of the system starting from an initial state in the positive orthant R 9 + , remains forever in R 9 + .The right-hand side is Lipschitz continuous: there exists an unique maximal solution.
On the other hand, from the first four equations of model system (3), it follows that From the last three equations of model system (3), it also follows that So that Therefore, all feasible solutions of model system (3) enter the region:

III. THE DISEASE-FREE EQUILIBRIA AND STABILITY ANALYSIS
Now let N the net reproductive number [25] given by We prove the following result Proposition 3.1: a) If N ≤ 1, then, system (3) has only one trivial disease-free equilibrium Proof: See Appendix A. In Proposition 3.1, we have shown that system (3) have at least two equilibria depending of the value of treshold N and the basic reproduction number R 0 .Precisely, we have proved that when N ≤ 1, model sytem (3) admits only one equilibrium called trivial equilibrium (T E := P 0 ).When N > 1, model sytem (3) admits additionally the disease-free equilibrium (DF E := P 1 ).We prove, in the following, that the trivial equilibrium is locally and globally asymptotically stable whenever N ≤ 1.Then, we prove that the trivial equilibrium is unstable when N > 1, and the disease-free equilibrium is locally asymptotically stable whenever R 0 < 1.Using Kamgang and Sallet approach [48], a necessary condition for the global stabilty of the disease-free equilibrium is also given.

A. Local stability analysis
The local stability of the trivial equilibrium and the disease-free equilibrium is given in the following result: Proposition 3.2: a) If N ≤ 1, then the trivial equilibrium TE is locally asymptotically stable.b) If N > 1, then the trivial equilibrium is unstable and the Disease Free Equilibrium P 1 is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1, where R 0 is the basic reproduction number [26], [82], given by Proof: See appendix B. Proof: See Appendix C. 2) Global asymptotic stability of the diseasefree equilibrium : Following [30], we prove that the disease-free equilibrium DF E := P 1 is globally asymptotically stable under a certain threshold condition.To this aim, we use a result obtained by Kamgang and Sallet [48], which is an extension of some results given in [82].Using the property of DFE, it is possible to rewrite (3) in the following manner ) where X S is the vector representing the state of different compartments of non transmitting individuals (S h , V h , R h , A v , S v ) and the vector X I represents the state of compartments of different transmitting individuals (E h , I h , E v , I v ).Here, we have A direct computation shows that the eigenvalues of A 1 (X) are real and negative.Thus the system ẊS = A 1 (X)(X S − X DF E ) is globally asymptotically stable at X DF E .Note also that A 2 (X) is a Metzler matrix, i.e. a matrix such that off diagonal terms are non negative [8], [47].
Following D, we now consider the bounded set G: Let us recall the following theorem [48] Theorem 3.1: (1) G is positively invariant relative to (11).
(2) The system ẊS = A 1 (X)(X S − X DF E ) is Globally asymptotically stable at X BRDF E .(3) For any x ∈ G, the matrix A 2 (X) is Metzler irreducible.(4) There exists a matrix Ā2 , which is an upper bound of the set Then the DFE is GAS in G. (See [48] for a proof).
Let us now verify the assumptions of the previous theorem: it is obvious that conditions (1-3) of the theorem are satisfied.An upper bound of the set of matrices M, which is the matrix Ā2 is given , To check condition (5) in theorem 3.1, we will use the following useful lemma [48] [48].In our case, we have Clearly, A is a stable Metzler matrix.Then, after some computations, we obtain D − CA −1 B is a stable Metzler matrix if and only if where Thus we claim the following result Theorem 3.2: This shows that with the establishment of an effective treatment, it is possible to have R 0 and R G less than 1.Theorem (3.2) means that for R 0 < R G < 1, the DFE is the unique equilibrium (no co-existence with an endemic equilibrium then it is possible to have co-existence with endemic equilibrium.To confirm whether or not the backward bifurcation phenomenon occurs in this case, one could use the approach developed in [19], [31], [82], which is based on the general center manifold theory [43].

A. Existence of endemic equilibria
We now turn to study the existence of an endemic equilibrium of model system (3).Let R 0 the basic reproduction number [26], [82] given by Eq. (10).
we claim the following Proposition 4.1:
Proof: See appendix D.
The backward bifurcation phenomenon, in epidemiological systems, indicate the possibility of existence of at least one endemic equilibrium when R 0 is less than unity.Thus, the classical requirement of R 0 < 1 is, although necessary, no longer sufficient for disease elimination [6], [14], [40], [75].In some epidemiological models, it has been shown that the phenomenon of backward bifurcation is caused by factors such as nonlinear incidence (the infection force), disease-induced death or imperfect vaccine [15], [16], [31], [40], [70], [75].
It is important to note that case (i) of Proposition 4.1 suggests the possibility of a pithcfork (Forward) bifurcation when R 0 = 1.Also, case (iv) of Proposition 4.1 suggests that the principal cause of occurence of backward bifurcation phenomenon is the disease-induced death in both humans and vectors.
In the following remark, we prove that, in absence of disease-induced death in both populations, the disease-free equilibrium is the unique equilibrium whenever N > 1 and R 0 | δ=0,µv=µ1=µ2 < 1.Using the direct Lyapunov method, we prove the global asymptotic stability of the disease-free equilibrium whenever In the absence of disease-induced death, i.e, δ = 0 and with Equation ( 15) have only one positive solution whenever R 1 > 1.If R 1 ≤ 1, then coefficients B 00 , B 01 , B 02 are always positive, and the diseasefree equilibrium is the unique equilibrium.From this we conclude that the disease-induced mortality is the possible cause for the occurrence of multiple endemic equilibria below the classical threshold R 1 = 1.
The global stability of the disease-free equilibrium may be achieved by Lyapunov method.To this aim, let us consider the following Lyapunov function [37], [40] are positive weights given by g 1 = 1; .
Along the solutions of (3) we have: After replacing the constants g i , i = 1, . . ., 4 by their value, and using the fact that in the first, second, fifth, sixth and seventh equations of Eq. ( 3) with It follows from the LaSalle's invariance principle [45] that every solution of (3) (when R 1 ≤ 1), with initial conditions in D 1 converges to P 1 , as t → ∞.Hence, the DFE (P 1 ), of model ( 3), is GAS in

B. Bifurcation analysis
Previous Analysis state that multiple endemic equilibria may occur althougt R 0 < 1.In order to better investigate the variation of model's prediction as R 0 varied, we perform a bifurcation analysis at the criticality, i. e. at the Diseasefree Equilibrium DF E := P 1 and R 0 = 1.On one hand, this will provide a rigorous proof that the occurrence of multiple endemic equilibria comes from a backward bifurcation.On the other hand, we will get also information on the stability of equilibria near the criticality.In particular, on the stability of the endemic equilibrium points emerging from the criticality.We study the center manifold near the criticality by using the approach developed in [19], [31], [82], which is based on the general centre manifold theory [43].In summary, this approach establishes that the normal form representing the dynamics of the system on the center manifold is given by u = a u 2 + b u, where, u represent a real function of real variable, and (17) Note that the symbol denotes a bifurcation parameter to be chosen, f i 's denotes the right hand side of system (3), x denotes the state vector, x 0 the Disease-free Equilibrium P 1 ; v and w denote the left and right eigenvectors, respectively, corresponding to the null eigenvalue of the Jacobian matrix of system (3) evaluated at the criticality.
In order to apply this approach, let us choose β hv as bifurcation parameter.From R 0 = 1 we get the critical value .
Note also that in f k (0, 0), the first zero corresponds to the disease-free equilibrium, P 1 , for the system (3).Since The Jacobian matrix of the system (4) evaluated at the disease-free equilibrium P 1 with β hv = β * hv is given by , 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 where we have set The eigenvalues of the Jacobian matrix J(P 1 ) are λ 1 = −µ h , λ 2 = −k 1 , and the other eigenvalues are the eigenvalue of the following matrix The characteristic polynomial of J is given by f (λ) = λ 6 +a 5 λ 5 +a 4 λ 4 +a 3 λ 3 +a 2 λ 2 +a 1 λ+a 0 (18) with The others coefficients a 5 , a 4 , a 3 , a 2 , and a 1 are obtained after computations on Maxima software [58], [89].Since at the criticality, we have R 0 = 1, then a 0 = 0, thus equation (18 Then, the Jacobian J(P 1 ) of the linearized system has a simple zero eigenvalue and therefore P 1 is a nonhyperbolic equilibrium for R 0 = 1.In order to get the coefficients ( 16) and ( 17), we need to calculate the right and the left eigenvectors corresponding to the zero eigenvalue.The right eigenvector of J(P 1 ) is given by w = (w 1 , w 2 , w 3 , w 4 , w 5 , w 6 , w 7 , w 8 , w 9 ) T where a) Computation of a : Using the non-zero components of v and the associated non-zero partial derivatives of f (at the DFE P 1 ), for system (3), we obtain We finally obtain (See the details in appendix E) Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241where and Since b > 0 according to the sign of w i , v i , for i ∈ {1 . . ., 9}, we conclude that the backward bifurcation phenomenon may occurs if a > 0.
We can summarize the results obtained above in the following theorem: Theorem 4.1: If a > 0, then model (3) exhibits backward bifurcation at R 0 = 1.If the reversed inequality holds, then the bifurcation at R 0 = 1 is forward.This is illustrated by simulating the model with different set of parameter values (it should be stated that these parameters are chosen for illustrative purpose only, and may not necessarily be realistic epidemiologically): -Using the parameters values in Table II -Using the parameters values in Table II -In the absence to disease induced death (δ = 0) and choosing βhv = 4.0188 and K = 1000 such that R 0 = 1, equation ( 44) admit only one solution λ * h = 0 which corresponds to the disease-free equilibrium.In this case, the backward bifurcation phenomenon does not occurs.
To conclude, depending to the values of parameters of model ( 3), the phenomenon of backward bifurcation may occurs when the classical basic reproduction number R 0 is less than unity.

IMPACT
Since a future dengue vaccine, for example, is expected to be imperfect, it is instructive to determine whether or not its widespread use in a community will be benefic (or not) [10], [40], [68].Now, consider the following model (model 3 without vaccination).
with λ h and λ v defined at (1) and (2), respectively.Following procedure in [26], [82], the corresponding basic reproduction number of model (19), R s , is given by So we deduce that From Eq. ( 21), it follows that, in the absence of vaccination (ξ = 0) or when the vaccine efficacy is very low ( → 0), we have R vac = R s .However, when humans vaccination comes to play, the basic reproductive number is reduced by a factor of (πξ + µ h ) (µ h + ξ) < 1.Since increasing vaccination efforts results in decreasing the magnitude of arboviruses infection, humans vaccination can contribute to control the spread of arboviral diseases.
In the following, we use the set of parameters values given in Table III, which refer to Dengue and Chikungunya.Figs.3-5 show several simulations, by varying the vaccine efficacy and the percentage of population that is vaccinated.Figure 3 shows simulations with different proportions of succeptible human vaccinated, using an imperfect vaccine, with a level of efficacy of 60%.Both total number of infected humans and infected vectors reache a peack after 25 days approximatively.However, when = 60%, the variation of vaccine coverage parameter have not a great impact in the number of infected humans and vectors.Figure 4 illustrates the effect of vaccine efficacy in the reduction of the total number of infected humans and vectors.It is clear that the effectiveness of the vaccine has a great impact when ≥ 90%.Thus, it is suitable to add to vaccination (when < 90%) another control, such as, treatment of infected individuals, personal protection, and vector control strategies to stop the spread of arboviral diseases.
Figure 5 shows the representation of the basic reproduction number R 0 plotted as function of the vaccine efficacy parameter and the proportion [68], [60], [61] of susceptible population vaccinated ξ.The use of a vaccine with level of efficacy greather than 90% approximatively, dramaticaly decrease the basic reproduction number, when the proportion of susceptible humans vaccinated are greather than 85%.We observe the same result at Figure 6.Thus, the use of a vaccine with a high level of efficacy and a wide vaccine coverage has an impact on stopping the spread of the disease.However, if the vaccine efficacy is not high, it is important to add another control strategies.Our sensitive analysis in later section will further support this result.

VI. SENSITIVITY ANALYSIS
To determine the best way to fight against arboviruses, it is necessary to know the relative importance of the various factors responsible for their transmission in both the human population than in the vector population, as well as effective means to fight these diseases.The transmission of the disease is directly related to R 0 , and the   prevalence of the disease is directly related to the infected states, especially for sizes of E h (t), I h (t), E v (t) and I v (t).These variables are relevant to  III. the individuals (humans and vectors) who have some life stage of arboviruses in their bodies.The number of infectious humans, I h , is especially important because it represents the people who may be clinically ill, and is directly related to the total number of arboviral deaths [22].We will perform a global sensitivity analysis.

A. Mean values of parameters and initial values of variables
Since we focus in this article, to a general model of arboviral diseases, we will, in this sec- tion, use the parameters values of two particular arboviruses, Dengue and Chikungunya.It is important to note that these two diseases are transmitted by the same mosquito: Aedes albopictus.However, dengue is also transmited by Aedes aegypti [30], [35], [36], [38], [40], [61], [68], [90].The mean values of parameters are listed in Table III, the range values of parameters are in Table IV and the initial conditions are given in Table V.

B. Uncertainty and sensitivity analysis
1) Sensitivity analysis of R 0 : We study the impact of each parameter of the model on the value of the basic reproduction number R 0 .Following the approach of Wu and colleagues [88], we perform the analysis by calculating the Partial Rank Correlation Coefficients (PRCC) between each parameter of our model and the basic reproduction number, R 0 .Table III  the mean value for each parameter.It is important to notice that, variations of these parameters in our deterministic model lead to uncertainty to model predictions since the basic reproductive number varies with parameters.Due to the absence of data on the distribution function, a uniform distribution is chosen for all parameters.The sets of input parameter values sampled using the Latin Hypercube Sampling (LHS) method were used to run 1,000 simulations.
With these 1,000 runs of Latin Hypercube Sampling, the derived sampling distribution of R 0 is shown in Figure 7. From this sampling we get that the mean of R 0 is 1.9304 and the standard deviation is 1.6128.Hence, statistically we are very confidential that model ( 3) is in an endemic state since R 0 > 1.From the previous sampling we compute the Partial Rank Correlation Coefficients between R 0 and each parameter of model ( 3).The result is displayed in Table VI.According to Boloye Gomero [13], the parameters with large PRCC values (> 0.5 or < −0.5) as well as corresponding small p-values (< 0.05) are most influential in model (3).
Table VI show that the parameter have the highest influence on the reproduction number R 0 .Although is the vaccine efficacy.This suggests that the development of a vaccine with high level of efficacy may potentially be an effective strategy to reduce R 0 .The other parameters with an important effect are θ, a, Λ h and µ 2 .The parameters which do not have almost any effect on R 0 are ξ, δ, µ A , m and µ b .In particular, the least sensitive parameters is µ b , the number of eggs at each deposit per capita.
2) Sensitivity analysis of Infected states of model (3): With 1,000 runs of Latin hypercube sampling, we compute the PRCC between infected states E h (t), I h (t), E v (t), and I v (t) and each parameter of model (3).The result is displayed in Tables VII-X.As in Table VI, the parameters with large PRCC values (> 0.5 or < −0.5) as well as corresponding small p-values (< 0.05) are most influential in model (3).
From Tables VII-X, we can observe the following facts: -For the value of E h , the parameters with more influence are θ, K, a, , Λ h and µ 2 .The parameters which do not have almost any effect on the variation of E h are µ h , δ, µ A , m and µ b .In particular, the least sensitive parameters is µ b , the number of eggs at each deposit per capita; -For the value of I h , the parameters with more influence are Λ h and γ h .The least sensitive pa-   rameters is µ b , the number of eggs at each deposit per capita; -For the value of E v , the parameters with more influence are the maturation rate from larvae to adult θ, and the capacity of breeding sites K.The other parameter is the natural mortality rate of vector µ v .The least sensitive parameters is m, the number of alternative source of blood; -For the value of I v , the parameters with more influence are θ, K, γ v , a and βvh .The least sensitive parameters is βhv , the probability of transmission of infection from an infectious vector to a susceptible human.
Although the model is sensitive to the variations of the vaccine efficacy parameter , there are other parameters (such as θ, a, K, µ v , µ 2 ) which have a considerable impact on the value of the basic reproduction number R 0 and the number of infected individuals.Thus, it is important to take into account other control strategies in the fight against arboviral diseases.

VII. NUMERICAL SIMULATION
In order to illustrate some of the results obtained in the previous sections, we provide here some numerical simulations.We use the nonstandard scheme given by (22).It is important to note that standard numerical methods may fail to preserve the dynamics of continuous models [4], [59], [81].Generally, compartmental models are solved using standard numerical methods, for example, Euler or Runge Kutta methods included in software package such as Scilab [76] or Matlab [57].These methods can sometimes lead to spurious behaviours which are not in adequacy with the continuous system properties that they aim to approximate.For example, they may lead to negative solutions, exhibit numerical instabilities, or even converge to the wrong equilibrium for certain values of the time discretization or the model parameters (see [3], [4], [5], [81] for further investigations).

A. A nonstandard scheme for the model (3)
Following [30], system (3) is discretized as follows: , and , which implies that the scheme is consistant with formulation (11).Rearranging (22), we obtain the foollowing new expression with ) .
where I 4 and I 5 are the identity matrix of order 4 and 5 respectively.Thus, we claim the following result: Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Lemma 7.1: Our non-standard numerical scheme (22) is positively stable, i.e. for X k ≥ 0 we obtain X k+1 ≥ 0, where Proof: We suppose X k ≥ 0. A k is a positive diagonal matrix and thus, A −1 k ≥ 0. B 12 is a positive matrix and we also have −A 1 (X k )X DF E ≥ 0. To show that B k is a positive matrix, it suffices to choose φ(h) such that where A 1 and A 2 are lower bounds for the sets {X ∈ D|A 1 (X)} and {X ∈ D|A 2 (X)} respectively.Following [30], to have B k ≥ 0, it suffices to consider the following time-step function with Concerning the equilibria of our numerical scheme, we have the following result Lemma 7.2: Our non-standard numerical scheme (22) and the continuous model (3) have the same equilibria.
Proof: See appendix F. The stability of the trivial equilibrium is given by the foollowing result Lemma 7.3: If φ(h) has choosen as equation (25), then the tivial equilibrium T E := P 0 is locally asymptotically stable for our numerical scheme (22) whenever N ≤ 1.
Proof: See appendix G. Now, we also have the following result concerning the stability of the disease-free equilibrium: Lemma 7.4: If φ(h) has choosen as equation ( 25) and R 0 < 1, then the disease-free equilibrium DF E := P 1 is locally asymptotically stable for our numerical scheme (22) .
Proof: The proof of Lemma 7.4 follows the proof of Proposition 3.4 in [30].See also [5] for a proof in a more general setting.

B. Simulation Results
We now provide some numerical simulations to illustrate the theoretical results (local stability, global stability and backward bifurcation).We use parameters values given in Table III with ξ = 0.475, = 0.60, K = 1000 and initial conditions given in Table V.
Figure 8 illustrates the asymptotic stability of the trivial equilbrium whenever the treshold N is less than unity.In Figure 9, when N > 1 the trivial equilibrium is unstable and the diseasefree equilibrium is stable (first panel).The phenomenon of backward bifurcation occurs in the second panel of figure 9, where the stable diseasefree equilibrium of the model co-exists with a stable endemic equilibrium when the associated reproduction number, R 0 , is less than unity.Figures 10-11 show the existence of at least one endemic equilibrium whenever R 0 is equal or greather than unity.
It is important to mentione that the simulation results discussed in this work are subject to the uncertainties (See section VI) in the estimates of the parameter values (tabulated in Table III) used in the simulations.The effect of such uncertainties on the results obtained can be assessed using a sampling technique, such as Latin Hypercube Sampling.

VIII. CONCLUSIONS
In this paper, we formulated a compartmental model which takes into account a future vaccination strategy in human population, the aquatic development stage of vector and the alternative sources of blood.
The analysis has been performed by means of stability, bifurcation and sensitivity analysis.We have obtained that the disease-induced mortality may be the main cause for the occurrence of the backward bifurcation (see remark 4.1).This means that relatively high values of diseaseinduced mortality rate may induce stable endemic states also when the basic reproduction number R 0 is below the classical threshold R 0 = 1.The stability analysis reveals that for N ≤ 1, the trivial equilibrium is globally asymptotically stable.When N > 1 and R 0 < 1, the diseasefree equilibrium is locally asymptotically stable.In  the absence of disease-induced death, the diseasefree equilibrium is also globally asymptotically stable.The reduced version of the model (3) (in the absence of disease-induced mortality in both human and vector populations) have a unique endemic equilibrium point whenever its associated reproduction number R 1 exceeds unity.
Taking as cases study the dengue and chikungunya transmission, we used parameter values from the literature to estimate statistically the basic reproduction number, R 0 , and to perform a global sensitivity analysis on the basic reproduction number and infected states (E h , I h , E v , I v ).Using Latin Hypercube Sampling, we obtain that the mean of R 0 is 1.9304.Hence, statistically we are very confident that our model ( 3) is in an endemic state.The global sensitivity analysis reveals that, apart from the parameters related to vaccination, particularly vaccine efficacy, other parameters ( such as parameters related to vector population) also have a great impact on the basic reproduction number (R 0 ) and on the number of infected humans and vectors (E Numerical simulations of the model (3), using a nonstandard qualitatively stable scheme, show that the use of a vaccine with high level of efficacy has a proponderant role in the reduction of the disease spread.However, since the efficacy of the proposed vaccine for dengue, for example, has been around 60%, it is suitable to combine vaccination with other mechanisms of control.Also, to be the best control strategy, the vaccination process must verify the following conditions: (a) The vaccine must be approved by the relevant agencies (such as WHO, CDC), before its use.(b) The vaccine efficacy should be high, as well as vaccine coverage.(c) The price of the vaccine must be low for countries affected by the disease.
There are already governments, affected by the diseases, willing to use the vaccine before it is approved, which can have unpredictable consequences, so condition (a) does not hold.Moreover, according to previous analysis and french laboratory SANOFI, the condition (b) does not hold.Thus it is important to know what happens when we combine vaccination with other mechanisms of control already studied in the literature, such as personal protection, chemical interventions and education campaigns [30], [40], [60], [61], [63], [64], [67], [68], [69].This is the perspective of our work.

ACKNOWLEDGMENT
Hamadjam ABBOUBAKAR and Léontine Nkague NKAMBA thank the Direction of UIT of Ngaoundéré and ENS of Yaoundé I, respectively, for their financial help granted under research missions in the year 2014.The authors are very grateful to two anonymous referees, for valuable remarks and comments, which significantly contributed to the quality of the paper.

APPENDIX
Appendix A: PROOF OF PROPOSITION 3.1 To find the equilibrium points of our system, we will solve the following system represents any arbitrary endemic equilibrium of (3).Further, let Solving the last three equations in (26) at steady state gives where Substituting (29) in the sixth equation of (26) gives The trivial solution of ( 30) is A * v = 0. Substituting this solution in (29) gives Then we obtain the trivial equilibrium , 0, 0, 0, 0, 0, 0, 0 .Now we suppose that A * v = 0.The possible solution(s) of ( 30) is the solution(s) of the following equation The direct resolution of equation (31) gives where Let us first compute the equilibrium without Disease, i.e. (32), we obtain Thus, the existence of vector is regulated by the threshold N .If N ≤ 1, the system (3) correspond to human population of free vectors and the trivial equilibrium in this case is P 0 .Now we suppose that N > 1.From ( 28) and ( 29) (with λ * v = λ * v = 0), we obtain the non trivial equilibrium or the disease-free equilibrium P 1 = S 0 h , V 0 h , 0, 0, 0, A 0 v , S 0 v , 0, 0 , where Appendix B: PROOF OF PROPOSITION 3.2 We consider the Jacobian matrix associated to model ( 3) at the equilibrium T E. we have were S 0 = S 0 h + πV 0 h .The eigenvalues of J (T E) are given by λ The characteristic polynomial of J is given by where Thus, it is clear that all coefficients are always positive since N < 1.Now we just have to verify that the Routh-Hurwitz criterion holds for polynomial P(λ).To this aim, setting The Routh-Hurwitz criterion of stability of the trivial equilibrium T E is given by We have Thus we conclude that the trivial equilibrium is locally asymptotically stable.Now we assume that N > 1.Following the procedure and the notation in [82], we may obtain the basic reproduction number R 0 as the dominant eigenvalue of the next-generation matrix [26], [82].Observe that model (3) has four infected populations, namely E h , I h , E v , I v .It follows that the matrices F and V defined in [82], which take into account the new infection terms and remaining transfer terms, respectively, are given by and the dominant eigenvalue of the nextgeneration matrix F V −1 is given by Eq. (10).
Appendix D: PROOF OF PROPOSITION 4.1.
We compute now the endemic equilibrium, i.e. we are looking for an equilibrium such that λ * h = 0 and λ * v = 0. We assume that N > 1.From the sixth equation of ( 26), at equilibrium, we have From the last third equations of ( 26), at equilibrium, we have we will observe the following two cases.a) Absence of disease-induced death in vector: The absence of disease-induced death in vector is traduce by the relation µ v = µ 1 = µ 2 , then equation (38) becomes Equalling Eqs. ( 37) and ( 39) gives like before Substituting A * v by A 0 v in equation ( 29) gives Replacing (41) in the expression of λ * h gives where Replacing (28) in the expression of λ * v gives where k 11 = k 3 η h + γ h .Substituting (43) in (42) gives the following equation where We consider λ * h = 0, otherwise we recover DFE.The positive endemic equilibria are obtained by solving Eq. ( 44) for λ * h .The coefficient B 4 is always positive and coefficient B 0 is negative (resp.positive) whenever R 0 > 1 (resp.R 0 < 1).The number of possible nonnegative real roots of the polynomial of Eq. ( 44) depends on the signs of B 3 , B 2 and B 1 .The various possibilities for the roots of f (λ * h ) are tabulated in Table XI and XII.From tables XI and XII , we deduce the following result which gives various possibilities of nonnegative solutions of (44).
Lemma A.2: Assume that N > 1.Then, the arboviral-disease model (3) 1. could have a unique endemic equilibrium whenever R 0 > 1. 2. could have more than one endemic equilibrium whenever R 0 > 1. 3. haven't endemic equilibrium whenever R 0 < 1. 4. could have one or more than one endemic equilibrium whenever R 0 < 1. Case 4 of Lemma A.2 suggests that co-existence of the disease-free equilibrium and endemic equilibrium for the arboviral-disease model (3) is possible, and hence the potential occurrence of a backward bifurcation phenomenon when R 0 < 1.Also, case 2 of Lemma A.2 suggests the possibility of a pithcfork (Forward) bifurcation when R 0 = 1.and

Appendix
T is an equilibrium of the discrete system (49), then we have after few simplifications ) which is equivalent to A 1 (X * )(X * S − X DF E ) + A 12 (X * )X * I = 0 A 2 (X * )X * I = 0 (51) where A 1 , A 12 and A 2 are given at Equation (11).

Fig. 2 .
Fig. 2. Time profile of infectious humans using different initial conditions showing that the equilibrium λ * 33h = 8.5310 is stable even if R0 = 1 .

Biomath 4 (Fig. 5 .
Fig. 5.The basic reproduction number R0 plotted as function of the vaccine efficacy parameter and the proportion of susceptible population vaccinated ξ.The set of parameter values is given in TableIII.

Fig. 6 .
Fig.6.Time profile of total number of infected human and vector without vaccination and with vaccination.

FrequenciesFig. 7 .
Fig. 7. Sampling distribution of R0 from 1,000 runs of Latin hypercube sampling.The mean of R0 is 1.9304 and the standard deviation is 1.6128.

TABLE VI PRCC
BETWEEN R0 AND EACH PARAMETER.

TABLE XI TOTAL
NUMBER OF POSSIBLE REAL ROOTS OF (44) WHEN R0 > 1.