ABSTRACT

In this paper, a deterministic mathematical model for Chikungunya virus (Chikv) transmission and control is developed and analyzed to underscore the effect of vaccinating a proportion of the susceptible human, and vertical transmission in mosquito population. The disease free, and endemic equilibrium states were obtained and the conditions for the local and global stability or otherwise were given. Sensitivity analysis of the effective reproductive number, $R_c$ (the number of secondary infections resulting from the introduction of a single infected individual into a population where a proportion is fairly protected) shows that the recruitment rate of susceptible mosquito ( ${\Lambda }_M$ ) and the proportion of infectious new births from infected mosquito ( $\beta$ ) are the most sensitive parameters. Bifurcation analysis of the model using center manifold theory reveals that the model undergoes backward bifurcation (coexistence of disease free and endemic equilibrium when $R_C<1$ ). Numerical simulation of the model shows that vaccination of susceptible human population with imperfect vaccine will have a positive impact and that vertical transmission in mosquito population has a negligible effect. To the best of our knowledge, our model is the first to incorporate vaccinated human compartment and vertical transmission in (Chikv) model.

Keywords:

Chikungunya Virus, Stability, Equilibrium, Vaccination, Endemic

1. Introduction

Chikungunya is a mosquito-borne viral disease that was first observed in Tanzania in 1952 [1] . In 1964, there was epidemic of Chikungunya in Vellore, Calcutta and Maharashtra state/provinces of India [2] . Ibadan, South Western Nigeria witnessed an epidemic of Chikungunya virus in 1969 when the virus was isolated from 49 patients [3] . The disease has been identified in over 60 countries in Asia, Africa, Europe and America, and the name describes the stooping appearance of the sufferers [4] . It is an RNA virus that belongs to the alphavirus genus and the family [5] . The symptoms include abrupt onset of fever accompanied by joint pain, muscle pain headache, nausea and rash [6] . Occasionally the infection may go unrecognized or be misdiagnosed and could be acute, sub-acute and chronic.

In recent years, the virus has risen from relatively obscurity to become a global public health menace affecting millions of persons throughout the tropical and subtropical regions of the world and as such has also become a frequent cause of travel associated febrile illness [7] . The virus is transmitted through the bite of female Aedes aegypti and Aedes albopictus mosquitoes. Aedes aegypti breeds in the ubiquitous small pools of water found around human habitation [8] . Unlike Aedes aegypti which exists in tropical and subtropical area, Aedes albopictus can also thrive in temperate regions, thus potentially introducing Chikungunya to new ecological niche [9] . These species of mosquitoes are found biting throughout the daylight hours. Mother to child transmission of Chikungunya virus has been reported [10] .

Diagnosis is by confirming the presence of anti-Chikungunya antibody in the patient. At the moment, there is no vaccine or treatment for the disease. Protection is by covering of exposed skin with long pants and long sleeved shirts, insect repellents and insecticide treated mosquito nets. Since the beginning of the 19th century, mathematical model has become a veritable tool in the study of vector-borne diseases [11] [12] [13] . For (Chikv), we cite the following work, Dumont and Domerg [14] , propose a model, including human and mosquito compartments that are associated with the time course of the first epidemic of Chikungunya in Reunion Island. Using entomological results, they investigated the links between the episode of 2005 and the outbreak of 2006. Moulay, Azziz and Cadivel [15] , developed a Chikungunya transmission model for the spread of the epidemic in both humans and mosquitoes, the model involves a temporal dynamics of vector (Aedes albopictus), depending on climatic factors. In the study, they provided estimates of the transmission potential of the virus and assessed the efficacy of the measures undertaken by public health authorities to control the epidemic spread in Italy. Ruiz et al. [16] , analyzed the potential risk of Chikungunya introduction into the US, their study combines a climate-based mosquito population dynamics stochastic model with an epidemiological model to identify temporal windows that have epidemic risk.

Pongsumpun and Sangsawang [17] , model studied theoretically an age-structured model for Chikungunya involving juvenile and adult human populations, giving conditions for the disease-free and endemic states respectively. They also suggested alternative way for controlling the 8disease. Yakob and Clements [18] , analysed a simple, deterministic mathematical model for the transmission of the virus between humans and mosquitoes. They fitted the model to the large Reunion epidemic data and estimated the type reproduction number for Chikungunya, their model provided a close approximation of both the peak incidence of the outbreak and the final epidemic size.

In this work, we proposed a deterministic mathematical model for the spread, and control of Chikv. Our model attempt to bridge identified gaps in the works cited above. Specifically, our model incorporated an imperfect vaccinated human compartment and vertical transmission in the mosquito population.

2. Model Formulation

The chic model is represented by nine non-linear ordinary differential equation consisting of human-sub population and mosquito sub-population. The human sub-population is divided into; susceptible human $S_H$ , vaccinated human $V_H$ , exposed human $E_H$ , infected symptomatic human $I_1$ , infected asymptomatic human $I_2$ , recovered Human R, such that the total human population, $N_H=S_H+V_H+E_H+I_1+I_2+R$ . While the mosquito sub-population is divided into; susceptible mosquito $S_M$ , exposed mosquito $E_M$ , and infected mosquito $I_3$ , such that the total mosquito population, $N_M=S_M+E_M+I_3$ .

The parameters of the model and their values are given in Table 1, while Figure 1 is the schematic diagram of the transmission dynamics.

The susceptible human sub-population is generated at a constant rate ${\Lambda }_H$ , which includes birth and immigration. The vaccinated population is generated as members of the susceptible population receive vaccination at the rate $\nu$ , a proportion of the vaccinated with time lose their immunity at the rate $\psi$ as their vaccine wanes and move back to the susceptible population. Member of the

