BIFURCATION ANALYSIS IN A SEROGROUP A MENINGOCOCCAL INFECTION MODEL

1Department of Mathematics, Faculty of Science and Technic, University Cheikh Anta Diop, Dakar, Senegal 2Department of Computer Engineering, National Advanced School of Engineering, University Cheikh Anta Diop, Dakar, Senegal 3Laboratory of Mathematics, Department of Mathematics and Computer Science, Faculty of Science, University of Douala, PO Box 24157 Douala, Cameroon 4IRD, Sorbonne Université, UMMISCO, F-93143, Bondy, France


INTRODUCTION
Neisseria meningitidis is a major cause of invasive bacterial infections globally [1,2]. A notable feature of the meningococcus is its fluid epidemiology Various different agents can cause meningitis. The most common ones are either bacterial or viral organisms. The presence of these organisms in the subarachnoid space (the space between the pia mater and the arachnoid mater occupied by the cerebrospinal fluid) causes an inflammatory reaction which involves all the structures that lie within or adjacent to these areas. Almost any bacteria gaining entrance to the meninges may produce meningitis. Bacterial meningitis is sometimes due to infection in another part of the body (lungs, ears, nose, throat, sinuses) that spreads into the meninges.
About 10% of the population carry the germs for days, weeks, or months without becoming ill. Meningococcal carriage is most prevalent in teenagers and young adults, where up to 60% of the carriers are often recorded. Instead, meningococcal carriage seems particularly low in young children and babies. Also, invasive meningitis is most common in children between ages of 1 month and 5 years. It is much less common in adults, unless they have a special risk factor, such as impaired immune system, chronic ear and nose infections, pneumococcal pneumonia, sickle cell disease, and others. If meningitis is not early treated, can lead to swelling of the fluid surrounding the brain and spinal column as well as severe and permanent disabilities, and even death (it is fatal in the 50%-80% of untreated cases). Even with an early diagnosis and an adequate treatment, 5% to 10% patients die usually within 24 to 48 hours after the onset of symptoms (see, for example, [2,3]).
Neisseria meningitidis serogroups B and C are the most common cause of disease in Europe, the Americas, Australia, and New Zealand. Cases tend to occur more frequently in winter and spring. Neisseria meningitidis serogroup A is the main cause of disease in Africa and Asia.
Worldwide, the highest rates of disease occur in the "meningitis belt" of sub-Saharan Africa [4][5][6][7][8][9][10][11]. This extends across the dry, savannah regions from Senegal in the west, across to Ethiopia in the east. During epidemics, this region can have an annual incidence rate of 1.000 cases per 100.000 population [2]. The largest recorded outbreak of meningococcal disease in history occurred in Africa in 1996 when over 20,000 deaths were reported to the World Health Organization (WHO) [4]. In the 2009 epidemic season 14 African countries reported 78,416 suspected cases and 4.053 deaths. This is the largest reported number since the 1996 epidemic [4]. More than 85% of these cases were reported in one epidemic foci, encompassing northern Nigeria and Niger, and were characterized by the predominance of MenA [4]. Once a successfully diagnosis is made, the antibiotic treatment must be started as soon as possible. Empirical therapy includes ceftriaxone or cefotaxime, and vancomycin for Streptococcus pneumoniae. There is a vaccine against meningococcal disease which is 85%-100% effective in preventing four kinds of bacteria (the serogroups A, C, Y, W-135). Efficient vaccines against universal group B are in late states of development (see [12]). After vaccination, immunity develops within 7-10 days and remains effective for approximately [3][4][5] years. An understanding of the epidemiology of meningococcal epidemics has been severely curtailed by the lack of an adequate typing scheme capable of distinguishing epidemic strains of MenA.
The dynamics of MenA are complex due to a combination of various factors of societal order such as population dynamics or their susceptibility to virulent serotypes, but also climatic such as temperature and humidity, or environmental such as the presence of desert aerosols [2,5,[7][8][9].
A deep understanding of the disease dynamics would have a significant impact on the effective prevention and control strategies. Mathematical modeling and numerical simulations have the potential, and offer a promising way, to achieve this [13][14][15][16][17]. Unfortunately there are few works dealing with the design of mathematical models to study the spreading of meningococcal meningitis [18][19][20][21][22][23][24]. However, all of them are based on differential or partial differential equations ( [18][19][20][21][22][23][24]) and they exhibit some drawbacks: they do not take into account the waning of the vaccine-induced immunity and the vaccine efficacy.
For many epidemiological models, a threshold condition that indicates whether an infection introduced into a population will be eliminated or become endemic was defined [25][26][27] . The basic reproduction number R 0 is defined as the average number of secondary infections produced by an infected individual in a completely susceptible population [25]. In models with only two steady states and a transcritical bifurcation, R 0 > 1 implies that the endemic state is stable (e.g., the infection persists), and R 0 ≤ 1 implies that the uninfected state is stable (e.g., the infection will die out). The co-existence of disease-free equilibrium and endemic equilibrium points when the basic reproduction number (R 0 < 1) is typically associated with the backward or subcritical bifurcation. This phenomenon was found in many epidemiological settings (see for instance, [28][29][30][31][32][33][34][35] and references therein). The epidemiological implication of is that the classical requirement of having the associated reproduction number less than unity, while necessary is not a sufficient condition for disease control. Results showed that a threshold level of reinfection exists in all cases of the model. Beyond this threshold, the dynamics of the model are described by a backward bifurcation. However, uncertainty analysis of the parameters showed that this threshold is too high to be attained in a realistic epidemic [28].
The aim of this work is to propose a realistic mathematical model in order to identify the contribution of the vaccination to the transmission of the MenA in an endemic situation. We formulate and analyze a mathematical model for MenA that incorporates the key epidemiological and biological features of the disease such as the vaccination against MenA, the vaccine-induced immunity and the vaccine efficacy. We provide a theoretical study of the model. We provide a theoretical study of the model. The sensitivity analysis of the model is carried out to identify the most influential parameters on the model output variables. We compute the disease-free equilibrium and derive the basic reproduction number R 0 that determines the extinction and the persistence of the disease. We compute equilibria and study their stability. More precisely, we show that there exists a threshold parameter ξ such that when R 0 ≤ ξ < 1, the disease-free equilibrium is globally asymptotically stable while when ξ ≤ R 0 < 1, the model exhibits the phenomenon of backward bifurcation on a feasible region. We use this result to investigate the impact of imperfect random mass vaccination with waning immunity within a community. The sensitivity analysis of the threshold parameter ξ has been performed in order to determine which parameter makes the threshold parameter ξ as large as possible. We found that the threshold parameter increases with the rate of waning vaccine-induced immunity, the vaccination coverage rate and the natural mortality. Thus, if the vaccination coverage is sufficiently large, the disease will die out in the host population. Some discussion about the MenA persistence condition was deduced. The theoretical results are illustrated using numerical simulations.

Model description.
We consider four distinct population, according to their disease status: susceptible S (individuals which are susceptible to the disease), vaccinated individuals V (healthy individuals who have been vaccinated acquiring temporal immunity), carriers C (healthy individuals who carry MenA and are infectious) and infectious I (infected individuals who show the symptoms of the infection). Thus, the total population at time t is The model is based on the following set of biological assumptions: (i) newborns who are protected from the infection for an average period of 3 months du to maternal antibody are assumed to be susceptible to the infection; (ii) the vaccine-induced immunity acquired by vaccinated individuals is supposed to wane (iii) vaccinated individuals who vaccination is unsuccessful can catch the MenA infection and (iv) the vaccination may reduce but not completely eliminate susceptibility to infection.
Individuals are recruited into the population at constant rate Λ. The vaccinated population increases by the recruitment of individuals who are protected by the vaccination at rate πΛ where π is the proportion of recruitment who has been vaccinated, while the complementary proportion (1 − π)Λ is non-protected and enters the class of susceptible individuals. The population of vaccinated individuals is increased by the vaccination of susceptible individuals at rate ϕ. Since the vaccine does not confer a total immunity to all vaccine recipients, vaccinated individuals loss their protection and return to the susceptible class S at rate θ .
Since, the vaccination may reduce but not completely eliminate susceptibility to infection [36], we consider a factor ν as the infection rate of vaccinated members. When ν = 1, the vaccine is perfectly effective, when ν = 0, the vaccine has no effect at all. The value 1 − ν can be understood as the inefficacy level of the vaccine.
Transmission of bacteria occurs due to adequate contacts between susceptible or vaccinated individuals and carriers and infectious Thus, susceptible and vaccinated individuals are infected at rates λ S and (1 − ν)λV , respectively where λ is the force of infection given by with β the effective contact rate for MenA transmission and ε ≥ 1 a modification parameter accounting the fact that carriers (those in the C class), are more infectious than infectious (those in the I class). A proportion p of susceptible individuals who are newly infected is assumed to develop the symptoms of the infection, while the complementary part 1 − p becomes carriers.
Also, we suppose that a fraction q of vaccinated individuals who are newly infected becomes infectious and the remainder 1 − q enters the class of carriers. Du to the vaccine efficacy, one has p > q. A schematic model flowchart is depicted in Fig. 1.
The transmission model is given by the following deterministic system of nonlinear ordinary differential equations:   The parameter values used for numerical simulations are given in Table 2.

Sensitivity analysis.
We carried out the sensitivity analysis to determine the model's robustness to parameter values. That is to help us identifying the parameters that are most influential in determining disease dynamics [37]. A Latin Hypercute Sampling (LHS) scheme [38,39] samples 1000 values for each input parameter using a uniform distribution over the range of biologically realistic values, listed in Table 3 with descriptions and references given in Table 2 are computed. An output is assumed sensitive to an input if the corresponding PRCC is less than −0.50 or greater than +0.50, and the corresponding p-value is less than 5%.

PRCCs and significance
Pars

PRCCs and significance
Pars

PRCCs and significance
Pars

PRCCs and significance
Pars

Basic properties.
In this section, we study the basic results of solutions of model system (3), which are essential for the well-posedness of the system and shall be useful in the investigation of the long run dynamics. We have the following result.
Lemma 1. The non-negative orthant R 4 + is positively invariant for model system (3). Precisely, nonnegative solutions of model system (3) stay non-negative for all times.
We must prove that τ = +∞. Assume by contradiction that this is not the case, then, τ < +∞. By a continuity argument, one of the functions S(t), V (t), C(t) and I(t) will vanish at first at the first given time τ and, either remains zero for good or becomes negative at least in a small interval (τ, τ + ), in the right hand side of τ.
Without loss of generality, one can choose the first vanishing function to be, said S(t). Therefore, if τ 1 ∈ (τ, τ + ) denotes the first real number such that S(τ 1 ) = 0, we have Consequently, Relation (5) implies that there exists a positive number τ + 1 > τ 1 such that Note that S(t) is continuous, and the relations (5) and (6) , imply that τ + 1 is an extremum (more precisely, a minimum) of S(t). Moreover, since S(t) is a differentiable function on R, one hasṠ(τ + 1 ) = 0. This contradicts (4). Therefore, τ = +∞, and the proof is achieved. Similarly, one can prove that V (t) > 0, C(t) > 0 and I(t) > 0 for all t > 0. This concludes the proof.
We can now show that model system (3) is mathematically and epidemiologically well posed. In this respect, the following theorem is relevant .
is a (dissipative) dynamical system in the positively invariant region: Proof: Firstly, note that the right hand side of model system (3)  done for any initial condition in Ω so that its invariance follows immediately.
This completes the proof of the ultimate bounds, and consequently the proof the positive invariance of Ω. Now, the local existence and the boundedness of a solution ensure the global existence (in R 4 + ) of a unique solution of any initial value problem of model system (3). This achieves the proof.

DISEASE-FREE EQUILIBRIUM AND ITS STABILITY
The disease-free equilibrium (DFE) for an epidemiological model is an equilibrium such that the disease is absent in the community. Thus, is the DFE of model system (3), then C 0 = I 0 = 0. As a consequence, the disease-free equilibrium of model system (3) is It can be shown that Q 0 attracts the region: The linear stability of Q 0 is governed by the basic reproductive number R 0 which is defined as the spectral radius of the next generation matrix [26,27]. The basic reproductive number is the expected number of secondary infections generated by a single infected individual during his/her entire infected period. Using the notations in van den Driessche and Watmough [27] for model system (3), the matrices F and W for the new infection terms and the remaining transfer terms are, respectively, given by The Jacobian matrices of F and W at the disease free equilibrium Q 0 are respectively, Then, the basic reproduction ratio is defined, following [27], as the spectral radius of the next generation matrix, FW −1 : where ρ represents the spectral radius. The threshold quantity R 0 is the basic reproduction number for meningitis infection. It measures the average number of new meningitis infections generated by a single infective in a completely susceptible population.
The relevance of the basic reproduction number is due to the following result established in [27].
Lemma 2. The disease-free equilibrium Q 0 of model system (3) is locally asymptotically stable if R 0 < 1, The biological implication of Lemma 2 is that, a sufficiently small flow of infectious individuals will not generate outbreak of the disease unless R 0 > 1. For a better control on the disease, the global asymptotic stability (GAS) of the DFE is needed. Actually, enlarging the basin of attraction of Q 0 to be a part or the entire Ω is, for the model under consideration a more challenging task involving relatively new result To this end, we shall use the two outstanding theories, namely: the comparison theory for monotone dynamical systems [41] (see Smith 1995) and the theory of asymptotically autonomous systems (see [42], Castillo-Chavez and Thieme 1995) to prove the global attractivity of the disease free equilibrium.
Consider the (C, I) sub-equations in model system (3) and using the fact that S(t)/N(t) ≤ 1 and V (t)/N(t) ≤ 1 for all t ≥ 0, one has the following comparison linear system in C and I: Note that (0, 0) is the unique equilibrium of the linear system (12). System (12) can be written as where X = (C, I) T and Note that the solution X = (0, 0) is the unique equilibrium of the linear comparison system (13) which is globally asymptotically stable if and only if ρ( FW −1 ) < 1, that is Therefore, all solutions of the linear comparison system (13) converge to the trivial solution X = (0, 0) when t → ∞, with t > T . It is obvious to see that F −W as the Jacobian of model (13) is a M-matrix and a irreducible matrix [43,44]. Thus, by the comparison theorem for monotone dynamical systems [40] , we can conclude that the C and I components of model (3) also converge to zero when t → ∞, with t > T .
Now, replacing C = 0 and I = 0 in the first two equations of model system (3) gives the following (asymptotic (limiting) autonomous linear system: By Gronwall Lemma [40], one has that Taking the limit of Eq. (17), when t tend to +∞, and using the comparison theorem [40], one can conclude that S(t) converges to S 0 . Next, after plugging S 0 (see as the limiting solution of (17) ) into the second equation (V ) and using Gronwall Lemma [40], on has that Taking the limit of Eq. (18), when t tend to +∞, and using the comparison theorem [40], one can conclude that V (t) converges to V 0 .
Thus, the equilibrium (S 0 ,V 0 ) of system (16) is a globally asymptotically stable equilibrium in Ω.
Finally, by the theory of asymptotically autonomous systems (see, Castillo-Chavez and Thieme 1995 [42]), one can conclude that the DFE Q 0 is attractif in Ω. Since the locally stability has been proven and Ω is absorbant, the disease-free equilibrium Q 0 is globally asymptotically stable in R 4 + .
Thus, we have proved the following result for the global stability of the DFE Q 0 .
With this mind, one has where F 01 and F 02 are defined as in Eq. (10). Equality is only achieved when ϕ = 0, i.e., when there is no vaccination. Observe that R 0 ≤ R 0 . The constraint R 0 ≤ ξ defines implicitly a critical vaccination portion ϕ ≥ ϕ c that must achieved for eradication: Since vaccination entails costs, to choose the smallest coverage that achieves eradication could be the best option. In this way, the entire population does not need to be vaccinated in order to eradicate the disease (this is the herd immunity phenomenon). Vaccinating at the critical level ϕ c does not instantly  Table 2. It is evident that a large coverage of vaccination may dramatically decreases the number of infected individuals with and without symptoms. Also, one can observe that the time evolution of infected individuals converge slowly to the disease free equilibrium.
This means that the condition ρ > ρ c is necessary but not sufficient for the eradication of the disease within a community.
Vaccinating at the critical level ϕ v does not instantly lead to the disease eradication. Thus, to better control the infection, the condition R 0 ≤ ξ is required or the threshold parameter ξ should be as large as possible such that ζ ≈ 1. Note that the constraint R 0 ≤ ξ defines implicitly a critical vaccine efficacy ν > ν c that must be achieved for the eradication of the infection. Then, we get by Replacing ξ by its expression:   Table 2. β = 0.5 (so that ν c = 0.635) for three different values of the efficacy level: ν = 0 (so that R 0 = 2.5549), ν = 0.35 (so that R 0 = 1.3489 and ν < ν c ) and ν = 0.75 (so that R 0 = 0.5549 and ν > ν c ). All other parameter values are as in Table 2.   Table 2. Now, we are interested in which parameter that enlarge the basin of attraction of the disease-free equilibrium Q 0 . Remind that since the DFE is GAS on Ω whenever R 0 ≤ ξ , to enlarge the basin of attraction of the DFE, the threshold parameter ξ should be as large as possible. We are thus interested in which parameter that makes the threshold parameter ξ as large as possible. To do so, we perform the analysis by calculating the sensitivity indices of the threshold parameter ξ . Sensitivity analysis is used to determine the relative importance of model parameters to MenA transmission and its prevalence.
Sensitivity analysis is commonly used to determine the robustness of model predictions to parameter values, since there are usually errors in data collection and estimated values [37]. We are thus interested in parameters that significantly affect the threshold parameter ξ since they are parameters that should be taken into consideration when considering an intervention strategy. Since the threshold parameter ξ is a differentiable function of the parameters, the sensitivity index may alternatively be defined using partial derivatives. For instance, the computation of the sensitivity index of ξ with respect to π using the parameter values in Table 2 is given by This shows that ξ is an decreasing function of π and the parameter π has any influence on the threshold parameter ξ . We tabulate the indices of the remaining parameters in Table 3. From Table 8, parameters whose sensitivity indices have negative signs decrease the value of ξ as their values increase, while those with positive signs increase the value of ξ as they increase. One can observe that the parameters θ , ϕ and µ are the most influential parameters of the threshold parameter ξ . This implies that when these parameters increase, the threshold parameter ξ will also increase and this will enlarge the basin of attraction of the disease free equilibrium.

ENDEMIC EQUILIBRIUM AND ITS STABILITY
Herein, we calculate the endemic equilibrium and study its stability. To do so, we first rewrite model system (3) in the following compact form: Let Q * = (x * , y * ) be the positive endemic equilibrium of model system (24). This positive endemic equilibrium (steady state with y(t) > 0) is obtained by setting the right hand side of Eq. (24) to zero, giving: where λ * is the force of infection at the endemic equilibrium, given by With the notations of model system (24), the basic reproduction number (11) may be written as Replacing the above expression of y * into the first equation of (25) yields A simple calculation gives With this in mind, Eq. (29) becomes Solving the above Eq. (31) yields From Eq. (32), one can deduce that On the other hand, from Eq. (28), one can deduce that Combining Eqs. (26) and (34) gives Also, from Eq. (28), one can deduce that Remind that at the steady state, the size of the total population satisfies the following equation: Then, using Eq. (36), one can deduce that Equaling Eqs. (35) and (38), and using Eq. (33), it can be shown that the non-zero equilibria of model system (24) satisfy the following quadratic equation in term of λ * : with R 0 defined as in Eq. (11). Thus, positive endemic equilibria Q * are obtained by solving for λ * from the quadratic equation (39)   (ii) more than one endemic equilibrium if R 0 < 1; (iii) one or more endemic equilibria if R 0 < 1; (iv) no endemic equilibria if R 0 < 1.

The proof of case (i) of Lemma 3 is straightforward and evident. Case (ii) of Lemma 3 indicates
the possibility of a backward bifurcation in model system (3) (where the locally asymptotically stable DFE co-exists with a locally asymptotically stable endemic equilibrium when R 0 < 1, (see for instance, [28][29][30][31][32][33][34][35] and references therein) in model system (3) when R 0 < 1. Now, using the center manifold theory, we are going to show that if ξ ≤ R 0 < 1 and for a certain set of model parameters, model system (3) has exactly two endemic equilibria, with one stable and and another one unstable. To do this, we use the theorem of Castillo-Chavez and Song [35].
We have the following result.
Theorem 3. Model system (3) undergoes a backward bifurcation at R 0 = 1 if the coefficient a defined as in Eq. (46) is negative, otherwise there exists a unique endemic equilibrium Q * which is locally asymptotically stable for R 0 > 1, but close to 1.
The proof of Theorem 3 is given in the Appendix A.
The backward and forward bifurcation phenomenon is illustrated by simulating model system (3) with the parameters values of Table 2. The associated backward bifurcation diagram which occur in different signs of the coefficient a defined as in Eq. (46) is depicted in Fig. 5.   (so that ξ = 0.3929 and ξ < R 0 = 0.8579 ≤ 1). It clearly appears that when ξ ≤ R 0 < 1, the profiles can converge to either the disease-free equilibrium or an endemic equilibrium point, depending on the initial sizes of the population of the model (owing to the phenomenon of backward bifurcation). The epidemiological significance of the phenomenon of backward bifurcation is that the classical requirement of ξ ≤ R 0 < 1 is, although necessary, no longer sufficient for disease eradication. In such a scenario, disease elimination would depend on the initial sizes of the population (state variables) of the model. That is, the presence of backward bifurcation in the Men A transmission model (3) suggests that the feasibility of controlling Men A when ξ ≤ R 0 < 1 always be dependent on the initial sizes of the population.  Further, as a consequence, it is instructive to try to determine the "cause" of the backward bifurcation phenomenon in model system (3). Note that if ν = 1 which means that the vaccine is perfect, then b 2 = 0.
The role of perfect vaccine on backward bifurcation is investigated in the following. Now, let us investigate the role of imperfect vaccination. This corresponds to the case where there the vaccine is perfect within in the host population, that is , ν = 1. Then, model system (3) becomes Now, let us investigate the role of perfect vaccine. This corresponds to the case where there the vaccine is perfect within the host population, that is, ν = 1. In this case, model system (24) becomes  Table   2.
The above model has the same disease-free equilibrium Q 0 . The basic reproduction number becomes where x 0 1 = S 0 . Apart this equilibrium state, the model can also have a unique positive endemic equilibrium state. In the presence of perfect vaccine (i.e. ν = 1), the coefficients b 0 , b 1 and b 2 in Eq. (39) reduce to It is worth noting that the coefficient b 0 is positive if R 0 < 1, and negative if R 0 > 1, while b 1 is negative.
Hence, when ν = 1, no endemic equilibrium exists whenever R 0 < 1. It follows then that, owing to the absence of multiple endemic equilibria for model system (3) with ν = 1 and R 0 < 1, a backward bifurcation is unlikely to happen for model system (3). The absence of multiple endemic equilibria suggests that the disease-free equilibrium of model system (3) is globally asymptotically stable when R 0 < 1.

THE MODEL WITH STOCHASTICITY
We performed Monte-Carlos simulation on the vaccine efficacy ν to see how it affects the disease free equilibrium when R 0 < 1 for a fixed population size. Using parameter values as in Table 2, it follows that for random generation of ν i , the disease could converge to an endemic equilibrium or a disease free equilibrium (see Fig. 8). It illustrates the existence of backward bifurcation in the presence of vaccine efficacy.    Table 2. ultimately to eradicate it. The waning of recover-induced immunity was also introduced in the model.
The sensitivity analysis of the model has been investigated. We found that the model variables are most sensitive to the natural death rate, the rate of waning recover-induced immunity, the MenA induced mortality, the proportion of carriers who clear the infection and the vaccination coverage rate. The model was rigorously analyzed to gain insight into its qualitative dynamics. We computed the diseasefree equilibrium and derived the basic reproduction number R 0 that determines the outcome of MenA within the community. We prove that there exists a threshold parameter ξ such that the disease-free equilibrium is globally asymptotically stable in a feasible region whenever R 0 ≤ ξ . This result has been used to investigate the impact of imperfect random mass vaccination with waning immunity within a community. It was shown that eradication success depends on the type of vaccine as well as on the vaccination coverage. The imperfect vaccines may not completely prevent infection but could reduce the probability of being infected, thereby reducing the disease burden. Through numerical simulations, we found that the best way to control the infection is to combine a large coverage of vaccination of susceptible individuals combined with the production of a vaccine with a high level of efficacy. Also, it was mainly found that the model exhibits the phenomenon of backward bifurcation, where the stable disease-free equilibrium co-exists with a stable endemic equilibrium, when ξ ≤ R 0 ≤ 1. However, the backward bifurcation dynamics feature is caused by the vaccine efficacy had mainly occur for values that over ten times than estimations from many countries of the African meningitis belt. We also established that if R 0 > 1, the disease-free equilibrium is unstable and there exists a unique endemic equilibrium that may be locally asymptotically stable but when the basic reproduction number is close to 1. Numerical simulation have been presented to illustrate the theoretical results.

APPENDIX A: PROOF OF THEOREM 3
Herein, we give the proof of Theorem 3 on the local stability of the endemic equilibrium point Q * of model system (3). To apply this theory, the following simplification and change of variables are made first of all. Let z 1 = S, z 2 = V , z 3 = C and z 4 = I so that N = z 1 + z 2 + z 3 + z 4 . Further, by using vector notation follows: System (42) has a DFE given by Q 0 = (z 0 1 , z 0 2 , 0, 0) where z 0 . The Jacobian of model system (42), at the DFE Q 0 , is given by where F 01 and F 02 are defined as in Eq. (10). Consider, next, the case when R 0 = 1.
It follows that the Jacobian (J(Q 0 )) of system (42) at the DFE Q 0 , with β = β * , denoted by J β * has a simple zero eigenvalue (with all other eigenvalues having negative real parts). Hence, the Centre Manifold theory [38] can be used to analyze the dynamics of system (42). In particular, the theorem in Castillo and Song [39], reproduced below for convenience, will be used to show that when R 0 > 1, there exists a unique endemic equilibrium of system (42) (as shown in Lemma 1) which is locally asymptotically stable for R 0 near 1 under certain condition.
Theorem 4. (Castillo-Chavez & Song [39]): Consider the following general system of ordinary differential equations with a parameter φ : where 0 is an equilibrium point of the system (that is, f (0, φ ) ≡ 0 for all φ ) and assume (1) A = D z f (0, 0) = ∂ f i ∂ z j (0, 0) is the linearization matrix of system (43) around the equilibrium 0 with φ evaluated at 0. Zero is a simple eigenvalue of A and other eigenvalues of A have negative real parts; (2) Matrix A has a right eigenvector u and a left eigenvector v (each corresponding to the zero eigenvalue).
Let f k be the k th component of f and then, the local dynamics of the system around the equilibrium point 0 is totally determined by the signs of a and b.
In order to apply the above theorem, the following computations are necessary (it should be noted that we use β * as the bifurcation parameter, in place of φ in Theorem (Castillo-Chavez & Song [39])).
Similarly, the components of the left eigenvectors of J β * (corresponding to the zero eigenvalue), denoted by v = (v 1 , v 2 , v 3 , v 4 ) T , are given by Computation of b: For the sign of b, it can be shown that the associated non-vanishing partial derivatives of f are Substituting the respective partial derivatives into the expression where N 0 = z 0 1 + z 0 2 = Λ/µ. Then, it follows that (46) a = v 3 4