www.biomathforum.org/biomath/index.php/biomath ORIGINAL ARTICLE Modeling the Dynamics of Arboviral Diseases with Vaccination Perspective Hamadjam Abboubakar∗, Jean C. Kamgang†, Léontine N. Nkamba‡, Daniel Tieudjo† and Lucas Emini¶ ∗Department of Computer Science, UIT–University of Ngaoundéré, Cameroon abboubakarhamadjam@yahoo.fr †Department of Mathematics and Computer Science, ENSAI–University of Ngaoundéré, Cameroon jckamgang@yahoo.fr, tieudjo@yahoo.com ‡Department of Mathematics, ENS–University of Yaoundé I, lnkague@gmail.com ¶Department of Mathematics, Polytechnic–St. Jerome Catholic University, Cameroon lemini@univ-catho-sjd.com Received: 31 December 2014, accepted: 24 July 2015, published: 27 August 2015 Abstract—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 impor- tance 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 com- munities. Keywords-Mathematical model; Arboviral dis- ease; 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 Cat- alogue 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, en- cephalitis, West Nile Fever and chikungunya. A wide range of arbovirus diseases are transmit- ted by mosquito bites and constitute a public health emergency of international concern. Ac- cording 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 af- fects all the present mechanisms of control directly Citation: Hamadjam Abboubakar, Jean C. Kamgang, Léontine N. Nkamba, Daniel Tieudjo, Lucas Emini, Modeling the Dynamics of Arboviral Diseases with Vaccination Perspective, Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 1 of 30 http://www.biomathforum.org/biomath/index.php/biomath http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... or indirectly. 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. As part of the necessary multi–disciplinary re- search approach, mathematical models have been extensively used to provide a framework for un- derstanding arboviral diseases transmission and control strategies of the infection spread in the host population. In the literature, considerable works can be found on the mathematical modeling of specific arboviral diseases, like West Nile Fever, yellow fever, dengue and chikungunya, see e.g. [2], [17], [24], [30], [35], [36], [38], [40], [56], [60], [61], [64], [68], [79]. Although these models highlight the means to fight against these ar- bovirus, few papers deal with models that consider vaccination [40], [68], [79]. 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 R0 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 R0 < 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 R0) 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 deter- mine 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 R0. 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 R0 < 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. Vac- cine 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. II. MODEL FORMULATION, INVARIANT REGION. In this section we describe the mathematical model that we shall study in this paper. The for- mulation 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 Page 2 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... are divided into compartments described by time– dependent state variables. This said, the compart- ments in which the populations are divided are the following ones: –For humans, we consider susceptible (denoted by Sh), vaccinated (Vh), exposed (Eh), infectious (Ih) and resistant or immune (Rh). Humans sus- ceptible 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 con- tact with an exposed or infectious vector at a rate λh (the incidence rate of infection for humans), given by λh = aβ̃hv Nv Nh + m (ηvEv + Iv) Nv = βhv (ηvEv + Iv) Nh + m , (1) 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 (Ev or Iv) to a susceptible human (Sh or Vh). 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 1 Nh + m . Thus, a human receives in average a Nv Nh + m bites per unit of times. Then, the infection rate per suscepti- ble human is given aβ̃hv Nv Nh + m (ηvEv + Iv) Nv . 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 stud- ies 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 vac- cination 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 com- partments: susceptible (Sv), exposed (Ev), infec- tious (Iv) and aquatic (Av). 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 popula- tion 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 λv = aβ̃vh (ηhEh + Ih) Nh Nh Nh + m = βvh (ηhEh + Ih) Nh + m (2) where β̃vh is the probability of transmission of infection from an infectious human (Eh or Ih) to a susceptible vector (Sv), 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 par- ticular 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 Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 3 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... TABLE I THE STATE VARIABLES OF MODEL (3). Humans Vectors Sh: Susceptible Av Aquatic Vh: Vaccined Sv: Susceptible Eh: Infected Ev: Exposed Ih: Infectious Iv: Infectious Rh: Resistant (immune) 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:  Ṡh = Λh −λhSh − ξSh −µhSh V̇h = ξSh − (1 − �)λhVh −µhVh Ėh = λh [Sh + (1 − �)Vh] − (µh + γh)Eh İh = γhEh − (µh + δ + σ)Ih Ṙh = σIh −µhRh Ȧv = µb ( 1 − Av K ) (Sv + Ev + Iv) − (θ + µA)Av Ṡv = θAv −λvSv −µvSv Ėv = λvSv − (µ1 + γv)Ev İv = γvEv −µ2Iv (3) In model (3) the upper dot denotes the time deriva- tive. 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 Chi- roleu [30], and consider, in our sensitive analysis, that K depend of the location, thus K = χNh, where χ is a positive integer which represent the location and Nh 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. TABLE II DESCRIPTION OF PARAMETERS OF MODEL (3). Parameter Description Λh Recruitment rate of humans ξ Vaccine coverage � The vaccine efficacy ηh, ηv Modification parameters β̃hv Probability of transmission of infection from an infectious vector to a susceptible human β̃vh Probability of transmission of infection from an infectious humans to a susceptible vector γh Progression rate from Eh to Ih γv Progression rate from Ev to Iv µh Natural mortality rate in humans µv Natural mortality rate of susceptible vectors µA Natural mortality of larvae µ−11 Average lifespan of exposed mosquitoes µ−12 Average lifespan of infected mosquitoes θ Maturation rate from larvae to adult δ Disease–induced death rate in humans σ Recovery rate for humans a Average number of bites m Number of alternative source of blood K Capacity of breeding sites µb Number of eggs at each deposit per capita We set π = 1 − �, k1 = µh + ξ, k2 = µh + γh, k3 = µh + δ + σ, k4 = µ1 + γv, k6 = µA + θ. Let Nh the total human population and Nv the total of adult vectors, i.e. Nh = Sh + Vh + Eh + Ih + Rh and Nv = Sv + Ev + Iv. System (3) can be rewritten in the following way dX dt = A(X)X + F (4) with X = (Sh,Vh,Eh,Ih,Rh,Av,Sv,Ev,Iv) T , A(X) = (Aij)1≤i,j≤9 were A1,1 = −(λh + k1), A2,1 = ξ, A2,2 = −(πλh + µh), A3,1 = λh, A3,2 = πλh, A3,3 = −k2, A4,3 = γh, A4,4 = −k3, A5,4 = σ, A5,5 = −µh, A6,7 = A6,8 = A6,9 = µb, A7,6 = θ, A7,7 = −(λv + µv), A8,7 = λv, A8,8 = −k4, A9,8 = γv, A9,9 = −µ2, A6,6 = − ( k6 + µb Sv + Ev + Iv K ) 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 ∈ R9+. Thus, using the fact that F ≥ 0, system (4) is positively invariant in R9+, Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 4 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... which means that any trajectory of the system starting from an initial state in the positive orthant R9+, 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 Ṅh(t) = Λh −µhNh −δIh ≤ Λh −µhNh. (5) So that 0 ≤ Nh(t) ≤ Λh µh + ( Nh(0) − Λh µh ) e−µht. (6) Thus, at t −→∞, 0 ≤ Nh(t) ≤ Λh µh . From the last three equations of model system (3), it also follows that Ṅv(t) = θAv −µvSv −µ1Ev −µ2Iv ≤ θAv −µvNv. (7) So that 0 ≤ Nv(t) ≤ θAv µv + ( Nv(0) − θAv µv ) e−µvt. (8) Thus, at t −→ ∞, 0 ≤ Nv(t) ≤ θK µv since Av ≤ K. Therefore, all feasible solutions of model system (3) enter the region: D = { (Sh,Vh,Eh,Ih,Rh,Av,Sv,Ev,Iv) ∈ R9 : Sh + Vh + Eh + Ih + Rh ≤ Λh µh ; Av ≤ K; Sv + Ev + Iv ≤ θK/µv} . III. THE DISEASE–FREE EQUILIBRIA AND STABILITY ANALYSIS Now let N the net reproductive number [25] given by N = µbθ µv(θ + µA) . (9) We prove the following result Proposition 3.1: a) If N ≤ 1, then, system (3) has only one trivial disease–free equilibrium TE := P0 = ( Λh k1 , ξΛh µhk1 , 0, 0, 0, 0, 0, 0, 0 ) . b) If N > 1, then, system (3) has a Disease–Free Equilibrium P1 = ( S0h,V 0 h , 0, 0, 0,A 0 v,S 0 v, 0, 0 ) , with S0h = Λh k1 , V 0h = ξΛh k1µh , A0v = K ( 1 − 1 N ) , S0v = θ µv K ( 1 − 1 N ) . 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 R0. Precisely, we have proved that when N ≤ 1, model sytem (3) admits only one equi- librium called trivial equilibrium (TE := P0). When N > 1, model sytem (3) admits additionally the disease–free equilibrium (DFE := P1). We prove, in the following, that the trivial equilib- rium 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 R0 < 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 P1 is locally asymptotically stable if R0 < 1 and unstable if R0 > 1, where R0 is the basic reproduction number [26], [82], given by R20 = βhvβvhKθ(πξ + µh)(k3ηh + γh) (µh + ξ)(µh + γh)(µh + δ + σ) × Λhµh(µ2ηv + γv) µvµ2(Λh + mµh) 2(µ1 + γv) ( 1 − 1 N ) . (10) Proof: See appendix B. Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 5 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... B. Global stability analysis 1) Global asymptotic stability of the trivial equilibrium TE := P0: Proposition 3.3: If N ≤ 1, then, TE := P0 is globally asymptotically stable on D. Proof: See Appendix C. 2) Global asymptotic stability of the disease– free equilibrium : Following [30], we prove that the disease–free equilibrium DFE := P1 is glob- ally 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{ ẊS = A1(X)(XS −XDFE) + A12(X)XI ẊI = A2(X)XI (11) where XS is the vector representing the state of different compartments of non transmitting indi- viduals (Sh, Vh, Rh, Av, Sv) and the vector XI represents the state of compartments of different transmitting individuals (Eh, Ih, Ev, Iv). Here, we have XS = (Sh,Vh,Rh,Av,Sv)T , XI = (Eh,Ih,Ev,Iv) T , X = (XS,XI) and XDFE = ( S0h,V 0 h , 0, 0, 0,A 0 v,S 0 v, 0, 0 )T , with A1(X) =   −k1 0 0 0 0 ξ −µh 0 0 0 0 0 −µh 0 0 0 0 0 −K6 K7 0 0 0 θ −µv   , A12(X) =   0 0 −a13 −a14 0 0 0 −a23 −a24 0 0 σ 0 0 0 0 0 κ κ 0 −a41 −a42 0 0 0   , A2(X) =   −k2 0 b13 b14 γh −k3 0 0 b31 b32 −k4 0 0 0 γv −µ2   , with K6 = ( k6 + µb S0v K ) , K7 = µb ( 1 − Av K ) , a13 = βhvηvSh Nh + m , a14 = βhvSh Nh + m , a23 = πβhvηvVh Nh + m , a24 = πβhvVh Nh + m , a41 = βvhηhSv Nh + m , a42 = βvhSv Nh + m , b13 = βhvηvH Nh + m , b14 = βhvH Nh + m , b31 = βvhηhSv Nh + m , b32 = βvhSv Nh + m , κ = µb ( 1 − Av K ) and H = (Sh + πVh). A direct computation shows that the eigenvalues of A1(X) are real and negative. Thus the system ẊS = A1(X)(XS−XDFE) is globally asymptot- ically stable at XDFE. Note also that A2(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: G = { (Sh,Vh,Eh,Ih,Rh,Av,Sv,Ev,Iv) ∈ R9 : Sh ≤ Nh,Vn ≤ Nh,Eh ≤ Nh,Ih ≤ Nh,Rh ≤ Nh, N̄h = Λh/(µh + δ) ≤ Nh ≤ N0h = Λh/µh; Av ≤ K; Sv + Ev + Iv ≤ θK/µv} . Let us recall the following theorem [48] Theorem 3.1: Let G ⊂ U = R5 × R4. The system (11) is of class C1, defined on U. If (1) G is positively invariant relative to (11). (2) The system ẊS = A1(X)(XS − XDFE) is Globally asymptotically stable at XBRDFE. (3) For any x ∈G, the matrix A2(X) is Metzler irreducible. (4) There exists a matrix Ā2 , which is an upper bound of the set M = {A2(x) ∈M4(R) : x ∈G} with the property that if A2 ∈ M, for any x̄ ∈ G, such that A2(x̄) = Ā2, then x̄ ∈ R5 ×{0}. (5) The stability modulus of Ā2, α(A2) = max λ∈sp(A2) Re(λ) satisfied α(A2) ≤ 0. 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 Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 6 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... by Ā2 =   −k2 0 b̄13 b̄14 γh −k3 0 0 βvhηhS̄v N̄h + m βvhS̄v N̄h + m −k4 0 0 0 γv −µ2   , where b̄13 = βhvηv(S̄h + πV̄h) N̄h + m , b̄14 = βhv(S̄h + πV̄h) N̄h + m , S̄h = S0h, V̄h = V 0 h , Āv = K, S̄v = θ µv K, and N̄h = Λh (µh + δ) . To check condition (5) in theorem 3.1, we will use the following useful lemma [48] which is the a characterization of Metzler stable matrices: Lemma 3.1: Let M be a square Metzler matrix written in block form ( A B C D ) with A and D square matrices. M is Metzler stable if and only if matrices A and D −CA−1B are Metzler stable. A necessary condition for a Metzler matrix to be stable is that the elements on the diagonal are negative. Note also that A is a Metzler stable matrix is equivalent to A is invertible and −A−1 ≥ 0. Lemma 3.1 allows to reduce the problem of Metzler stability, by induction, to the stability of 2 × 2 Metzler matrices [48]. In our case, we have A = ( −k2 0 γh −k3 ) , B =   βhvηv(S̄h + πV̄h)N̄h + m βhv(S̄h + πV̄h)N̄h + m 0 0  , C =   βvhηhS̄vN̄h + m βvhS̄vN̄h + m 0 0  , and D = ( −k4 0 γv −µ2 ) . Clearly, A is a stable Metzler matrix. Then, after some computations, we obtain D − CA−1B is a stable Metzler matrix if and only if R2G ≤ 1 (12) where R2G = βhvβvhKθΛh(ηvµ2 + γv)(k3ηh + γh) µvµ2µhk1k2k3k4 × (µh + πξ)(µh + δ) 2 (Λh + m(µh + δ)) 2 Thus we claim the following result Theorem 3.2: If N > 1 and R2G ≤ 1, then the disease–free equilibrium P1 is globally asymptot- ically stable in G. Remark 3.1: Note that R2G = R 2 0 (µh + δ) 2(Λh + mµh) 2 µ2h(Λh + m(µh + δ)) 2 ( N N − 1 ) and condition (12) is equivalent to R20 ≤ ( N − 1 N ) µ2h (µh + δ) 2 (Λh + m(µh + δ)) 2 (Λh + mµh) 2 (13) In absence of disease–induced death in human (δ = 0), inequality (13) becomes R20 ≤ ( N − 1 N ) < 1. (14) This shows that with the establishment of an effective treatment, it is possible to have R0 and RG less than 1. Theorem (3.2) means that for R0 < RG < 1, the DFE is the unique equilibrium (no co-existence with an endemic equilibrium). If R0 ∈ [RG, 1], then it is possible to have co-existence with en- demic 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]. IV. THE ENDEMIC EQUILIBRIA AND BIFURCATION ANALYSIS A. Existence of endemic equilibria We now turn to study the existence of an en- demic equilibrium of model system (3). Let R0 the basic reproduction number [26], [82] given by Eq. (10). we claim the following Proposition 4.1: Let N > 1 and µv ≤ µ1 ≤ µ2. Then Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 7 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... (i) There exists at least one endemic equilibrium whenever R0 > 1. (ii) The backward bifurcation phenomenon may occurs when R0 ≤ 1. (iii) The disease–induced death is responsible of the backward bifurcation phenomenon. (iv) In the absence of the disease–induced death (δ = 0 and µv = µ1 = µ2), system (4) have a unique endemic equilibrium whenever R0 > 1, and the backward bifurcation phe- nomenon not occurs whenever R0 ≤ 1 (See remark 4.1). Proof: See appendix D. The backward bifurcation phenomenon, in epi- demiological systems, indicate the possibility of existence of at least one endemic equilibrium when R0 is less than unity. Thus, the classical requirement of R0 < 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 (For- ward) bifurcation when R0 = 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 pop- ulations, the disease–free equilibrium is the unique equilibrium whenever N > 1 and R0|δ=0,µv=µ1=µ2 < 1. Using the direct Lya- punov method, we prove the global asymptotic stability of the disease–free equilibrium whenever R0|δ=0,µv=µ1=µ2 < 1. Remark 4.1: Assumed that N > 1. Let k7 = Λh + mµh, k8 = πξ + µh, k11 = k3ηh + γh and R1 = R0|δ=0,µv=µ1=µ2 . In the absence of disease-induced death, i.e, δ = 0 and µv = µ1 = µ2, Eq. (44) (see appendix D) becomes λ∗h [ B02(λ ∗ h) 2 + B01λ ∗ h + B00 ] = 0 (15) with B02 = k2k3k27πµv + βvhk7k11Λhµhπ > 0, B00 = k1k2k3k 2 7µhµv ( 1 −R21 ) and B01 = k1k2k3k 2 7µvπ ( 1 −µhR21 ) + k2k3k 2 7µhµv +βvhk7k8k11Λhµh. Equation (15) have only one positive solution whenever R1 > 1. If R1 ≤ 1, then coefficients B00, B01, B02 are always positive, and the disease- free equilibrium is the unique equilibrium. From this we conclude that the disease–induced mor- tality is the possible cause for the occurrence of multiple endemic equilibria below the classical threshold R1 = 1. The global stability of the disease–free equilib- rium may be achieved by Lyapunov method. To this aim, let us consider the following Lyapunov function [37], [40] Y = 4∑ i=1 giIi where I = (Eh,Ih,Ev,Iv) and gi, i = 1, . . . , 4 are positive weights given by g1 = 1; g2 = k2 (k3ηh + γh) , g3 = k2k3(N 0 h + m) βvhS 0 v(k3ηh + γh) , g4 = βhv [ S0h + πV 0 h ] µ2(N 0 h + m) . Along the solutions of (3) we have: Ẏ = 4∑ i=1 giİi = g1Ėh + g2İh + g3Ėv + g4İv = g1 [λh [Sh + (1 − �)Vh] − (µh + γh)Eh] +g2 [γhEh − (µh + δ + σ)Ih] +g3 [λvSv − (µ1 + γv)Ev] + g4 (γvEv −µ2Iv) = ( g1 βhvηv [Sh + πVh] N0h + m + g4γv −g3k4 ) Ev + ( g1 βhv [Sh + πVh] N0h + m −g4µ2 ) Iv + ( g3 βvhηhSv N0h + m + g2γh −g1k2 ) Eh + ( g3 βvhSv N0h + m −g2k3 ) Ih After replacing the constants gi, i = 1, . . . , 4 by their value, and using the fact that Sh ≤ S0h, Vh ≤ V 0h , Av ≤ A 0 v, and Sv ≤ S0v in D1 = {(Sh,Vh,Eh,Ih,Rh,Av,Sv,Ev,Iv) ∈D : Sh ≤ S0h,Vh ≤ V 0 h ,Av ≤ A 0 v,Sv ≤ S0v } , Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 8 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... it follows that Ẏ ≤ ( g1 βhvηv [ S0h + πV 0 h ] N0h + m + g4γv −g3k4 ) Ev = k2k3k4(N 0 h + m) βvhS 0 v(k3ηh + γh) ( R21 − 1 ) Ev We have Ẏ ≤ 0 if R1 ≤ 1, with Ẏ = 0 if R1 = 1 or Ev = 0. Whenever Ev = 0, we also have Eh = 0, Ih = 0 and Iv = 0 . Substituting Eh = Ih = Ev = Iv = 0 in the first, second, fifth, sixth and seventh equations of Eq. (3) with δ1 = 0 gives Sh(t) → S0h, Vh(t) → V 0 h , Rh(t) → 0, Av(t) → A0v, Sv(t) → S0v as t →∞. Thus [Sh(t),Vh(t),Eh(t),Ih(t),Rh(t),Av(t),Sv(t),Ev(t) ,Iv(t)] → (S0h,V 0 h , 0, 0, 0,A 0 v,S 0 v, 0, 0) as t →∞. It follows from the LaSalle’s invariance principle [45] that every solution of (3) (when R1 ≤ 1), with initial conditions in D1 converges to P1, as t → ∞. Hence, the DFE (P1), of model (3), is GAS in D1 if R1 ≤ 1. B. Bifurcation analysis Previous Analysis state that multiple endemic equilibria may occur althougt R0 < 1. In order to better investigate the variation of model’s pre- diction as R0 varied, we perform a bifurcation analysis at the criticality, i. e. at the Disease– free Equilibrium DFE := P1 and R0 = 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?u2 + b?$u, where, u represent a real function of real variable, a? = v 2 ·Dxxf(x0,$)w2 ≡ ≡ n∑ k,i,j=1 vkwiwj ∂2fk(0, 0) ∂xi∂xj (16) and b? = v ·Dx$f(x0,$)w ≡ n∑ k,i=1 vkwi ∂2fk(0, 0) ∂xi∂$ . (17) Note that the symbol $ denotes a bifurcation parameter to be chosen, fi’s denotes the right hand side of system (3), x denotes the state vector, x0 the Disease–free Equilibrium P1; v and w denote the left and right eigenvectors, respectively, cor- responding 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 R0 = 1 we get the critical value β∗hv = µvµ2k1k2k3k4(Λh + mµh) 2 ( N N − 1 ) βvhΛhµhKθ(πξ + µh)(µ2ηv + γv) [ηhk3 + γh]) . Note also that in fk(0, 0), the first zero corresponds to the disease–free equilibrium, P1, for the system (3). Since βhv = β∗hv is the bifurcation parameter, it follows from $ = βhv − β∗hv that $ = 0 when βhv = β ∗ hv which is the second component in fk(0, 0). The Jacobian matrix of the system (4) evaluated at the disease–free equilibrium P1 with βhv = β∗hv is given by J(P1) =   −k1 0 0 0 0 ξ −µh 0 0 0 0 0 −k2 0 0 0 0 γh −k3 0 0 0 0 σ −µh 0 0 0 0 0 0 0 − βvhηhS 0 v H − βvhS 0 v H 0 0 βvhηhS 0 v H βvhS 0 v H 0 0 0 0 0 0 Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 9 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... 0 0 − β∗hvηvS 0 h H − β∗hvS 0 h H 0 0 − β∗hvπηvV 0 h H − β∗hvπV 0 h H 0 0 β∗hvηvS0 H β∗hvS0 H 0 0 0 0 0 0 0 0 −K1 K2 K2 K2 θ −µv 0 0 0 0 −k4 0 0 0 γv −µ2   , where we have set H = N0h + m, K1 = µbθ µv and K2 = k6µv θ . The eigenvalues of the Jacobian matrix J(P1) are λ1 = −µh, λ2 = −k1, and the other eigenvalues are the eigenvalue of the following matrix J̄ =  −k2 0 0 0 β∗ hv ηvS0 H β∗ hv S0 H γh −k3 0 0 0 0 0 0 −K1 K2 K2 K2 − βvhηhS 0 v H − βvhS 0 v H θ −µv 0 0 βvhηhS 0 v H βvhS 0 v H 0 0 −k4 0 0 0 0 0 γv −µ2   . The characteristic polynomial of J̄ is given by f(λ) = λ6 +a5λ 5 +a4λ 4 +a3λ 3 +a2λ 2 +a1λ+a0 (18) with a0 = − k1k2k3k4k 2 7µ2µbµv(k6µv −µbθ) k1k 2 7µbµv ( 1 −R20 ) . The others coefficients a5, a4, a3, a2, and a1 are obtained after computations on Maxima software [58], [89]. Since at the criticality, we have R0 = 1, then a0 = 0, thus equation (18) becomes f(λ) = λ ( λ5 + a5λ 4 + a4λ 3 + a3λ 2 + a2λ + a1 ) . Then, the Jacobian J(P1) of the linearized system has a simple zero eigenvalue and therefore P1 is a nonhyperbolic equilibrium for R0 = 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(P1) is given by w = (w1,w2,w3,w4,w5,w6,w7,w8,w9) T where w1 = − β∗hvk9µhΛh k21γvk7 w9 < 0, w2 = − ξΛhβ ∗ hvk9(µh + k1π) k21µhγvk7 w9 < 0, w3 = β∗hvΛhk9(µh + k1π) k1k2k7γv w9 > 0, w4 = β∗hvΛhγhk9(µh + k1π) k1k2k3k7γv w9 > 0, w5 = β∗hvΛhσγhk9(µh + k1π) k1k2k3k7µhγv w9 > 0, w7 = − βvhµhKθ µvk7 ( 1 − 1 N ) (ηhw3 + w4) < 0, w8 = βvhµhKθ µvk4k7 ( 1 − 1 N ) (ηhw3 + w4) > 0, w6 = µb k6N2 (w7 + w8 + w9) and w9 > 0. Similarly, J(P1) has a left eigenvector v = (v1,v2,v3,v4,v5,v6,v7,v8,v9) where v1 = v2 = v5 = v6 = v7 = 0, v3 = µvk1k7 βvhΛhk8 v9, v4 = βvhKθµh(ηvµv + γv) k3k4k7µv ( 1 − 1 N ) v9, v8 = (ηvµv + γv) k4 v9 and v9 > 0. a) Computation of a?: Using the non–zero components of v and the associated non–zero partial derivatives of f (at the DFE P1), for system (3), we obtain a? = 1 2 v3 9∑ i,j=1 wiwj ∂2f3(0, 0) ∂xi∂xj + 1 2 v8 9∑ i,j=1 wiwj ∂2f8(0, 0) ∂xi∂xj . We finally obtain (See the details in appendix E) a? = φ1 −φ2 Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 10 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... where φ1 = 1 2 v3 { β∗ hv µh k1(Λh + mµh) 2 [(�ξΛh + mµh) (w1 + π w2) ×(w8ηv + w9 + ηv + 1) −2Λh(µh + πξ) (w3 + w4 + w5) (w8ηv + w9)]} − 1 2 v8 βvhµ 2 h Kθ µv(Λh + µhm) 2 ( 1 − 1 N ) w5 (ηhw3 + w4) + 1 2 v8 βvhµh (Λh + µhm) (ηhw3 + w4) w7 − 1 2 v8 βvhµ 2 h Kθ ( 1 − 1 N ) µv(Λh + µhm) 2 × [ 2(ηh + 1)w3w4 + 2 ( ηhw 2 3 + w 2 4 )] < 0 and φ2 = 1 2 v8 βvhµ 2 h Kθ µv(Λh + µhm) 2 ( 1 − 1 N ) × [(ηhw3 + w4) (w1 + w2)] < 0 b) Computation of b?: b? = v3 Λh(µh + πξ) k1(Λh + µhm) (ηvw8 + w9) > 0. Since b? > 0 according to the sign of wi,vi, 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 R0 = 1. If the reversed inequality holds, then the bifurcation at R0 = 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 illus- trative purpose only, and may not necessarily be realistic epidemiologically): —Using the parameters values in Table II, ex- cept µv = µ1 = µ2 = 1/14, Λh = 200, � = 0.80, ξ = 0.475, δ = 0.6, β̃hv = 6, β̃vh = 50 and K = 1000 such that R0 = 0.6095 < 1 and a? = 1.0348 × 10−5 > 0, the numerical resolu- tion of equation (44) (see appendix A), gives the following solution: λ∗1h = 0, λ ∗ 2h = 0.0083, λ ∗ 3h = 10.9412, λ∗4h = −0.0080 and λ ∗ 5h = −0.0001; Note that the first solution λ∗1h = 0 corresponds to the disease free equilibrium. The second, and third solution, λ∗2h = 0.0083, λ ∗ 3h = 10.9412, correspond to endemic equilibria; λ∗2h = 0.0083 correspond to unstable endemic equilibrium and λ∗3h = 10.9412 corresponds to the stable endemic equilibrium. The fourth and fifth solution λ∗4h = −0.0080 and λ∗5h = −0.0001 are not biologically meaningful. 0 20 40 60 80 100 0 100 200 300 400 time (days) In fe c ti o u s H u m a n s , Ih (t ) Fig. 1. Time profile of infectious humans using different initial conditions showing that the equilibrium λ∗3h = 10.9412 is stable even if R0 = 0.6095 < 1 . —Using the parameters values in Table II, ex- cept µv = µ1 = µ2 = 1/14, Λh = 100, � = 0.80, ξ = 0.475, δ = 0.6, β̃hv = 4.0385, β̃vh = 100 and K = 1000 such that R0 = 1 and a? = 2.3665×10−4 > 0 , the numerical resolution of equation (44), gives the following solution: λ∗11h = 0, λ ∗ 22h = 0.0114, λ ∗ 3h = 8.5310, and λ∗44h = −0.0111; The first solution λ ∗ 1h = 0 corre- sponds to the disease free equilibrium. The second, and third solution, λ∗2h = 0.0083, λ ∗ 33h = 8.5310, correspond to endemic equilibria; λ∗22h = 0.0114 correspond to unstable endemic equilibrium and λ∗33h = 8.5310 corresponds to the stable endemic equilibrium. The fourth solution λ∗4h = −0.0111 is not biologically meaningful. —In the absence to disease induced death (δ = 0) and choosing β̃hv = 4.0188 and K = 1000 such that R0 = 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. —Choosing β̃hv = 10 and K = 1000 such that R0 = 1.630976 > 1 and a? = −1.8011 < 0, equation (44) admit only one positive solution given by λ∗1h = 0.0001, which correspond to the Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 11 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... 0 20 40 60 80 100 0 50 100 150 200 250 300 time (days) In fe c ti o u s H u m a n s , Ih (t ) Fig. 2. Time profile of infectious humans using different initial conditions showing that the equilibrium λ∗33h = 8.5310 is stable even if R0 = 1 . endemic equilibria when the basic reproduction number, R0, is greater than 1. To conclude, depending to the values of param- eters of model (3), the phenomenon of backward bifurcation may occurs when the classical basic reproduction number R0 is less than unity. V. THRESHOLD ANALYSIS AND VACCINE 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). Ṡh = Λh −λhSh −µhSh Ėh = λhSh − (µh + γh)Eh İh = γhEh − (µh + δ + σ)Ih Ṙh = σIh −µhRh Ȧv = µb ( 1 − Av K ) (Sv + Ev + Iv) − (θ + µA)Av Ṡv = θAv −λvSv −µvSv Ėv = λvSv − (µ1 + γv)Ev İv = γvEv −µ2Iv (19) with λh and λv defined at (1) and (2), respectively. Following procedure in [26], [82], the correspond- ing basic reproduction number of model (19), Rs, is given by R2s = βhvβvhKθ(k3ηh + γh)(γv + ηvµ2) µvµ2(µh + γh)(µh + δ + σ)(µ1 + γv) × Λhµh (Λh + mmuh) 2 ( 1 − 1 N ) (20) So we deduce that Rvac := R0 = Rs √ (πξ + µh) (µh + ξ) . (21) From Eq. (21), it follows that, in the absence of vaccination (ξ = 0) or when the vaccine efficacy is very low (� → 0), we have Rvac = Rs. 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 ar- boviruses infection, humans vaccination can con- tribute 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 simu- lations, 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 R0 plotted as function of the vaccine efficacy parameter � and the proportion Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 12 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... TABLE III BASELINE VALUES OF PARAMETERS OF MODEL (3) AND THEIR SOURCES. Parameter Baseline value Sources Λh 2.5 day −1 [40] ξ variable � Variable ηh, ηv 0.35 Assumed ˜βhv 0.75 day −1 [40] ˜βhv 0.75 day −1 [40] γh 1/3 day −1 [30] γv 1/2 day−1 [30] µh day−1 (71 × 365) [68] µv (1/14) day−1 [40] µA 1/5 day−1 [30] µ−11 10 days [30] µ−12 5 days [30] θ 0.08 day−1 [30], [68] δ 10−3 day−1 [40] σ 0.1428 day−1 [2], [40] a 1 day−1 [40], [61] m 100 Assumed K 2 × 5000 Assumed µb 6 day −1 [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 R0, and the 0 50 100 150 200 0 200 400 600 800 1000 time (days) T o ta l n u m b e r o f In fe c te d h u m a n s ,E h (t )+ I h (t ) ξ=0.05 ξ=0.25 ξ=0.5 xi=0.75 xi=1 0 50 100 150 200 0 500 1000 1500 2000 2500 time (days) T o ta l n u m b e r o f In fe c te d v e c to rs ξ=0.050 ξ=0.25 ξ=0.5 ξ=0.75 ξ=1 Fig. 3. Total number of infected humans and vectors varying the proportion of susceptible huamans vaccinated ξ = (0.05; 0.25; 0.5; 0.75; 1) with a vaccine simulating 60% of effectiveness (i.e. � = 0.60 or π = 1 − � = 0.4). 0 50 100 150 200 0 500 1000 1500 2000 2500 time (days) T o ta l n u m b e r o f In fe c te d v e c to rs ε=0.25 ε=0.50 ε=0.60 ε=0.70 ε=0.80 ε=0.90 ε=1 0 20 40 60 80 100 0 200 400 600 800 time (days) T o ta l n u m b e r o f In fe c te d h u m a n s ε=0.25 ε=0.50 ε=0.60 ε=0.70 ε=0.80 ε=0.90 ε=1 Fig. 4. Infected humans and Vector varying the efficacy level of the vaccine � = (0.25; 0.50; 0.80; 0.90; 1) and considering that 85% of susceptible humans is vaccinated. prevalence of the disease is directly related to the infected states, especially for sizes of Eh(t), Ih(t), Ev(t) and Iv(t). These variables are relevant to Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 13 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... 0 0.5 1 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 ξ ε R 0 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 Table III. 0 50 100 150 200 250 300 350 400 0 200 400 600 800 time (days) T o ta l n u m b e r o f In fe c te d h u m a n s Without vaccination With vaccination:ξ=0.85, ε=0.60 With vaccination:ξ=0.85, ε=0.80 With vaccination:ξ=0.85, ε=1 0 50 100 150 200 250 300 350 400 0 500 1000 1500 2000 2500 time (days) T o ta l n u m b e r o f In fe c te d v e c to rs Without vaccination With vaccination:ξ=0.85, ε=0.60 With vaccination:ξ=0.85, ε=0.80 With vaccination:ξ=0.85, ε=1 Fig. 6. Time profile of total number of infected human and vector without vaccination and with vaccination. the individuals (humans and vectors) who have some life stage of arboviruses in their bodies. The number of infectious humans, Ih, 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- TABLE IV PARAMETER VALUE RANGES OF MODEL (3) USED AS INPUT FOR THE LHS METHOD. Parameter Range Parameter Range Λh [1 , 6 ] µA [1/10,1/4] ξ [0.05,1] µ1 [1/21,1/3] � [0.5,0.9] µ2 [1/21,1/3] ηh, ηv [0.1,0.8] θ [0.01,0.17] ˜βhv [0.375,1] δ [10 −5,10−2] ˜βvh [0.375,1] σ [0.1428,1/3] γh [1/12,1/2] a [1,3] γv [1/21,1/2] m [1,201] µh [ 1 78 × 365 , 1 45 × 365 ] K 103×[10,15] µv [1/21,1/10] µb [6,18] TABLE V INITIAL CONDITIONS. Human Initial value vector Initial value Sh: 1000 Av 1000 Vh: 0 Sv: 500 Eh: 20 Ev: 20 Ih: 10 Iv: 40 Rh: 0 tion, use the parameters values of two particular arboviruses, Dengue and Chikungunya. It is impor- tant 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 R0: We study the impact of each parameter of the model on the value of the basic reproduction number R0. Fol- lowing 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 repro- duction number, R0. Table III troughly estimates Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 14 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... 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 pa- rameter values sampled using the Latin Hypercube Sampling (LHS) method were used to run 1,000 simulations. With these 1,000 runs of Latin Hypercube Sam- pling, the derived sampling distribution of R0 is shown in Figure 7. From this sampling we get that the mean of R0 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 R0 > 1. 0 2 4 6 8 10 12 14 16 18 0 100 200 300 400 500 600 R 0 F re q u e n c ie s 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. From the previous sampling we compute the Partial Rank Correlation Coefficients between R0 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 R0. 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 R0. The other parameters with an im- portant effect are θ, a, Λh and µ2. The parameters TABLE VI PRCC BETWEEN R0 AND EACH PARAMETER. Parameter Correlation Coefficients P–values 1 Λh *–0.6067 1.4578E−99 2 ξ 0.0529 0.0977 3 � ***–0.8043 2.6732E−223 4 ηh 0.2879 4.0576E−20 5 ˜βhv 0.4354 1.3609E−46 6 γh –0.2598 1.4099E−16 7 µh 0.2526 9.9492E−16 8 δ –0.0386 0.2274 9 σ –0.3269 7.7785E−26 10 ηv 0.2039 1.1635E−10 11 ˜βvh 0.4215 1.7130E−43 12 γv 0.2117 2.1787E−11 13 µv –0.3029 3.0015E−22 14 µA –0.0121 0.7049 15 µ1 –0.2948 4.2501E−21 16 µ2 *–0.5087 1.2669E−65 17 θ **0.7626 3.0823E−187 18 a **0.7134 3.4096E−153 19 m –0.0436 0.1724 20 K 0.3880 1.4683E−36 21 µb 0.0082 0.7973 which do not have almost any effect on R0 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 Eh(t), Ih(t), Ev(t), and Iv(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 follow- ing facts: –For the value of Eh, the parameters with more influence are θ, K, a, �, Λh and µ2. The parameters which do not have almost any effect on the variation of Eh 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 Ih, the parameters with more influence are Λh and γh. The least sensitive pa- Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 15 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... TABLE VII PRCC BETWEEN INFECTED HUMANS Eh AND EACH PARAMETER. Parameter Correlation Coefficients P–values 1 Λh **0.6842 3.2080E−136 2 ξ 0.4115 2.4590E−41 3 � ***0.7177 6.8762E−156 4 ηh –0.2457 6.1306E−15 5 β̃hv –0.4215 1.7187E−43 6 γh 0.2172 6.2865E−12 7 µh 0.0086 0.7879 8 δ -0.0259 0.4176 9 σ 0.3395 7.4246E−28 10 ηv –2378 4.5858E−14 11 β̃vh –0.4232 7.4972E−44 12 γv –0.2311 2.4083E−13 13 µv 0.2906 1.5881E−20 14 µA 0.0210 0.5122 15 µ1 0.3340 5.8090E−27 16 µ2 *0.5747 3.1691E−87 17 θ ***–0.7599 3.7832E−185 18 a ***–0.7597 4.9923E−185 19 m 0.0537 0.0931 20 K ***–0.7477 4.2124E−176 21 µb –0.0068 0.8328 TABLE VIII PRCC BETWEEN INFECTIOUS HUMANS Ih AND EACH PARAMETER. Parameter Correlation Coefficients P–values 1 Λh ***0.8727 9.1342E−307 2 ξ 0.0078 0.8062 3 � –0.2887 2.8614E−20 4 ηh 0.0711 0.0261 5 β̃hv 0.0850 0.0078 6 γh ***–0.8722 5.9181E−306 7 µh –0.0363 0.2566 8 δ 0.0412 0.1978 9 σ –0.0531 0.0965 10 ηv 0.0310 0.3316 11 β̃vh 0.1297 4.6364E−5 12 γv –0.0179 0.5764 13 µv –0.0544 0.0886 14 µA –0.0222 0.4877 15 µ1 –0.0580 0.0697 16 µ2 –0.0423 0.1855 17 θ 0.1312 3.7931E−5 18 a 0.1428 2.8933E−6 19 m –0.0017 0.9586 20 K 0.1783 1.9260E−8 21 µb –0.0054 0.8648 TABLE IX PRCC BETWEEN INFECTED VECTORS Ev AND EACH PARAMETER. Parameter Correlation Coefficients P–values 1 Λh –0.0186 0.5603 2 ξ –0.0111 0.7280 3 � 0.0135 0.6723 4 ηh –0.1086 6.6203E−4 5 β̃hv –0.0664 0.0375 6 γh 0.0560 0.0798 7 µh –0.0295 0.3563 8 δ 0.0116 0.7170 9 σ 0.0734 0.0215 10 ηv –0.0273 0.3928 11 β̃vh –0.0913 0.0043 12 γv 0.0069 0.8282 13 µv **-0.5923 7.6830E−94 14 µA 0.0157 0.6235 15 µ1 0.0331 0.3006 16 µ2 0.0043 0.8933 17 θ ***0.9225 0 18 a –0.0822 0.0100 19 m 0.0027 0.9324 20 K ***0.9199 0 21 µb 0.1125 4.1594E−4 TABLE X PRCC BETWEEN INFECTIOUS VECTORS Iv AND EACH PARAMETER. Parameter Correlation Coefficients P–values 1 Λh 0.2254 9.3729E−13 2 ξ –0.0090 0.7785 3 � –0.1228 1.1697E−4 4 ηh 0.3126 1.1661E−23 5 β̃hv 0.0031 0.9216 6 γh –0.3233 2.7921E−25 7 µh 0.0381 0.2338 8 δ –0.0215 0.5015 9 σ –0.3869 2.4025E−36 10 ηv 0.0196 0.5402 11 β̃vh 0.5584 2.0585E−109 12 γv –0.6287 6.0859E−109 13 µv –0.4856 4.0722E−59 14 µA 0.0294 0.3583 15 µ1 –0.4380 3.3922E−47 16 µ2 –0.0103 0.7470 17 θ **0.8728 7.6088E−307 18 a *0.6011 2.5895E−97 19 m –0.0640 0.0451 20 K *0.8600 5.9602E−288 21 µb –0.0770 0.0159 Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 16 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... rameters is µb, the number of eggs at each deposit per capita; –For the value of Ev, 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 Iv, 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 R0 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 ob- tained in the previous sections, we provide here some numerical simulations. We use the nonstan- dard 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 Mat- lab [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:  Xk+1S −X k S φ(h) = A1(Xk)(XkS −XDFE) −D12(XkI )X k+1 S + B12(X k)XkI Xk+1I −X k I φ(h) = A2(Xk+1S )X k I (22) such that −D12(XI)XS + B12(X)XI = A12(X)XI (23) with D12(XI) =   λh 0 0 0 0 0 πλh 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 λv   , and B12(X) =  0 0 0 0 0 0 0 0 0 0 0 σ 0 0 0 0 0 µb ( 1 − A0v K ) µb ( 1 − Av K ) 0 0 0 0 0 0   , which implies that the scheme is consistant with formulation (11). Rearranging (22), we obtain the foollowing new expression { AkXk+1 = Bk Xk ≥ 0 (24) with Ak = ( I5 + φ(h)D12(X k I ) 0 0 I5 ) and Bk =( XkS + φ(h) [ A1(Xk)(XkS −XDF E) + B12(X k)XkI ] XkI [ I4 + φ(h)A2(Xk+1S ) ] ) . where I4 and I5 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 Page 17 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... Lemma 7.1: Our non-standard numerical scheme (22) is positively stable, i.e. for Xk ≥ 0 we obtain Xk+1 ≥ 0, where Xk = ( Skh,V k h ,E k h,I k h,R k h,A k v,S k v ,E k v ,I k v )T . Proof: We suppose Xk ≥ 0. Ak is a positive diagonal matrix and thus, A−1k ≥ 0. B12 is a pos- itive matrix and we also have −A1(Xk)XDFE ≥ 0. To show that Bk is a positive matrix, it suffices to choose φ(h) such that Id + φ(h)A1(X) ≥ Id + φ(h)A1 ≥ 0, Id + φ(h)A2(X) ≥ Id + φ(h)A2 ≥ 0 where A1 and A2 are lower bounds for the sets {X ∈ D|A1(X)} and {X ∈ D|A2(X)} respec- tively. Following [30], to have Bk ≥ 0, it suffices to consider the following time-step function φ(h) = 1 −e−Qh Q (25) with Q ≥ max (k1,k2,k3,µh,k4,k6,µv,µ2). We have proved that Xk ≥ 0 implies Xk+1 ≥ 0. 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 TE := P0 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 R0 < 1, then the disease–free equilibrium DFE := P1 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 disease– free equilibrium is stable (first panel). The phe- nomenon of backward bifurcation occurs in the second panel of figure 9, where the stable disease- free equilibrium of the model co–exists with a stable endemic equilibrium when the associated re- production number, R0, is less than unity. Figures 10–11 show the existence of at least one endemic equilibrium whenever R0 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 vacci- nation 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 mortal- ity may be the main cause for the occurrence of the backward bifurcation (see remark 4.1). This means that relatively high values of disease– induced mortality rate may induce stable endemic states also when the basic reproduction number R0 is below the classical threshold R0 = 1. The stability analysis reveals that for N ≤ 1, Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 18 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... 0 200 400 600 800 1000 1200 1400 1600 1800 2000 0 1000 2000 3000 4000 5000 6000 times (days) N o n i n fe c te d H u m a n S h V h R h 0 100 200 300 400 0 200 400 600 800 1000 times (days) V e c to r p o p u la ti o n A v S v E v I v Fig. 8. Time profile of both population without vector (with θ = 0.0008, so N = 0.2679 < 1. In this case the trivial equilibrium is globally asymptotically stable. 0 500 1000 1500 2000 0 200 400 600 800 1000 1200 1400 time (days) In fe c te d h u m a n s a n d v e c to rs E h (t) I h (t) E v (t) I v (t) 0 500 1000 1500 2000 0 200 400 600 800 1000 time (days) In fe c te d h u m a n s a n d v e c to rs E h (t) I h (t) E v (t) I v (t) Fig. 9. Time profile infected humans and vectors. First panel:R0 =: 0.2377 < 1 and Second panel: R0 = 0.9405 < 1. The backward bifurcation phenomenon is illustrate in second panel. the trivial equilibrium is globally asymptotically stable. When N > 1 and R0 < 1, the disease– free equilibrium is locally asymptotically stable. In 0 50 100 150 200 250 300 350 400 0 100 200 300 400 500 600 700 time (days) In fe c te d h u m a n s E h (t) I h (t) 0 50 100 150 200 250 300 350 400 0 200 400 600 800 time (days) In fe c te d h u m a n s E h (t) I h (t) Fig. 10. Time profile of infected humans with β̃hv = 42.9631, Λh = 20, so that R0 = 1 (first panel) and β̃hv = β̃vh = 20, Λh = 20, so that R0 = 3.5233 > 1 (second panel). 0 2000 4000 6000 8000 10000 0 200 400 600 800 1000 time (days) In fe c te d v e c to rs E v I v 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10 4 0 200 400 600 800 1000 1200 1400 1600 time (days) In fe c te d v e c to rs E v I v Fig. 11. Time profile of infected vectors with β̃hv = 42.9631, Λh = 20, so that R0 = 1 (first panel) and β̃hv = β̃vh = 20, Λh = 20, so that R0 = 3.5233 > 1 (second panel). the absence of disease-induced death, the disease– free equilibrium is also globally asymptotically Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 19 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... 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 R1 exceeds unity. Taking as cases study the dengue and chikun- gunya transmission, we used parameter values from the literature to estimate statistically the basic reproduction number, R0, and to perform a global sensitivity analysis on the basic repro- duction number and infected states (Eh, Ih, Ev, Iv). Using Latin Hypercube Sampling, we obtain that the mean of R0 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 (R0) and on the number of infected humans and vectors (Eh,Ih,Ev,Iv). 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 vaccina- tion 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 conse- quences, so condition (a) does not hold. Moreover, according to previous analysis and french labo- ratory 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 Nk- ague 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 con- tributed 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   Λh −λhSh − (ξ + µh)Sh = 0 ξSh − (1 − �)λhVh −µhVh = 0 λh [Sh + (1 − �)Vh] − (µh + γh)Eh = 0 γhEh − (µh + δ + σ)Ih = 0 σIh −µhRh = 0 µb ( 1 − Av K ) (Sv + Ev + Iv) − (θ + µA)Av = 0 θAv −λvSv −µvSv = 0 λvSv − (µ1 + γv)Ev = 0 γvEv −µ2Iv = 0 (26) To this aim, let P∗ = (S∗h,E ∗ h,I ∗ h,R ∗ h,A ∗ v,S ∗ v,E ∗ v,I ∗ v ) represents any arbitrary endemic equilibrium of (3). Further, let λ∗h = βhv(ηvE ∗ v + I ∗ v ) (N∗h + m) , λ∗v = βvh(ηhE ∗ h + I ∗ h) (N∗h + m) , (27) be the forces of infection of humans and vectors at steady state, respectively. Solving the first five Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 20 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... equations in (26) at steady state gives S∗h = Λh k1 + λ ∗ h , V ∗h = ξΛh (k1 + λ ∗ h)(πλ ∗ h + µh) , E∗h = Λhλ ∗ h(πξ + µh + πλ ∗ h) k2(k1 + λ ∗ h)(πλ ∗ h + µh) , I∗h = γhΛhλ ∗ h(πξ + µh + πλ ∗ h) k2k3(k1 + λ ∗ h)(πλ ∗ h + µh) , R∗h = σγhΛhλ ∗ h(πξ + µh + πλ ∗ h) µhk2k3(k1 + λ ∗ h)(πλ ∗ h + µh) . (28) where π = 1 − �, k1 = µh + ξ, k2 = µh + γh and k3 = µh + σ + δ. Solving the last three equations in (26) at steady state gives S∗v = θA∗v (µv + λ∗v) , E∗v = θA∗vλ ∗ v k4(µv + λ∗v) , I∗v = γvθA ∗ vλ ∗ v µ2k4(µv + λ∗v) . (29) where k4 = µ1 + γv. Substituting (29) in the sixth equation of (26) gives A∗v { µbθ µ2k4 ( 1 − A∗v K )( µ2k4 + k5λ ∗ v µv + λ∗v ) −k6 } = 0 (30) with k5 = µ2 + γv and k6 = θ + µA. The trivial solution of (30) is A∗v = 0. Substitut- ing this solution in (29) gives S∗v = E ∗ v = I ∗ v = 0. When E∗v = I ∗ v = 0, we also have λ ∗ h = 0, thus E∗h = I ∗ h = R ∗ h = 0, S ∗ h = Λh k1 and V ∗h = ξΛh µhk1 . Then we obtain the trivial equilibrium P0 =( Λh k1 , ξΛh µhk1 , 0, 0, 0, 0, 0, 0, 0 ) . Now we suppose that A∗v 6= 0. The possible so- lution(s) of (30) is the solution(s) of the following equation µbθ µ2k4 ( 1 − A∗v K )( µ2k4 + k5λ ∗ v µv + λ∗v ) −k6 = 0 (31) The direct resolution of equation (31) gives A∗v = K   µ2µbθk4 ( 1 − 1 N ) + αλ∗v µbθ(µ2k4 + k5λ ∗ v)   (32) where N = µbθ µvk6 and α = µbθk5 −µ2k4k6. Let us first compute the equilibrium without Disease, i.e. λ∗h = λ ∗ v = 0 or Eh = Ih = Ev = Iv = 0. From (32), we obtain A0v := K ( 1 − µvk6 µbθ ) := K ( 1 − 1 N ) (33) 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 P0. 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 P1 =( S0h,V 0 h , 0, 0, 0,A 0 v,S 0 v, 0, 0 ) , where S0h = Λh k1 , V 0h = ξΛh k1µh , A0v = K ( 1 − 1 N ) , S0v = θ µv A0v. Appendix B: PROOF OF PROPOSITION 3.2 We consider the Jacobian matrix associated to model (3) at the equilibrium TE. we have J(TE) =   −k1 0 0 0 0 0 ξ −µh 0 0 0 0 0 0 −k2 0 0 0 0 0 γh −k3 0 0 0 0 0 σ −µh 0 0 0 0 0 0 −k6 0 0 0 0 0 θ 0 0 0 0 0 0 0 0 0 0 0 0 0 − βhvηvS 0 h N0h + m − βhvS 0 h N0h + m 0 − πβhvηvV 0 h N0h + m − πβhvV 0 h N0h + m 0 βhvηvS0 N0h + m βhvS0 N0h + m 0 0 0 0 0 0 µb µb µb −µv 0 0 0 −k4 0 0 γv −µ2   , were S0 = S0h + πV 0 h . The eigenvalues of J(TE) are given by λ1 = λ2 = −µh, λ3 = −k1, λ4 = −k2, λ5 = −k3, and λ6, λ7, λ8, λ9 are eigenvalues of the submatrix Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 21 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... J̄ =   −k6 µb µb µb θ −µv 0 0 0 0 −k4 0 0 0 γv −µ2   . The characteristic polynomial of J̄ is given by P(λ) = λ4 +A1λ3 +A2λ2 +A3λ+A4 = 0 (34) where A1 = µv + µ2 + k4 + k6, A2 = k6µv (1 −N) + (k4 + µ2) (µv + k6) + µ2k4, A3 = k6µv (1 −N) (k4 + µ2) + µ2k4 (µv + k6) and A4 = µ2k4 (1 −N). 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 H1 = A1, H2 = ∣∣∣∣A1 1A3 A2 ∣∣∣∣, H3 = ∣∣∣∣∣∣ A1 1 0 A3 A2 A1 0 A4 A3 ∣∣∣∣∣∣, H4 = ∣∣∣∣∣∣∣∣ A1 1 0 0 A3 A2 A1 1 0 A4 A3 A2 0 0 0 A4 ∣∣∣∣∣∣∣∣ = A4H3. The Routh-Hurwitz criterion of stability of the trivial equilibrium TE is given by  H1 > 0 H2 > 0 H3 > 0 H4 > 0 ⇔   H1 > 0 H2 > 0 H3 > 0 A4 > 0 (35) We have H1 = A1 = µv + µ2 + k4 + k6 > 0, H2 = A1A2 −A3 = (k6 + k4 + µ2) µ 2 v + ( µ2k6 ( 1 − µbθ µ2k6 ) +k26 + 2k4k6 + µ2k6 + k 2 4 + 2µ2k4 + µ 2 2 ) µv + µ2k 2 6 ( 1 − µbθ µ2k6 ) + k4k 2 6 + (k4 + µ2) 2 k6 + µ2k 2 4 + µ 2 2k4, H3 = A1A2A3 −A21A4 −A 2 3 = (k4 + µ2) (µv + k6) × ( k6µv (1 −N) + µ2µv + µ2k6 + µ22 ) × ( k6µv (1 −N) + k4µv + k4k6 + k24 ) , Using inequality 1/µ2 ≤ 1/µ1 ≤ 1/µv, we obtain H2 > 0. H3 > 0 if N < 1; A4 > 0 if and only if N < 1. 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 R0 as the dominant eigenvalue of the next–generation matrix [26], [82]. Observe that model (3) has four infected populations, namely Eh, Ih, Ev, Iv. 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 F = 1 N0h + m ×   0 0 βhvηvS0 βhvS0 0 0 0 0 βvhηhS 0 v βvhS 0 v 0 0 0 0 0 0   , with N0h = Λh µh , V =  (µh + γh) 0 0 0 −γh (µh + δ + σ) 0 0 0 0 (µ1 + γv) 0 0 0 −γv µ2   , and the dominant eigenvalue of the next– generation matrix FV −1 is given by Eq. (10). The local stability of the disease–free equilib- rium P1 is a direct consequence of Theorem 2 of [82]. This ends the proof. Appendix C: PROOF OF PROPOSITION 3.3 Setting Y=X-TE with X = (Sh,Vh,Eh,Ih,Rh,Av,Sv,Ev,Iv) T , we can rewrite (3) in the following manner dY dt = B(Y )Y (36) Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 22 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... where B(Y ) =   −λh −k1 0 0 0 0 ξ −πλh −µh 0 0 0 λh πλh −k2 0 0 0 0 γh −k3 0 0 0 0 σ −µh 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 − βhvηvS 0 h Nh + m − βhvS 0 h Nh + m 0 0 − πβhvηvV 0 h Nh + m − πβhvV 0 h Nh + m 0 0 βhvηvS0 Nh + m βhvS0 Nh + m 0 0 0 0 0 0 0 0 −A66 µb µb µb θ −λv −µv 0 0 0 λv −k4 0 0 0 γv −µ2   , with S0 = (S0h + πV 0 h ), A66 =( k6 + µb Sv + Ev + Iv K ) . It is clear that Y = (0, 0, 0, 0, 0, 0, 0, 0, 0) is the only equilibrium. Then it suffices to consider the following Lyapunov function L(Y ) =< g,Y > were g = ( 1, 1, 1, 1, 1, 1, k6 θ , k6 θ , k6 θ , k6 θ ) . Straightforward computations lead that L̇(Y ) =< g,Ẏ >def= < g,B(Y )Y > = −µh(Y1 + Y2 + Y3 + Y4 + Y5) −δY4 + k6µv θ (N − 1)Y7 + k6µ1 θ ( µbθ k6µ1 − 1 ) Y8 −µb Y6 K (Y7 + Y8 + Y9) + k6µ2 θ ( µbθ k6µ2 − 1 ) Y9 Using the fact that 1/µ2 ≤ 1/µ1 ≤ 1/µv, we have µbθ k6µ1 − 1 ≤ 0 and µbθ k6µ2 − 1 ≤ 0, which implies that L̇(Y ) ≤ 0 if N ≤ 1. Moreover, the maximal invariant set contained in { L|L̇(Y ) = 0 } is {(0, 0, 0, 0, 0, 0, 0, 0, 0)}. Thus, from Lyapunov theory, we deduce that (0, 0, 0, 0, 0, 0, 0, 0, 0) and thus, TE := P0, is GAS if N ≤ 1. Appendix D: PROOF OF PROPOSITION 4.1. We compute now the endemic equilibrium, i.e. we are looking for an equilibrium such that λ∗h 6= 0 and λ∗v 6= 0. We assume that N > 1. From the sixth equation of (26), at equilibrium, we have S∗v + E ∗ v + I ∗ v = Kk6A ∗ v µb(K −A∗v) (37) From the last third equations of (26), at equilib- rium, we have µvS ∗ v + µ1E ∗ v + µ2I ∗ v = θA ∗ v (38) we will observe the following two cases. a) Absence of disease–induced death in vec- tor: The absence of disease–induced death in vector is traduce by the relation µv = µ1 = µ2, then equation (38) becomes S∗v + E ∗ v + I ∗ v = θ µv A∗v (39) Equalling Eqs. (37) and (39) gives like before A0v := K ( 1 − µvk6 µbθ ) = K ( 1 − 1 N ) . (40) Substituting A∗v by A 0 v in equation (29) gives S∗v = ( 1 − 1 N ) Kθ (µv + λ∗v) , E∗v = ( 1 − 1 N ) Kθλ∗v k4(µv + λ∗v) , I∗v = ( 1 − 1 N ) Kθγvλ ∗ v µvk4(µv + λ∗v) . (41) Replacing (41) in the expression of λ∗h gives λ ∗ h = βhv(ηvE ∗ v + I ∗ v ) (N∗h + m) = k10 λ∗v (µv + λ∗v) ×( βhvµhk2k3(k1 + λ ∗ h)(πλ ∗ h + µh) k2k3k7(k1 + λ ∗ h)(πλ ∗ h + µh) − δγhΛhλ ∗ h(k8 + πλ ∗ h) ) (42) where k7 = (Λh + mµh), k8 = πξ + µh, k9 = µ2ηv + γv = ηvµv + γv and k10 = k9Kθ µvk4 ( 1 − 1 N ) . Replacing (28) in the expression of λ∗v gives λ ∗ v =( βvhΛhµhk11λ ∗ h(k8 + πλ ∗ h) k2k3k7(k1 + λ ∗ h)(πλ ∗ h + µh) − δγhΛhλ ∗ h(k8 + πλ ∗ h) ) (43) Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 23 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... where k11 = k3ηh + γh. Substituting (43) in (42) gives the following equation f(λ ∗ h) := λ ∗ h [ B4(λ ∗ h) 4 + B3(λ ∗ h) 3 +B2(λ ∗ h) 2 + B1λ ∗ h + B0 ] = 0 (44) where B4 = π 2 [k7(µhk3 + γh(µh + σ)) + δγhmµh] ×{µv [k7(µhk3 + γh(µh + σ)) + δγhmµh] +βvhk11Λhµh} B3 = 2πX{k2k3k7µh(1 + π) + δγhΛhµh +πξX}µv + βvhπΛhµhk11 {πk2k3Y +k7 [k8(k2k3 − 2δγh) + µhk2k3]} B2 = µv [k1k2k3k7π −δΛhγhπξ + Xµh] 2 + 2k1k2k3k7πµhµvX + βvhΛhµ 2 hπk2k3k11 {πk1k7 −βhvk10 [π(k8 + k1) + µh]} + βvhk8k11Λhµh [k2k3k7πµh + k8X] B1 = 2k1k2k3k7µhµv [k8X + πµhk2k3k7] + k2k3k11βvhΛhµ 2 h{k1k7k8 −βhvk10(µhk8 + k1π(k8 + µh))} with X = k2k3k7−δγhΛh, Y = k1k7−βhvµhk10; and B0 = µ 2 hµvk 2 1k 2 2k 2 3k 2 7 ( 1 −R20 ) We consider λ∗h 6= 0, otherwise we recover DFE. The positive endemic equilibria P∗ = (S∗h,V ∗ h ,E ∗ h,I ∗ h,R ∗ h,A ∗ v,S ∗ v,E ∗ v,I ∗ v ) are obtained by solving Eq. (44) for λ∗h. The coefficient B4 is always positive and coefficient B0 is negative (resp. positive) whenever R0 > 1 (resp. R0 < 1). The number of possible nonnegative real roots of the polynomial of Eq. (44) depends on the signs of B3, B2 and B1. The various possibilities for the roots of f(λ∗h) are tabulated in Table XI and XII. From tables XI and XII , we deduce the fol- lowing result which gives various possibilities of nonnegative solutions of (44). Lemma A.1: Assume that N > 1 and µv = µ1 = µ2. Then, the arboviral-disease model (3) TABLE XI TOTAL NUMBER OF POSSIBLE REAL ROOTS OF (44) WHEN R0 > 1. Number of Cases B0 B1 B2 B3 B4 sign changes – + + + + 1 1 – – + + + 1 – – – + + 1 – – – – + 1 – + + – + 3 – + – + + 3 2 – + – – + 3 – – + – + 3 TABLE XII TOTAL NUMBER OF POSSIBLE REAL ROOTS OF (44) WHEN R0 < 1. Number of Cases B0 B1 B2 B3 B4 sign changes 1 + + + + + 0 + + + – + 2 + + – + + 2 + + – – + 2 2 + – + + + 2 + – – + + 2 + – – – + 2 3 + – + – + 4 1. has a unique endemic equilibrium when Case 1 of Table XI is satisfied and whenever R0 > 1. 2. could have more than one endemic equilib- rium when Case 2 of Table XI is satisfied whenever R0 > 1. 3. could have more than one endemic equilib- rium when Case 2, 3 of Table XII are satisfied and whenever R0 < 1. 4. has no endemic equilibrium when Case 1 of Table XII is satisfied and whenever R0 < 1. Case 3 of Lemma A.1 suggests that co-existence of the disease–free equilibrium and the endemic equilibrium for the arboviral-disease model (3) is possible, and hence the potential occurrence of the backward bifurcation phenomenon when R0 < 1. Also, case 2 of Lemma A.1 suggests the possibility of a pithcfork (Forward) bifurcation when R0 = 1. Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 24 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... b) Presence of disease–induced death in vec- tor:: Here, we will consider µv < µ1 < µ2 with µv 6= µ2. Equation (27) becomes λ∗h = βhvµhµvk2k3k10(k1 + λ ∗ h)(µh + πλ ∗ h) k2k3k7(k1 + λ ∗ h)(µh + πλ ∗ h) − δγhΛhλ ∗ h(k8 + πλ ∗ h) × ( µ2µbθk4N1 + αλ∗v µbθ(µ2k4 + k5λ∗v) )( λv µv + λv ) (45) with k12 = µvk10 K , N1 = ( 1 − 1 N ) and λ∗v = βvhµhΛhk11λ ∗ h(k8 + πλ ∗ h) k2k3k7(k1 + λ ∗ h)(µh + πλ ∗ h) − δγhΛhλ ∗ h(k8 + πλ ∗ h) (46) Substituting (46) in (45) gives the following equa- tion λ∗h 6∑ i=0 Ci(λ ∗ h) i = 0 (47) where C0 = k31k 3 2k 3 3k4k 3 7θµ2µvµbµ 3 h ( R20 − 1 ) and C6 = −µbπ3θX (µ2k4X + βvhk5k11Λhµh) ×(µvX + βvhk11Λhµh) , with X = (k2k3k7 −δΛhγh) > 0. The others coefficients C5, C4, C3, C2, and C1 are obtained after computations on Maxima software. We also obtain the following result which gives various possibilities of solutions of Eq. (47). Lemma A.2: Assume that N > 1. Then, the arboviral-disease model (3) 1. could have a unique endemic equilibrium whenever R0 > 1. 2. could have more than one endemic equilib- rium whenever R0 > 1. 3. haven’t endemic equilibrium whenever R0 < 1. 4. could have one or more than one endemic equilibrium whenever R0 < 1. Case 4 of Lemma A.2 suggests that co-existence of the disease–free equilibrium and endemic equi- librium for the arboviral-disease model (3) is possible, and hence the potential occurrence of a backward bifurcation phenomenon when R0 < 1. Also, case 2 of Lemma A.2 suggests the possibility of a pithcfork (Forward) bifurcation when R0 = 1. Appendix E: COMPUTATION OF a? OF THEOREM 4.1. a? = 1 2 v3 9∑ i,j=1 wiwj ∂2f3(0, 0) ∂xi∂xj + 1 2 v8 9∑ i,j=1 wiwj ∂2f8(0, 0) ∂xi∂xj . (48) Let a?3 = ∑9 i,j=1 wiwj ∂2f3(0, 0) ∂xi∂xj and a?8 = ∑9 i,j=1 wiwj ∂2f8(0, 0) ∂xi∂xj . After few computa- tions, we obtain a ? 3 = β∗hvµh(�ξΛh + mµh) k1(Λh + mµh)2 w1 (ηvw8 + w9) + β∗hvπµh(�ξΛh + mµh) k1(Λh + mµh)2 w2 (ηvw8 + w9) − β∗hvµhΛh(µh + πξ) k1(Λh + mµh)2 w3 (w8ηv + w9) − β∗hvµhΛh(µh + πξ) k1(Λh + mµh)2 w4 (w8ηv + w9) − β∗hvµhΛh(µh + πξ) k1(Λh + mµh)2 w5 (w8ηv + w9) + β∗hvηvµh k1(Λh + mµh)2 [(�ξΛh + mµh) (w1 + π w2) −Λh(µh + πξ)w8 (w3 + w4 + w5)] + β∗hvµh k1(Λh + mµh)2 [(�ξΛh + mµh) (w1 + π w2) −Λh(µh + πξ)w9 (w3 + w4 + w5)] = β∗hvµh k1(Λh + mµh)2 {(�ξΛh + mµh) (w1 + π w2) −Λh(µh + πξ) (w3 + w4 + w5)}(w8ηv + w9) + β∗hvµh k1(Λh + mµh)2 {(�ξΛh + mµh) (w1 + π w2) (ηv + 1) −Λh(µh + πξ) (w3 + w4 + w5) (ηvw8 + w9)} = β∗hvµh k1(Λh + mµh)2 × {(�ξΛh + mµh) (w1 + π w2) (w8ηv + w9 + ηv + 1) −2Λh(µh + πξ) (w3 + w4 + w5) (w8ηv + w9)} , Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 25 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... a ? 8 = w3 9∑ j=1 wj ∂2f8 ∂x3∂xj (x0, 0) + w4 9∑ j=1 wj ∂2f8 ∂x4∂xj (x0, 0) + w7 9∑ j=1 wj ∂2f8 ∂x5∂xj (x0, 0) + w8 9∑ j=1 wj ∂2f8 ∂x8∂xj (x0, 0) + w9 9∑ j=1 wj ∂2f8 ∂x9∂xj (x0, 0) = − βvhµ 2 hKθ ( 1 − 1 N ) µv(Λh + µhm)2 [ ηh ( w3 + 5∑ i=1 wi ) + w4 ] w3 − βvhµ 2 hKθ ( 1 − 1 N ) µv(Λh + µhm)2 [( w4 + 5∑ i=1 wi ) + ηhw3 ] w4 + βvhµh (Λh + µhm) (ηhw3 + w4) w7 = − βvhµ 2 hKθ ( 1 − 1 N ) µv(Λh + µhm)2 [(ηhw3 + w4) (w1 + w2 + w5) +2(ηh + 1)w3w4 + 2 ( ηhw 2 3 + w 2 4 )] + βvhµh (Λh + µhm) (ηhw3 + w4) w7 Using above results, Eq. (48) becomes a? = φ1 −φ2 where φ1 = 1 2 v3 { β∗ hv µh k1(Λh + mµh) 2 [(�ξΛh + mµh) (w1 + π w2) ×(w8ηv + w9 + ηv + 1) −2Λh(µh + πξ) (w3 + w4 + w5) (w8ηv + w9)]} − 1 2 v8 βvhµ 2 h Kθ µv(Λh + µhm) 2 ( 1 − 1 N ) w5 (ηhw3 + w4) + 1 2 v8 βvhµh (Λh + µhm) (ηhw3 + w4) w7 − 1 2 v8 βvhµ 2 h Kθ ( 1 − 1 N ) µv(Λh + µhm) 2 [2(ηh + 1)w3w4 +2 ( ηhw 2 3 + w 2 4 )] < 0 and φ2 = 1 2 v8 βvhµ 2 hKθ µv(Λh + µhm) 2 ( 1 − 1 N ) × [(ηhw3 + w4) (w1 + w2)] < 0 Appendix F: PROOF OF LEMMA 7.2 The Kamgang-Sallet approach used for (22) ensures that the trivial equilibrium(TE := P0) and the disease–free equilibrium (DFE := P1) are the fixed point of (22). Indeed, rewriting (22) gives  Sk+1h = φ(h)Λh + (1 −φ(h)k1)Skh 1 + φ(h)λkh V k+1h = φ(h)ξSkh + (1 −φ(h)µh)V k h 1 + φ(h)πλkh Ek+1h = (1 −φ(h)k2)E k h + φ(h)λ k h ×(Sk+1h + πV k+1 h ) Ik+1h = φ(h)γhE k h + (1 −φ(h)k3)I k h Rk+1h = φ(h)σI k h + (1 −φ(h)µh)R k h Ak+1v = [ 1−φ(h) ( k6+µb Skv+E k v+I k v K )] Akv + φ(h)µb(S k v + E k v + I k v ) Sk+1v = φ(h)θAkv + (1 −φ(h)µv)Skv 1 + φ(h)λkv Ek+1v = (1 −φ(h)k4)Ekv + φ(h)λkvSk+1v Ik+1v = φ(h)γvE k v + (1 −φ(h)µ2)Ikv (49) If X∗ = (Sh,V ∗ h ,E ∗ h,I ∗ h,R ∗ h,A ∗ v,S ∗ v,E ∗ v,I ∗ v ) T is an equilibrium of the discrete system (49), then we have after few simplifications  Λh −λ∗hS ∗ h −k1S ∗ h = 0 ξS∗h −πλ ∗ hV ∗ h −µhV ∗ h = 0 λ∗h(S ∗ h + πV ∗ h ) −k2E ∗ h = 0 γhE ∗ h −k3I ∗ h = 0 σI∗h −µhR ∗ h = 0 µb(S ∗ v + E ∗ v + I ∗ v ) ( 1 − A∗v K ) −k6A∗v = 0 θA∗v −λ∗vS∗v + µvS∗v = 0 k4E ∗ v −λ∗vS∗v = 0 γvE ∗ v −µ2I∗v = 0 (50) which is equivalent to{ A1(X∗)(X∗S −XDFE) + A12(X ∗)X∗I = 0 A2(X∗)X∗I = 0 (51) where A1, A12 and A2 are given at Equation (11). Appendix G: PROOF OF LEMMA 7.3 The Jacobian matrix associated with the right- hand side of the numerical scheme (22) at the Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 26 of 30 http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... tivial equilibrium TE := P0 is given by JTE = (Jij)1≤i,j≤9 with J1,1 = 1 −k1φ(h); J1,8 = − φ(h)βhvηvΛhµh k1(Λh + µhm) ; J1,9 = − φ(h)βhvΛhµh k1(Λh + µhm) ; J2,1 = φ(h)ξ; J2,2 = 1 −µhφ(h); J2,8 = − φ(h)πβhvηvξΛh k1(Λh + µhm) , J2,9 = − φ(h)πβhvξΛh k1(Λh + µhm) , J3,3 = 1 −k2φ(h); J3,8 = φ(h)βhvηvΛh(µh + πξ) k1(Λh + µhm) ; J3,9 = φ(h)βhvΛh(µh + πξ) k1(Λh + µhm) ; J4,3 = φ(h)γh, J4,4 = 1 −k3φ(h); J5,4 = φ(h)σ; J5,5 = 1 −µhφ(h); J6,6 = 1 −φ(h)k6; J6,7 = J6,8 = J6,9 = φ(h)µb; J7,6 = φ(h)θ, J7,7 = 1 −φ(h)µv, J8,8 = 1 −φ(h)k4; J9,8 = φ(h)γv; J9,9 = 1 −µ2φ(h) The eigenvalues of JTE are given by λ1 = λ2 = 1 −µhφ(h), λ3 = 1 −k1φ(h), λ4 = 1 −k2φ(h), λ5 = 1−k3φ(h), and λ6, λ7, λ8, λ9 are eigenval- ues of the submatrix J̄ =   J6,6 φ(h)µb φ(h)µb φ(h)µb J7,6 J7,7 0 0 0 0 J8,8 0 0 0 J9,8 J9,9   Since φ(h) > 0, it is clear that |λi| < 1, for i = 1, 2, . . . , 5. We need also to show that the characteristic polynomial associated with J̄ is Schur polynomials, i.e. polynomials such that all roots λi verify |λi| < 1. The characteristic polynomial associated with J̄ is given by P(λ) = (λ + µ2φ(h) − 1) (λ + k4φ(h) − 1) H(λ) where H(λ) = λ 2 + (φ(h)(µv + k6) − 2) λ + 1 + φ(h) 2 (k6µv −µbθ) −φ(h)(µv + k6) The roots of P(λ) are λ6 = 1 − µ2φ(h), λ7 = 1 − k4φ(h) and the others roots are the roots of H(λ). Note that |λ6| < 1 and |λ7| < 1. Now, we need to show that H(λ) is a Schur polynomial. To this aim, let q1 = (φ(h)(µv + k6) − 2) and q2 = 1 + φ(h)2(k6µv −µbθ)−φ(h)(µv + k6). Using Lemma 11 in [29], we just show that the following conditions hold: 1 +q1 +q2 > 0, 1−q1 +q2 > 0, 1−q2 > 0 (52) We compute 1 + q1 + q2 = φ(h)2k6µv(1 −N), 1 − q1 + q2 = 2 [(1 −φ(h)µv) + (1 −φ(h)k6)] + φ(h)2k6µv(1 −N) and 1 −q2 = φ(h) [µv + k6(1 −φ(h)µv) + φ(h)µbθ] If φ(h) has choosen as equation (25), then conditions (52) hold whenever N ≤ 1. Thus, H(λ) is a Schur polynomial. This ends the proof. REFERENCES [1] Ahmed Abdelrazec, Suzanne Lenhart, Huaiping Zhu, Transmission dynamics of West Nile virus in mosquitoes and corvids and non-corvids, J. Math. Biol. (2014) 68:1553–1582. DOI10.1007/s00285-013-0677-3 [2] Dipo Aldila, Thomas Gtz, Edy Soewono, An optimal con- trol problem arising from a dengue disease transmission model, Mathematical Biosciences 242 (2013) 9–16. [3] R. Anguelov, J. M. S. Lubuma, M. Shillor, Topologi- cal dynamic consistency of nonstandard finite difference schemes for dynamical systems, J. Difference Eq. Appl., 17, 1769–1791 (2011) [4] R. Anguelov, Y. Dumont, J. M. S. Lubuma, On nonstan- dard finite difference schemes in biosciences, Proceeding of the fourth international conference on application of mathematics in technical and natural sciences. American Institute of Physics conference, Proceedings AIP, 1487, 212–223 (2012) [5] R. Anguelov, Y. Dumont, J.M.-S. Lubuma, and M. Shillor, Dynamically consistent nonstandard finite dif- ference schemes for epidemiological Models. Journal of Computational and Applied Mathematics, 255, pp.161- 182, 2014. [6] Arino J., McCluskey C.C., van den Driessche P., Global results for an epidemic model with vaccination that exhibits backward bifurcation, SIAM Journal on Applied Mathematics 64 (2003) 260–276. [7] Bayer 360◦ Vector Control. Available at http://www.vectorcontrol.bayer.com/bayer/ cropscience/bes vectorcontrol.nsf/id/EN mosquitoes (Retrieved January 2014) [8] A. Berman, R.J. Plemmons, Nonnegative matrices in the mathematical sciences, SIAM., 1994. [9] N. P. Bhatia, G. P. Szego, Stability Theory of Dynamical Systems, Springer-Verlag, 1970 Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 27 of 30 DOI 10.1007/s00285-013-0677-3 http://www.vectorcontrol.bayer.com/bayer/ cropscience/bes_vectorcontrol.nsf/id/EN_mosquitoes http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... [10] Jr J.E. Blaney, N.S. Sathe, C.T. Hanson, C.Y. Firestone, B.R. Murphy, S.S. Whitehead, Vaccine candidates for dengue virus type 1 (Den 1) generated by replacement of the structural genes of rDen 4 and rDen4D30 with those of Den 1, Virology Journal 4 (23) (2007) 1–11. [11] Jr J.E. Blaney, J.M. Matro, B.R. Murphy, S.S. White- head, Recombinant, live-attenuated tetravalent dengue virus vaccine formulations induce a balanced, broad, and protective neutralizing anitibody response against each of the four serotypes in rhesus monkeys, Journal of Virology 79 (9) (2005) 5516–5528. [12] S. Blower, H. Dowlatabladi, Sensitivity and uncertainty analysis of complex models of disease transmission: an HIV model, as an example. Int. Stat. Rev., 62, 229–243 (1994) [13] Boloye Gomero, Latin Hypercube Sampling and Partial Rank Correlation Coefficient Analysis Applied to an Optimal Control Problem, Master Thesis, University of Tennessee, Knoxville, 2012. [14] Brauer F., Backward bifurcations in simple vaccination models, J. Math. Anal. Appl., 298, 418-431 (2004) [15] Buonomo B., D. Lacitignola, A note on the direction of the transcritical bifurcation in epidemic models. Nonlin- ear Analysis: Modelling and Control, 2011, Vol. 16, No. 1, 30–46. [16] Buonomo B., On the backward bifurcation of a vaccina- tion model with nonlinear incidence. Nonlin. Anal. Mod. Control, 20, 38–55 (2015). [17] J. R. Cannon, D. J. Galiffa, An epidemiology model suggested by yellow fever, Math. Methods Appl. Sci., 35, 196–206 (2012) [18] V. Capasso, Mathematical Structures of Epidemics Sys- tems. Lecture Notes in Biomathematics, 97. Springer- Verlag, Berlin, 2008 [19] C. Castillo–Chavez, B. Song, Dynamical models of tuberculosis and their applications, Math. Biosci. Eng, 1, 361–404 (2004) [20] A. Chippaux, Généralités sur arbovirus et arboviroses, overview of arbovirus and arbovirosis, Med. Maladies Infect., 33, 377–384 (2003) [21] N. Chitnis, D. Hardy, T. Smith, A Periodically–Forced Mathematical Model for the Seasonal Dynamics of Malaria in Mosquitoes. Bull. Math. Biol. 74, 1098–1124 (2012) [22] N. Chitnis, J. M. Hyman, J. M. Cushing, Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model, Bull. Math. Biol., 70, 1272–1296 (2008) [23] F.A.B. Coutinho, M.N. Burattini, L.F. Lopez, E. Massad, Threshold conditions for a non-autonomous epidemic system describing the population dynamics of dengue, Bulletin of Mathematical Biology 68 (2006) 2263–2282. [24] G. Cruz-Pacheco, L. Esteva, C. Vargas, Seasonality and outbreaks in West Nile Virus infection, Bull. Math. Biol., 71, 1378–1393 (2009) [25] Cushing J. M. and Yicang Z., The net reproductive value and stability in matrix population models, Nat. Resour. Model., 8 (1994), pp. 297–333. [26] O. Diekmann, J. A. P. Heesterbeek, Mathematical Epi- demiology of Infectious Diseases. Model building, anal- ysis and interpretation. John Wiley & Sons, Chichester, 2000 [27] K. Dietz, Transmission and control of arbovirus dis- eases, Epidemiology (ed. D. Ludwig, K. L. Cooke), pp. 104–121, SIAM, Philadelphia, 1975 [28] M. Dubrulle, L. Mousson, S. Moutailler, M. Vazeille and A.-B. Failloux, Chikungunya virus and aedes mosquitoes: Saliva is infectious as soon as two days after oral infection, PLoS One, 4 (2009). [29] Y. Dumont, F. Chiroleu, C. Domerg, On a temporal model for the Chikungunya disease: Modeling, theory and numerics, Mathematical Biosciences 213 (2008) 80– 91. [30] Y. Dumont, F. Chiroleu, Vector control for the chikun- gunya disease, Math. Biosci. Eng., 7, 313–345 (2010) [31] J. Dushoff, W. Huang, C. Castillo–Chavez, Backward bifurcations and catastrophe in simple models of fatal diseases, J. Math. Biol., 36, 227–248 (1998) [32] M. Dubrulle, L. Mousson, S. Moutailler, M. Vazeille, A.-B. Failloux, Chikungunya virus and aedes mosquitoes: Saliva is infectious as soon as two days after oral infection, PLoS One, 4 (2009). [33] K.H. Eckels, R. Putnak, Formalin-inactivated whole virus and recombinant subunit flavivirus vaccines, Ad- vances in Virus Research 61 (2003) 395–418. [34] Y. Eshita, T. Takasaki, I. Takashima, N. Komalamisra, H. Ushijima and I. Kurane, Vector competence of Japanese mosquitoes for dengue and West Nile viruses, Pesticide Chemistry (2007) doi:10.1002/9783527611249. ch23. [35] L. Esteva, C. Vargas, Analysis of a dengue disease transmission model, Math. Biosci., 150, 131–151 (1998) [36] L. Esteva, C. Vargas, A model for dengue disease with variable human population, J. Math. Biol., 38, 220–240 (1999) [37] F. Forouzannia, A.B. Gumel, Mathematical Analysis of an Age-structured Model for Malaria Transmission Dynamics, Mathematical Biosciences (2013), doi:http: //dx.doi.org/10.1016/j.mbs.2013.10.011 [38] Z. Feng, V. Velasco–Hernandez, Competitive exclusion in a vector–host model for the dengue fever, J. Math. Biol., 35, 523–544 (1997) [39] Figaro, Un vaccin contre la dengue disponible dès la mi-2015. http://www.lefigaro.fr/societes/2014/11/04/20005 -20141104ARTFIG00301-un-vaccin-contre-la-dengue-disponible -des-la-mi-2015.php (Accessed April 2015) [40] S.M. Garba, A.B. Gumel, M.R. Abu Bakar, Backward bifurcations in dengue transmission dynamics, Math. Biosci., 215, 11–25 (2008) [41] B. S. Goh, Global stability in a class of prey-predator models. Bull. Math. Biol. 40, 525-533 (1978) Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 28 of 30 doi:10.1002/9783527611249.ch23 doi:10.1002/9783527611249.ch23 http://dx.doi.org/10.1016/j.mbs.2013.10.011 http://dx.doi.org/10.1016/j.mbs.2013.10.011 http://www.lefigaro.fr/societes/2014/11/04/20005 -20141104ARTFIG00301-un-vaccin-contre-la-dengue-disponible -des-la-mi-2015.php http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... [42] D. J. Gubler, Human arbovirus infections worldwide, Ann. N. Y. Acad. Sci., 951, 13–24 (2001) [43] J. Guckenheimer, P. Holmes, Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields, Springer-Verlag, Berlin, 1983. [44] L. Guillaumot, Les moustiques et la dengue (in French). Institut Pasteur de Nouvelle–Caledonie, 2005 [45] J. K. Hale, Ordinary Differential Equations. John Wiley and Sons, New York. 1969 [46] Institut Pasteur, Chikungunya. http://www.pasteur.fr/fr/institut-pasteur/presse/ fiches-info/chikungunya#Traitement (Accessed August 2014) [47] J.A. Jacquez, C.P. Simon,Qualitative theory of compart- mental systems. SIAM Rev. 35 (1993), 43–79. [48] J. C. Kamgang, G. Sallet, Computation of threshold conditions for epidemiological models and global stabil- ity of the disease-free equilibrium (DFE), Mathematical Biosciences 213 (2008) 1–12. [49] N. Karabatsos, International Catalogue of Arboviruses, including certain other viruses of vertebrates. American Society of Tropical Medicine and Hygiene. San Antonio, TX. 1985, 2001 update [50] P. Koraka, S. Benton, G. van Amerongen, K.J. Stitelaar, A.D.M.E. Osterhaus, Efficacy of a live attenuated tetrava- lent candidate dengue vaccine in naive and previously infected cynomolgus macaques, Vaccine 25 (2007) 5409– 5416. [51] J. P. LaSalle, Stability theory for ordinary differential equations, J. Differ. Equ., 57–65 (1968) [52] J. P. LaSalle, The stability of dynamical systems, Society for Industrial and Applied Mathematics, Philadelphia, Pa., 1976. [53] Le Figaro, Chikungunya: Des vaccins bientôt testés chez l’homme, http://www.lefigaro.fr/sciences (Accessed August 2014) [54] Le Monde, Test prometteur pour un nouveau vaccin contre le chikungunya, http://www.lemonde.fr/sante/article (Accessed August 2014) [55] Lee-Jah Chang et al., Safety and tolerability of chikungunya virus-like particle vaccine in healthy adults: a phase 1 dose–escalation trial. http://dx.doi.org/10.1016/S0140-6736(14)61185-5. [56] N. A. Maidana, H. M. Yang, Dynamic of West Nile Virus transmission considering several coexisting avian populations, Math. Comput. Modelling, 53, 1247–1260 (2011) [57] MATLAB c©. Matlab release 12. The mathworks Inc., Natich, MA, 2000. [58] Michel Gosse, Faq Maxima, Version 0.93 du 13 juillet 2010. [59] R. E. Mickens, Nonstandard Finite Difference Models of Differential Equations. World Scientific, Singapore, 1994 [60] D. Moulay, M. A. Aziz–Alaoui, M. Cadivel, The Chikungunya disease: Modeling, vector and transmission global dynamics, Math. Biosci., 229,50–63 (2011) [61] Moulay D., Aziz–Alaoui M. A., Hee-Dae Kwon, Op- timal control of Chikungunya disease: larvae reduc- tion,treatment and prevention, Mathemtical Biosciences and Engineering, volume 9, Number 2, April 2012. [62] S. Naowarat, W. Tawarat, I. Ming Tang, Control of the transmission of chikungunya fever epidemic through the use of adulticide, Am. J. Appl. Sci., 8, 558–565 (2011) [63] S. Naowarat, P. Thongjaem, I. Ming Tang, Effect of mosquito repellent on the transmission model of chikun- gunya fever, Am. J. Appl. Sci., 9, 563–569 (2012) [64] P. Poletti, G. Messeri, M. Ajelli, R. Vallorani, C. Rizzo, S. Merler, Transmission potential of chikungunya virus and control measures: the case of Italy, PLoS One, 6, e18860 (2011) [65] Richard Taylor. Interpretation of the correlation coef- ficient: A basic review, Journal of Diagnostic Medical Sonography, 6(1):35–39, 1990. [66] R Development Core Team: a language and environ- ment for statistical computing. R Foundation for Sta- tistical Computing, Vienna, http://www.r-395project.org/ foundation [67] H.S. Rodrigues, M.T.T. Monteiro, D.F.M. Torres, A. Zinober, Dengue disease, basic reproduction number and control, Int. J. Comput. Math. 89 (3) (2012) 334. [68] H. S. Rodrigues, M. Teresa T. Monteiro, Delfim F.M. Torres, Vaccination models and optimal control strategies to dengue, Mathematical Biosciences 247 (2014) 1–12. [69] H.S. Rodrigues, M.T.T. Monteiro, D.F.M. Torres, In- secticide control in a dengue epidemics model, in: T.E. Simos, et al. (Eds.), Numerical analysis and applied mathematics. International conference on numerical anal- ysis and applied mathematics, Rhodes, Greece. American Institute of Physics Conf. Proc., no. 1281 in American Institute of Physics Conf. Proc., 2010, pp. 979–982. [70] M. Safan, M. Kretzschmar, K. P. Hadeler, Vaccination based control of infections in SIRS models with reinfec- tion: special reference to pertussis, J. Math. Biol., 67, 1083–1110 (2013) [71] M. Sanchez, S. Blower, Uncertainty and sensitivity analysis of the basic reproductive rate. American Journal of Epidemiology, 145: 1127 - 1137 (1997) [72] SANOFI PASTEUR, Dengue vaccine, a priority for global health, (2013) [73] J.F. Saluzzo, Empirically derived live attenuated vac- cines against dengue and Japanese encephalitis, Ad- vances in Virus Research 61 (2003), 419–443. [74] A. J. Saul, P. M. Graves, B. H. Kay, A cyclical feeding model for pathogen transmission and its application to determine vectorial capacity from vector infection. J. Appl. Ecology, 27, 123–133 (1990) [75] O. Sharomi, C.N. Podder, A.B. Gumel, E.H. Elbasha, J. Watmough, Role of incidence function in vaccine-induced backward bifurcation in some HIV models, Mathematical Biosciences 210 (2007) 436–463. [76] SCILAB. Open source software for numerical compu- tation. http://www.scilab.org [77] Z. Shuai, P. van den Driessche, Global dynamics of a Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 29 of 30 http://www.pasteur.fr/fr/institut-pasteur/presse/ fiches-info/chikungunya#Traitement http://www.lefigaro.fr/sciences http://www.lemonde.fr/sante/article http://dx.doi.org/10.1016/S0140-6736(14)61185-5 http://www.r- 395 project.org/foundation http://www.r- 395 project.org/foundation http://www.scilab.org http://dx.doi.org/10.11145/j.biomath.2015.07.241 H. Abboubakar et al., Modeling the Dynamics of Arboviral Diseases ... disease model including latency with distributed delays. Can. Appl. Math. Q., 19, 235–253 (2012). [78] M. Stein, Large Sample Properties of Simulations Using Latin Hypercube Sampling. Technometrics, 29, 143–151 (1987) [79] A.K. Supriatna, E. Soewono, S.A. van Gils, A two- age-classes dengue transmission model, Mathematical Biosciences 216 (2008) 114–121. [80] J. J. Tewa, J. L. Dimi, S. Bowong, Lyapunov functions for a dengue disease transmission model, Chaos Solit. Fract., 39, 936–941 (2009) [81] Valaire Yatat, Yves Dumont, Jean Jules Tewa, Pierre Courteron, Samuel Bowong, Mathematical Analysis of a Size Tree-Grass Competition Model for savanna Ecosys- tems, BIOMATH 3 (2014),1404214, http://dx.doi.org/10. 11145/j.biomath.2014.04.212 [82] P. van den Driessche,J. Watmough, Reproduction num- bers and the sub-threshold endemic equilibria for com- partmental models of disease transmission, Math. Biosci., 180, 29–48 (2002) [83] C. Vargas, L. Esteva, G. Cruz–Pacheco, Mathematical modelling of arbovirus diseases, 7th International Confer- ence on Electrical Engineering, Computing Science and Automatic Control, CCE 2010 [84] World Health Organization. Dengue and severe dengue. Fact sheet n.117. Updated September 2013. Available at www.who.int/mediacentre/factsheets/fs117/en (Retrieved January 2014) [85] World Health Organization, Immunological correlates of protection induced by dengue vaccines, Vaccine. 25 (2007) 4130–4139. [86] World Health Organization, Dengue and dengue haemorhagic fever, factsheet No.117, 2009. [87] A. Wilder–Smith, W. Foo, A. Earnest, S. Sremulanathan, N.I. Paton, Seroepidemiology of dengue in the adult pop- ulation of Singapore, Tropical Medicine and International Health 9 (2) (2004) 305–308. [88] L. Wu, B. Song, W. Du, J. Lou, Mathematical modelling and control of echinococcus in Qinghai province, China, Math. Biosci. Eng., 10, 425–444 (2013) [89] wxMaxima 11.08.0 c© 2004–2011 Andrej Vodopivec, http://andrejv.github.com/wxmaxima/ [90] A. Yebakima et al., Genetic heterogeneity of the dengue vector Aedes aegypti in Martinique, Tropical Medicine and International Health 9 (5) (2004) 582–587. Biomath 4 (2015), 1507241, http://dx.doi.org/10.11145/j.biomath.2015.07.241 Page 30 of 30 http://dx.doi.org/10.11145/j.biomath.2014.04.212 http://dx.doi.org/10.11145/j.biomath.2014.04.212 www.who.int/mediacentre/factsheets/fs117/en http://andrejv.github.com/wxmaxima/ http://dx.doi.org/10.11145/j.biomath.2015.07.241 Introduction Model formulation, Invariant region. The disease–free equilibria and stability analysis Local stability analysis Global stability analysis Global asymptotic stability of the trivial equilibrium TE:=P0 Global asymptotic stability of the disease–free equilibrium The endemic equilibria and bifurcation analysis Existence of endemic equilibria Bifurcation analysis Threshold analysis and vaccine impact Sensitivity analysis Mean values of parameters and initial values of variables Uncertainty and sensitivity analysis Sensitivity analysis of R0 Sensitivity analysis of Infected states of model (3) Numerical simulation A nonstandard scheme for the model (3) Simulation Results Conclusions Appendix References