susceptible and vaccinated populations acquire infection at the rate $\frac{{\alpha }_1{b}_m{I}_M}{N_H}$ and $\frac{{\alpha }_1{b}_m{I}_M\left(1-\epsilon \right)}{N_H}$ respectively and move to the exposed population, where

${\alpha }_1$ is the probability of infection, $b_m$ biting rate of mosquito and $\epsilon$ (where $0<\epsilon <1$ ) is the efficacy of the imperfect vaccine. Members of the exposed population move to either symptomatic infectious population at the rate ${\sigma }_1$ or to asymptomatic infectious population at the rate $\left(1-{\sigma }_1\right)$ . The recovered population is generated as both symptomatic and asymptomatic infected populations recover with lifelong immunity at the rate $\gamma$ . All human population are decreased by natural death at the rate ${\mu }_1$ , except the two infected populations that are decreased by disease induced death at the rate $\delta$ .

The susceptible mosquito population is generated by ${\Lambda }_M$ , this population is decreased by birth from infected mosquito (vertical transmission) at the rate $\beta {\Lambda }_M$ ; and as its members take a blood meal from either symptomatic or asymptomatic infected human (horizontal transmission) at the rate ${\alpha }_2$ . The exposed mosquito population progresses to infected mosquito population at the rate ${\sigma }_2$ . It is assumed that births from infected mosquito do not pass through the exposed class. All sub-populations of mosquito die naturally at the rate ${\mu }_2$ .

2.1. The Model Equation

From the model formulation, and schematic diagram Figure 1, we hereby present the model equations.

$\frac{\text{d}{S}_H}{\text{d}t}={\Lambda }_H+\psi V_H-\frac{{\alpha }_1{b}_M{S}_H{I}_M}{N_H}-\left(\nu +{\mu }_1\right)S_H$ , (1)

$\frac{\text{d}{V}_H}{\text{d}t}=\nu S_H-\frac{{\alpha }_1{b}_M\left(1-\epsilon \right)V_H{I}_M}{N_H}-\left(\psi +{\mu }_1\right)V_H$ , (2)

$\frac{\text{d}{E}_H}{\text{d}t}=\frac{{\alpha }_1{b}_M{I}_M}{N_H}\left(S_H+\left(1-\epsilon \right)V_H\right)-\left({\sigma }_1+{\mu }_1\right)E_H$ , (3)

$\frac{\text{d}{I}_1}{\text{d}t}={\sigma }_1{E}_H-\left(\gamma +{\mu }_1+\delta \right)I_1$ , (4)

$\frac{d{I}_2}{dt}=\left(1-{\sigma }_1\right)E_H-\left(\gamma +{\mu }_1+\delta \right)I_2$ (5)

$\frac{\text{d}R}{\text{d}t}=\gamma I_1+\gamma I_2-{\mu }_{1}R$ , (6)

$\frac{\text{d}{S}_M}{\text{d}t}={\Lambda }_M-\frac{{\alpha }_2{b}_M{S}_M\left(I_1+I_2\right)}{N_H}-\beta {\Lambda }_M{I}_M-{\mu }_2{S}_M$ , (7)

$\frac{\text{d}{E}_M}{\text{d}t}=\frac{{\alpha }_2{S}_M{b}_M\left(I_1+I_2\right)}{N_H}+\beta {\Lambda }_M{I}_M-\left({\sigma }_2+{\mu }_2\right)E_M$ , (8)

$\frac{\text{d}{I}_M}{\text{d}t}={\sigma }_2{E}_M-{\mu }_2{I}_M$ . (9)

Adding (1) to (6) gives

$\frac{\text{d}{N}_H}{\text{d}t}={\Lambda }_H-\delta \left(I_1+I_2\right)-{\mu }_1{N}_H$ . (10)

Also adding (7) to (9), gives

$\frac{\text{d}{N}_M}{\text{d}t}={\Lambda }_M-{\mu }_2{N}_M$ . (11)

where

$N_H\left(t\right)=S_H\left(t\right)+V_H\left(t\right)+E_H\left(t\right)+I_1\left(t\right)+I_2\left(t\right)+R\left(t\right)$ , (12)

$N_M\left(t\right)=S_M\left(t\right)+E_M\left(t\right)+I_M\left(t\right)$ . (13)

(12) and (13) are the total human population and Aides mosquito population respectively.

2.2. Basic Properties

For the Chikungunya model (1) to (9) to be epidemiological meaningful, it is necessary to prove that all its state variables are non-negative for all time. This means that the solution of the model Equations (1) to (9) with non-negative initial data will remain non-negative for all time $t>0$ .

Lemma 1.

The closed set

$D=\left\{\begin{array}{l}\left(S_H,V_H,E_H,I_1,I_2,R,S_M,E_M,I_M\right)\in {\Re }_{+}^9:\\ S_H,V_H,E_H,I_1,I_2,R\le \frac{{\Lambda }_H}{{\mu }_1}；S_M,E_M,I_M\le \frac{{\Lambda }_M}{{\mu }_2};\end{array}\right\}$ . (14)

is positively-invariant and attracting with respect to the basic model Equations (1) to (9).

Proof

From Equations (10) and (11);

$\frac{\text{d}{N}_H}{\text{d}t}\le {\Lambda }_H-{\mu }_1{N}_H$ , $\frac{\text{d}{N}_M}{\text{d}t}\le {\Lambda }_A-{\mu }_2{N}_M$ .

It follows that $\frac{\text{d}{N}_H}{\text{d}t}<0$ and $\frac{\text{d}{N}_M}{\text{d}t}<0$ if $N_H\left(t\right)>\frac{{\Lambda }_H}{{\mu }_1}$ and $N_A\left(t\right)>\frac{{\Lambda }_M}{{\mu }_2}$ respectively. Thus a standard comparison theorem as in Lakshmikantham and Martynyuk, [25] can be used to show that

$N_H\left(t\right)\le N_H\left(0\right){\text{e}}^{-{\mu }_1\left(t\right)}+\frac{{\Lambda }_M}{{\mu }_1}\left(1-{\text{e}}^{-{\mu }_1\left(t\right)}\right)$ and $N_M\left(t\right)\le N_M\left(0\right){\text{e}}^{-{\mu }_2\left(t\right)}+\frac{{\Lambda }_M}{{\mu }_2}\left(1-{\text{e}}^{-{\mu }_2\left(t\right)}\right)$ . In particular $N_H\left(t\right)\le \frac{{\Lambda }_H}{{\mu }_1}$ and $N_M\left(t\right)\le \frac{{\Lambda }_M}{{\mu }_2}$ if $N_H\left(0\right)\le \frac{{\Lambda }_H}{{\mu }_1}$ and $N_A\left(0\right)\le \frac{{\Lambda }_M}{{\mu }_2}$ respectively. Thus D is positively-invariant. Further, if $N_H\left(0\right)>\frac{{\Lambda }_H}{{\mu }_1}$ , and $N_M\left(0\right)>\frac{{\Lambda }_M}{{\mu }_2}$ , then either the solution enters D in finite time or $N_H\left(t\right)$ approaches $\frac{{\Lambda }_H}{{\mu }_1}$ , and $N_M\left(t\right)$ approaches $\frac{{\Lambda }_M}{{\mu }_2}$ , and the infected variables $E_H,I_1,I_2,E_A,I_3$ approaches 0.

Hence D is attracting, that is all solutions in ${\Re }_{+}^9$ eventually enters D. Thus in D, the basic model Equations (1) to (9) is well posed epidemiologically and mathematically according to [26] . Hence it is sufficient to study the dynamics of the basic model Equations (1) to (9).

Lemma 2. Let the initial data $F\left(0\right)\ge 0$ ,

where

$F\left(t\right)=\left(S_H,V_H,E_H,I_1,I_2,R,S_M,E_M,I_M\right)$ .

Then the solution $F\left(t\right)$ of the Chikungunya virus model (1) to (9) are non-negative for all $t\ge 0$ . Furthermore form (10) and (11),

$\underset{t\to \infty }{\lim }\sup N_H\left(t\right)=\frac{{\Lambda }_H}{{\mu }_1+\delta }$ and $\underset{t\to \infty }{\lim }\sup N_M\left(t\right)=\frac{{\Lambda }_M}{{\mu }_2}$ .

Proof

$t_1=\sup \left\{t>0:F\left(t\right)>0\in \left[0,t\right]\right\}$ . Thus $t_1>0$ . It follows from (1) that

$\begin{array}{l}\frac{\text{d}}{\text{d}t}\left\{S_H\left(t\right)\exp \left[{\alpha }_1{b}_m{\displaystyle \underset{0}{\overset{t_1}{\int }}\frac{I_M}{N_H}\left(\xi \right)\text{d}\xi +\left(\nu +{\mu }_1\right)t}\right]\right\}\\ =\left({\Lambda }_H+\psi V_H\right)\exp \left[{\alpha }_1{b}_M{\displaystyle \underset{0}{\overset{t_1}{\int }}\frac{I_M}{N_H}\left(\xi \right)\text{d}\xi +\left(\nu +{\mu }_1\right)t}\right],\end{array}$ (15)

So that,

$\begin{array}{l}\frac{\text{d}}{\text{d}t}{S}_H\left(t_1\right)\exp \left[{\alpha }_1{b}_m{\displaystyle \underset{0}{\overset{t_1}{\int }}\frac{I_M}{N_H}\left(\xi \right)\text{d}\xi +\left(\nu +{\mu }_1\right)t_1}\right]-S_H\left(0\right)\\ ={\displaystyle \underset{0}{\overset{t_1}{\int }}\left({\Lambda }_H+\psi V_H\right)}\exp \left[{\alpha }_1{b}_m{\displaystyle \underset{0}{\overset{P}{\int }}\frac{I_M}{N_H}\left(\xi \right)\text{d}\xi +\left(\nu +{\mu }_1\right)p}\right]\text{d}p\end{array}$ (16)

Hence,

$\begin{array}{l}{S}_H\left(t_1\right)=S_H\left(0\right)\exp \left[-{\alpha }_1{b}_m{\displaystyle \underset{0}{\overset{t_1}{\int }}\frac{I_M}{N_H}\left(\xi \right)\text{d}\xi +\left(\nu +{\mu }_1\right)t_1}\right]\\ +\exp \left[-{\alpha }_1{b}_m{\displaystyle \underset{0}{\overset{t_1}{\int }}\frac{I_M}{N_H}\left(\xi \right)\text{d}\xi +\left(\nu +{\mu }_1\right)t_1}\right]\\ {\displaystyle \underset{0}{\overset{t_1}{\int }}\left({\Lambda }_H+\psi V_H\right)}\exp \left[{\alpha }_1{b}_M{\displaystyle \underset{0}{\overset{P}{\int }}\frac{I_M}{N_H}\left(\xi \right)\text{d}\xi +\left(\nu +{\mu }_1\right)p}\right]\text{d}p>0.\end{array}$ (17)

Similarly, it can be shown that $F>0$ , for all $t>0$ .

For the second part of the proof, note that,

$\begin{array}{l}0<V_H\left(t\right)\le N_H\left(t\right),0<E_H\left(t\right)\le N_H\left(t\right),0<I_1\left(t\right)\le N_H\left(t\right),\\ 0<I_2\left(t\right)\le N_H\left(t\right),0<R\left(t\right)\le N_H\left(t\right),0<S_M\left(t\right)\le N_M\left(t\right),\\ 0<E_M\left(t\right)\le N_M\left(t\right),0<I_M\left(t\right)\le N_M(t)\end{array}$

From Equations (10) and (11),

$\frac{{\Lambda }_H}{{\mu }_1+\delta }\le \underset{t\to \infty }{\lim }\inf N_H\left(t\right)\le \underset{t\to \infty }{\lim }\sup N_H\left(t\right)=\frac{{\Lambda }_H}{{\mu }_1+\delta },$ (18)

and

$\frac{{\Lambda }_M}{{\mu }_2}\le \underset{t\to \infty }{\lim }\inf N_M\left(t\right)\le \underset{t\to \infty }{\lim }\sup N_M\left(t\right)=\frac{{\Lambda }_M}{{\mu }_2}$ . (19)

as required.

3. Results

3.1. Local Stability of Disease Free Equilibrium (DFE)

The basic model (1) to (9) has a DFE, $E_0$ obtained by setting the right-hand sides of the model equations to zero, which gives:

$\begin{array}{c}{E}_0=\left(S_H^{*},V_H^{*},E_H^{*}{,}^{*}{I}_1^{*},I_2^{*},S_M^{*},E_M^{*},I_M^{*}\right)\\ =\left(\frac{{\Lambda }_H\left(\psi +{\mu }_1\right)}{\left(\psi +{\mu }_1+\nu \right){\mu }_1},\frac{\nu {\Lambda }_H}{\left(\psi +{\mu }_1+\nu \right){\mu }_1},0,0,0,0,\frac{{\Lambda }_M}{{\mu }_2},0,0\right)\end{array}$ (20)

The linear stability of $E_0$ can be established using the next generation Matrix operator method on the system (I) to (9). Using the notation in [23] , the matrices F and V, for the new infection terms and the remaining transfer terms, are, respectively, given by,

$F=\left[\begin{array}{ccccc}0& 0& 0& 0& \frac{{\alpha }_1{b}_m\left[S_H^{*}+\left(1-\epsilon \right)V_H^{*}\right]}{N_H^{*}}\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& \frac{{\alpha }_2{b}_m{S}_M^{*}}{N_H^{*}}& \frac{{\alpha }_2{b}_m{S}_M^{*}}{N_H^{*}}& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right]$ , (21)

and,

$V=\left[-\begin{array}{ccccc}{K}_3& 0& 0& 0& 0\\ -{\sigma }_1& K_4& 0& 0& 0\\ \left(1-{\sigma }_1\right)& 0& K_4& 0& 0\\ 0& 0& 0& K_5& 0\\ 0& 0& 0& -{\sigma }_2& K_6\end{array}\right]$ . (22)

where,

$\begin{array}{l}{K}_1={\nu }_1+{\mu }_1,K_2=\psi +{\mu }_1,K_3={\sigma }_1+\mu ,\\ K_4=\gamma +{\mu }_1+\delta ,K_5={\sigma }_2+{\mu }_2,K_6={\mu }_1-\beta {\Lambda }_M\end{array}$ (23)

$R_c=\frac{1}{2}\left(\frac{M_1+\sqrt{M_2+M_3}}{M_4}\right)$ , (24)

$M_1={\Lambda }_M\beta K_3{K}_4{K}_5{N}_H^{*}$ , (25)

$M_2={\Lambda }_M\beta K_3{K}_4{K}_5{N}_H^{*}$ , (26)

$M_3=4K_3{K}_4{K}_5{K}_6{\alpha }_2{\sigma }_2{b}_m{S}_M^{*}\left({\alpha }_1{b}_m{S}_H^{*}+\left(1-\epsilon \right)\right)N_H^{*}$ , (27)

$M_4=\beta K_3{K}_4{K}_5{K}_6{N}_H^{*}$ . (28)

Hence using theorem 2 of [23] the following results are established.

Theorem 1 The disease free equilibrium, $E_0$ of the model (2.1) to (2.9) is locally asymptotically stable (LAS) if $R_c<1$ , and unstable if $R_c>1$ .

3.2. Global Stability of Disease Free Equilibrium

Consider the feasible region:

$D_1=\left\{X\in D_1:S_H\le S_H^{*},V_H\le V_H^{*},R\le R^{*},S_M\le S_M^{*}\right\},$ (29)

$X=\left\{S_H,V_H,E_H,I_1,I_2,R,S_M,E_M,I_M\right\}.$ (30)

Lemma 3. The region $D_1$ is positively invariant for the Chikungunya model Equations (1) to (9).

Proof

From Equations (1) to (9) and (20),

we have that, the only non-zero compartments at disease free equilibrium are;

$\begin{array}{l}\frac{d{S}_H}{dt}={\Lambda }_H+\psi V_H-\frac{{\alpha }_1{b}_M{S}_H{I}_M}{N_H}-\left(\nu +{\mu }_1\right)S,\\ \frac{d{V}_H}{dt}=\nu S_H-\frac{{\alpha }_1{b}_M\left(1-\epsilon \right)V_H{I}_M}{N_H}-\left(\psi +{\mu }_1\right)V_H,\\ \frac{d{S}_M}{dt}={\Lambda }_M-\frac{{\alpha }_2{b}_M{S}_M\left(I_1+I_2\right)}{N_H}-\beta {\Lambda }_M{I}_M-{\mu }_2{S}_M\end{array}$ (31)

Such that,

$\begin{array}{c}\frac{\text{d}{S}_H}{\text{d}t}={\Lambda }_H+\psi V_H-\frac{{\alpha }_1{b}_m{S}_H{I}_M}{N_H}-\left(\nu +{\mu }_1\right)S_H,\\ \le {\Lambda }_H+\psi V_H-\left(\nu +{\mu }_1\right)S_H\\ \le \left(\nu +{\mu }_1\right)\left[\frac{{\Lambda }_H\left(\psi +{\mu }_1\right)}{\left(\psi +{\mu }_1+\nu \right){\mu }_1}+\psi \frac{\nu {\Lambda }_H}{\left(\psi +{\mu }_1+\nu \right){\mu }_1}-S_H\right]\\ =\left(\nu +{\mu }_1\right)\left(S_H^{*}+\psi V_H^{*}-S_H\right)\end{array}$ , (32)

Hence,

$S_H\left(t\right)\le S_H^{*}+\psi V_H^{*}-\left[S_H^{*}-\psi V_H^{*}-S_H\left(0\right)\right]{\text{e}}^{-\left(\nu +{\mu }_1\right)t}.$ (33)

Thus if $N_H^{*}=\frac{{\Lambda }_H}{{\mu }_1}$ and $S_H\left(0\right)\le S_H^{*}+\psi V_H^{*}$ for all $t\ge 0$ , then $S_H\left(t\right)\le S_H^{*}+\psi V_H^{*}$ for all $t\ge 0$ .

Similarly, it follows from Equation (7) of our model and (20) where $S_M^{*}=\frac{{\Lambda }_M}{{\mu }_2}$ .

We have that,

$\begin{array}{c}\frac{\text{d}{S}_M}{\text{d}t}={\Lambda }_M-\frac{{\alpha }_2{b}_M{S}_M\left(I_1+I_2\right)}{N_H}-\beta {\Lambda }_M{I}_M-{\mu }_2{S}_M\\ \le {\Lambda }_M-{\mu }_2{S}_M\le {\mu }_2\left[\frac{{\Lambda }_M}{{\mu }_2}-S_M\right]={\mu }_2\left(S_M^{*}-S_M\right)\end{array}$ . (34)

Hence,

$S_M\left(t\right)\le S_M^{*}-\left[S_M^{*}-S_M\left(0\right)\right]{\text{e}}^{-{\mu }_{2}t}$ . (35)

Thus if $N_M^{*}=\frac{{\Lambda }_M}{{\mu }_2}$ and $S_M\left(0\right)\le S_M^{*}$ for all $t\ge 0$ , then $S_M\left(t\right)\le S_M^{*}$ for all $t\ge 0$ .

In summary, we have shown that $D_1$ is positively invariant and attracting with respect to the solutions of our model Equations (1) to (9).

Theorem 2

The DFE of the basic model (1) to (9) is Global Asymptotical Stability (GAS) in $D_1$ , whenever $R_C\le 1$ .

Proof

To prove the GAS of the DFE we adopt the approach in [27] .

Let $X=\left(S_H,V_H,R,S_M\right)$ and $Z=\left(E_H,I_1,I_2,E_M,I_M\right)$ and group our model Equations (1) to (8) into:

$\begin{array}{l}\frac{\text{d}X}{\text{d}t}=F\left(X,0\right),\\ \frac{\text{d}Z}{\text{d}t}=G\left(X,Z\right).\end{array}$ (36)

where $F\left(X,0\right)$ is the right hand side of ${\stackrel{\dot{}}{S}}_H,{\stackrel{\dot{}}{V}}_H,\stackrel{\dot{}}{R},{\stackrel{\dot{}}{S}}_M$ with $E_H=I_1=I_2=E_M=I_M=0$ and $G\left(X,Z\right)$ , the right hand side of ${\stackrel{\dot{}}{E}}_H,{\stackrel{\dot{}}{I}}_1,{\stackrel{\dot{}}{I}}_2,{\stackrel{\dot{}}{E}}_M,{\stackrel{\dot{}}{I}}_M$ . Next we consider the reduced system:

$\frac{\text{d}X}{\text{d}t}=F\left(X,0\right)$ given as,

$\begin{array}{l}\frac{\text{d}{S}_H}{\text{d}t}={\Lambda }_H-{\mu }_1{S}_H,\\ \frac{\text{d}{V}_H}{\text{d}t}=\nu S_H-\left(\psi +{\mu }_1\right)V_H,\\ \frac{\text{d}R}{\text{d}t}=-{\mu }_{1}R,\\ \frac{\text{d}{S}_M}{\text{d}t}={\Lambda }_M-{\mu }_2{S}_M.\end{array}$ (37)

Let $X^{*}=\left(S_H^{*},V_H^{*},R^{*},S_M^{*}\right)=\left(\frac{{\Lambda }_H\left(\psi +{\mu }_1\right)}{\left(\psi +{\mu }_1+\nu \right){\mu }_1},\frac{\nu {\Lambda }_H}{\left(\psi +{\mu }_1+\nu \right){\mu }_1},0,\frac{{\Lambda }_M}{{\mu }_2}\right)$ (38)

be an equilibrium of (37) we show that $X^{*}$ is a global stable equilibrium in $D_1$ .

To do this, we solve the Equations (37), which gives

$\begin{array}{l}{S}_H\left(t\right)=\left(\frac{{\Lambda }_H\left(\psi +{\mu }_1\right)}{\left(\psi +{\mu }_1+\nu \right){\mu }_1}+\psi V_H^{*}\right)-\left(\frac{{\Lambda }_H\left(\psi +{\mu }_1\right)}{\left(\psi +{\mu }_1+\nu \right){\mu }_1}+\psi V_H^{*}\right){\text{e}}^{-\left(\left(\psi +{\mu }_1+\nu \right){\mu }_1\right)t}\\ +S_H\left(0\right){\text{e}}^{\left(\left(\psi +{\mu }_1+\nu \right){\mu }_1\right)t},\\ S_H\left(t\right)\to \frac{{\Lambda }_H\left(\psi +{\mu }_1\right)}{\left(\psi +{\mu }_1+\nu \right){\mu }_1}+\psi V_H^{*},\end{array}$ (39)

as $t\to \infty$ .

$\begin{array}{l}{V}_H\left(t\right)=\frac{\nu {\Lambda }_H}{\left(\psi +{\mu }_1+\nu \right){\mu }_1}-\frac{\nu {\Lambda }_H}{\left(\psi +{\mu }_1+\nu \right){\mu }_1}{\text{e}}^{-\left(\left(\psi +{\mu }_1+\nu \right){\mu }_1\right)t}+V_H\left(0\right){\text{e}}^{\left(\left(\psi +{\mu }_1+\nu \right){\mu }_1\right)t},\\ V_H\left(t\right)\to \frac{\nu {\Lambda }_H}{\left(\psi +{\mu }_1+\nu \right){\mu }_1},\end{array}$ (40)

as $t\to \infty$ .

$R\left(t\right)=R\left(0\right){\text{e}}^{-{\mu }_{1}t},R\left(t\right)\to 0$ , (41)

as $t\to \infty$ .

$S_M\left(t\right)=\frac{{\Lambda }_M}{{\mu }_2}-\frac{{\Lambda }_M}{{\mu }_2}{\text{e}}^{-{\mu }_{2}t}+S_M\left(0\right){\text{e}}^{-{\mu }_{2}t},S_M\left(t\right)\to \frac{{\Lambda }_M}{{\mu }_2}$ , (42)

as $t\to \infty$ .

This asymptotic dynamics is independent of initial conditions in D. Hence the solution of xxx converges globally in $D_1$ .

Next we are required to show that $G\left(X,Z\right)$ satisfies the following two conditions in [19] pp246 namely;

$\begin{array}{l}G\left(X,0\right)=0,\\ G\left(X,Z\right)=D_Z\hat{G}\left(X^{*},0\right)Z-\hat{G}\left(X,Z\right),\hat{G}\left(X,Z\right)\ge 0,\end{array}$ (43)

where,

$\left(X^{*},0\right)=\left(\frac{{\Lambda }_H\left(\psi +{\mu }_1\right)}{\left(\psi +{\mu }_1+\nu \right){\mu }_1},\frac{\nu {\Lambda }_H}{\left(\psi +{\mu }_1+\nu \right){\mu }_1},0,\frac{{\Lambda }_M}{{\mu }_2}\right)$ . (44)

and $D_{Z}G\left(X^{*},0\right)$ is the Jacobian of $G\left(X,Z\right)$ taken with respect to $\left(E_H,I_1,I_2,E_M,I_M\right)$ and evaluated at $\left(X^{*},0\right)$ , which is an M-Matrix (the off diagonal elements are non-negative).

Thus,

$D_{Z}G\left(X^{*},0\right)=\left(\begin{array}{ccccc}-k_3& 0& 0& 0& Q_1\\ {\sigma }_1& -k_4& 0& 0& 0\\ 1-{\sigma }_1& 0& -k_4& 0& 0\\ 0& \frac{{\alpha }_2{b}_m{S}_M^{*}}{N_H^{*}}& \frac{{\alpha }_2{b}_m{S}_M^{*}}{N_H^{*}}& -k_5& 0\\ 0& 0& 0& {\sigma }_2& -k_6\end{array}\right)$ , (45)

$\hat{G}\left(X,Z\right)=\left(\begin{array}{ccccc}0& 0& 0& 0& Q_2{I}_M\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& {\alpha }_2{b}_M\frac{N_H^{*}}{S_M^{*}}{Q}_3& {\alpha }_2{b}_M\frac{N_H^{*}}{S_M^{*}}{Q}_3& 0& 0\\ 0& 0& 0& 0& \frac{1}{\beta {\Lambda }_M}\end{array}\right)$ , (46)

where,

$\begin{array}{l}{Q}_1=\frac{{\alpha }_1{b}_m{S}_H^{*}+\left(1-\epsilon \right)+V_H^{*}}{N_H^{*}},\\ Q_2=\frac{N_H^{*}}{S_H^{*}+\left(1-\epsilon \right)+V_H^{*}}\left(1-\frac{N_H^{*}}{S_H^{*}+\left(1-\epsilon \right)+V_H^{*}}\frac{S_H+\left(1-\epsilon \right)V_H}{N_H}\right),\\ Q_3=\left(1-\frac{N_H^{*}}{S_M^{*}}\frac{S_M}{N_H}{I}_1\right),\\ Q_4=\left(1-\frac{N_H^{*}}{S_M^{*}}\frac{S_M}{N_H}{I}_2\right).\end{array}$ (47)

Further $S_H\le S_H^{*}$ , $V_H\le V_H^{*}$ and $S_M\le S_M^{*}$ in $D_1$ . Thus, it follows that $\left(1-\frac{S_H}{S_H^{*}}\right)>0$ , $\left(1-\frac{V_H}{V_H^{*}}\right)>0$ and $\left(1-\frac{S_M}{S_M^{*}}\right)>0$ . Hence $\hat{G}\left(X,Z\right)\ge 0$ .

Therefore, by the theorem 2 in [28] , the disease-free equilibrium is globally asymptotically stable since in the absence of disease induced mortality the human population is constant.

3.3. Sensitivity Analysis

Here we present the sensitivity index of the parameters of the effective reproductive number $\left(R_C\right)$ . Sensitivity tells us how important each parameter is to disease transmission. Such information, is crucial not only to experimental design, but also to data assimilation and reduction of complex nonlinear model [29] . Sensitivity Analysis is commonly used to determine the robustness of model prediction to parameter values, since there are usually errors in data collection and presumed parameter values. It is used to determine parameters that have high impact on the $\left(R_C\right)$ and should be targeted by intervention strategies. Sensitivity indexes allows us to measure the relative changes in a variable when a parameter changes. The normalized forward sensitivity index of a variable with respect to a parameter is the ratio of relative changes in the parameter when the variable is a differentiable function of the parameter. The sensitivity index may be alternatively defined using partial derivatives. The sensitivity index of our model is given in Table 2.

From Table 2, the most sensitive parameter of $R_C$ is the recruitment rate of susceptible mosquito ( ${\Lambda }_M$ ) followed by the proportion of infectious new birth from infected mosquito ( $\beta$ ) while the natural birth rate of mosquito ( ${\mu }_2$ ) and the rate at which exposed mosquito become infectious ( ${\sigma }_2$ ) are equally sensitive to the $R_C$ according to the model. This means that any policy or practice capable of reducing these parameters will go a long way in reducing the menace of Chikungunya and at the long run, result to eradication.

Endemic Equilibrium

Let $E_1=\left(S_H^{**},V_H^{**},E_H^{**},R^{**},I_1^{**},I_2^{**},R^{**},S_M^{**},E_M^{**},I_M^{**}\right)$ , (48)

represents any arbitrary endemic equilibrium of the model (1) to (9). Further, let

${\lambda }_H^{**}=\frac{{\alpha }_1{b}_m{I}_M^{**}}{N_H^{**}},{\lambda }_M^{**}=\frac{{\alpha }_2{b}_M\left(I_1^{**}+I_2^{**}\right)}{N_H^{**}}.$ (49)

be the forces of infection of humans and vectors at steady state, respectively. Solving (1) to (9) in terms of ${\lambda }_H^{**}$ and ${\lambda }_M^{**}$ , we have;

$\begin{array}{l}{S}_H^{**}=\frac{{\Lambda }_H\left({\lambda }_H^{**}+k_2\right)}{\left({\lambda }_H^{**}+k_2\right)\left({\lambda }_H^{**}+k_1\right)+\psi \nu },V_H^{**}=\frac{\nu {\Lambda }_H}{\left({\lambda }_H^{**}+k_2\right)\left({\lambda }_H^{**}+k_1\right)+\psi \nu },\\ E_H^{**}=\frac{{\Lambda }_H{\lambda }_H^{**}\left({\lambda }_H^{**}+k_2\right)+\left(1-\epsilon \right)}{\left(\left({\lambda }_H^{**}+k_2\right)\left({\lambda }_H^{**}+k_1\right)-\psi \nu \right)k_3},I_1^{**}=\frac{{\sigma }_1{\Lambda }_H{\lambda }_H^{**}\left(\left({\lambda }_H^{**}+k_2\right)+\left(1-\epsilon \right)\right)\nu }{\left(\left({\lambda }_H^{**}+k_2\right)\left({\lambda }_H^{**}+k_1\right)-\psi \nu \right)k_3{k}_4},\\ I_2^{**}=\frac{\left(1-{\sigma }_1\right){\Lambda }_H{\lambda }_H^{**}\left(\left({\lambda }_H^{**}+k_2\right)+\left(1-\epsilon \right)\right)\nu }{\left(\left({\lambda }_H^{**}+k_2\right)\left({\lambda }_H^{**}+k_1\right)-\psi \nu \right)k_3{k}_4},R^{**}=\frac{\gamma {\Lambda }_H{\lambda }_H^{**}\left(\left({\lambda }_H^{**}+k_2\right)+\left(1-\epsilon \right)\right)\nu }{\left(\left({\lambda }_H^{**}+k_2\right)\left({\lambda }_H^{**}+k_1\right)-\psi \nu \right)k_3{k}_4{\mu }_1},\\ S_M^{**}=\frac{k_5{k}_6{\Lambda }_M}{{\lambda }_M^{**}\left(k_5{k}_6+\beta {\Lambda }_M{\sigma }_2\right)+{\mu }_2{k}_5{k}_6},E_M^{**}=\frac{k_6{\Lambda }_M{\lambda }_M^{**}}{{\lambda }_M^{**}\left(k_5{k}_6+\beta {\Lambda }_M{\sigma }_2\right)+{\mu }_2{k}_5{k}_6},\\ I_M^{**}=\frac{{\sigma }_2{\Lambda }_M{\lambda }_M^{**}}{{\lambda }_M^{**}\left(k_5{k}_6+\beta {\Lambda }_M{\sigma }_2\right)+{\mu }_2{k}_5{k}_6}.\end{array}$ (50)

Substituting (20) into (19) we have;

$\begin{array}{l}{\lambda }_M^{**}=\frac{{\alpha }_2{b}_m{\Lambda }_M{\lambda }_H^{**}\left(\left({\lambda }_H^{**}+k_2\right)+\left(1-\epsilon \right)\right)\nu }{\left(\left({\lambda }_H^{**}+k_2\right)\left({\lambda }_H^{**}+k_1\right)-\psi \nu \right)k_3{k}_4},\\ {\lambda }_H^{**}=A{\left({\lambda }_H^{**}\right)}^4+B{\left({\lambda }_H^{**}\right)}^3+C{\left({\lambda }_H^{**}\right)}^2+D\left({\lambda }_H^{**}\right)-E.\end{array}$ (51)

where,

$A=\left({\alpha }_2{b}_m{\Lambda }_H{k}_5{k}_6+\beta {\Lambda }_M{\sigma }_2\right)\left({\mu }_1+k_4{\mu }_1+\gamma \right)$ , (52)

$B=\left(T_2{k}_5{k}_6+\beta {\Lambda }_M{\sigma }_2\right)\left(k_3{k}_4\left({\mu }_1{k}_2+\psi \right)\right)-\left({\alpha }_1{\alpha }_2{\left(b_m\right)}^2{\Lambda }_H{\sigma }_2\right)k_3{k}_4{\mu }_1$ , (53)

$\begin{array}{c}C=[{\alpha }_2{b}_m{\Lambda }_H{\sigma }_2\left(k_2+\left(1-\epsilon \right)\nu \right)k_5{k}_6\\ +\beta {\Lambda }_M{\sigma }_2\left(k_3{k}_4{\mu }_1{k}_2+k_1{\mu }_1+k_4{\mu }_1+\gamma \left(1-\epsilon \right)\nu \right)]\\ -{\alpha }_1{\alpha }_2{\left(b_m\right)}^2{\Lambda }_H{\sigma }_2{k}_3{k}_4{\mu }_1\left(k_1+k_2\right)\\ +{\alpha }_1{\alpha }_2{\left(b_m\right)}^2{\Lambda }_H{\sigma }_2\left(k_2+\left(1-\epsilon \right)\nu \right)k_3{k}_4{\mu }_1,\end{array}$ (54)

$\begin{array}{l}D=\left({\alpha }_2{b}_m{\Lambda }_H{\sigma }_2\left(k_2+\left(1-\epsilon \right)\nu \right)\left({\alpha }_2{b}_m{\Lambda }_H{k}_5{k}_6+\beta {\Lambda }_M{\sigma }_2\right)\left(k_3{k}_4\left({\mu }_1{k}_2+\nu \right)\right)\right)\\ -\left({\alpha }_1{\alpha }_2{\left(b_m\right)}^2{\Lambda }_H{\sigma }_2\left(k_2+\left(1-\epsilon \right)\nu \right)\left(k_1+k_2\right)k_3{k}_4{\mu }_1\right)\\ +\left(k_1{k}_2{k}_3{k}_4-\psi \nu k_3{k}_4{\mu }_1\right)\left({\alpha }_1{\alpha }_2{\left(b_m\right)}^2{\Lambda }_H{\sigma }_2\right),\end{array}$ (55)

$E={\alpha }_1{\alpha }_2{\left(b_m\right)}^2{\Lambda }_H{\sigma }_2\left(k_2+\left(1-\epsilon \right)\nu \right)k_1{k}_2{k}_3{k}_4-\psi \nu k_3{k}_4{\mu }_1$ (56)

Theorem 3.6. The Chikungunya basic model (1) to (9) undergoes backward bifurcation whenever the coefficient a in equation is positive.

Proof. To prove this theorem, we use the Centre Manifold theory as in Castillo-Chavez and songs [30] [31] see the theorem in Appendix A.

Let $S_H=x_1$ , $V_H=x_2$ , $E_H=x_3$ , $I_1=x_4$ , $I_2=x_5$ , $R=x_6$ , $S_M=x_7$ , $E_M=x_8$ , $I_M=x_9$ so that $N_H=x_1+x_2+x_3+x_4+x_5+x_6$ and $N_M=x_7+x_8+x_9$ . Further by using vector notation

$X={\left(x_1+x_2+x_3+x_4+x_5+x_6+x_7+x_8+x_9\right)}^{\text{T}}$ Equations (1) to (9) can be written

as $\frac{\text{d}X}{\text{d}t}={\left(f_1+f_2+f_3+f_4+f_5+f_6+f_7+f_8+f_9\right)}^{\text{T}}$ as follow:

$\begin{array}{l}\frac{\text{d}{x}_1}{\text{d}t}={\Lambda }_H+\psi V_H-\frac{{\alpha }_1{b}_M{x}_1{x}_9}{x_1+x_2+x_3+x_4+x_5+x_6}-k{x}_1,\\ \frac{\text{d}{x}_2}{\text{d}t}=\nu x_2-\frac{{\alpha }_1{b}_M\left(1-\epsilon \right)x_2{x}_9}{x_1+x_2+x_3+x_4+x_5+x_6}-k_2{x}_2,\\ \frac{\text{d}{x}_3}{\text{d}t}=\frac{{\alpha }_1{b}_M{x}_9}{x_1+x_2+x_3+x_4+x_5+x_6}\left(x_1+\left(1-\epsilon \right)x_2\right)-k_3{x}_3,\\ \frac{\text{d}{x}_4}{\text{d}t}={\sigma }_1{x}_3-k_4{x}_4,\\ \frac{\text{d}{x}_5}{\text{d}t}=\left(1-{\sigma }_1\right)x_3-\left(\gamma +{\mu }_1+\delta \right)x_5,\\ \frac{\text{d}{x}_6}{\text{d}t}=\gamma x_4+\gamma x_5-{\mu }_1{x}_6,\\ \frac{\text{d}{x}_7}{\text{d}t}={\Lambda }_M-\frac{{\alpha }_2{b}_M{x}_7\left(x_4+x_5\right)}{x_1+x_2+x_3+x_4+x_5+x_6}-\beta {\Lambda }_M{x}_9-{\mu }_2{x}_7,\\ \frac{\text{d}{x}_8}{\text{d}t}=\frac{{\alpha }_2{S}_M{b}_M\left(x_4+x_5\right)}{x_1+x_2+x_3+x_4+x_5+x_6}+\beta {\Lambda }_M{x}_9-k_5{x}_8,\\ \frac{\text{d}{x}_9}{\text{d}t}={\sigma }_2{x}_8-k_6{x}_9.\end{array}\}$ (57)

Because it’s not always convenient to use $R_C=1$ as bifurcation parameter, we choose $P=P^{*}$ where $P^{*}={\alpha }_2{b}_m$ as the bifurcation parameter such that,

$P^{*}=\frac{1}{2}\left(\frac{M_4}{M_1+\sqrt{M_2+M_5}}\right)$ , (58)