Modelling the Effects of Vertical Transmission in Mosquito and the Use of Imperfect Vaccine on Chikungunya Virus Transmission Dynamics

Show more

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}$ .

Table 1. Parameters of the model Equations (1) to (9).

Figure 1. Schematic diagram of Chikungunya virus transmission dynamics, Equations (1) to (9).

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}}\uff1b{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}{\mathrm{lim}}\mathrm{sup}{N}_{H}\left(t\right)=\frac{{\Lambda}_{H}}{{\mu}_{1}+\delta}$ and $\underset{t\to \infty}{\mathrm{lim}}\mathrm{sup}{N}_{M}\left(t\right)=\frac{{\Lambda}_{M}}{{\mu}_{2}}$ .

Proof

${t}_{1}=\mathrm{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)\mathrm{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)\mathrm{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)\mathrm{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)}\mathrm{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)\mathrm{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]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\mathrm{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]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\displaystyle \underset{0}{\overset{{t}_{1}}{\int}}\left({\Lambda}_{H}+\psi {V}_{H}\right)}\mathrm{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}{\mathrm{lim}}\mathrm{inf}{N}_{H}\left(t\right)\le \underset{t\to \infty}{\mathrm{lim}}\mathrm{sup}{N}_{H}\left(t\right)=\frac{{\Lambda}_{H}}{{\mu}_{1}+\delta},$ (18)

and

$\frac{{\Lambda}_{M}}{{\mu}_{2}}\le \underset{t\to \infty}{\mathrm{lim}}\mathrm{inf}{N}_{M}\left(t\right)\le \underset{t\to \infty}{\mathrm{lim}}\mathrm{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}=4{K}_{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}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.05em}}+{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}\stackrel{^}{G}\left({X}^{*},0\right)Z-\stackrel{^}{G}\left(X,Z\right),\stackrel{^}{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)

$\stackrel{^}{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 $\stackrel{^}{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.

Table 2. Sensitivity analysis index for the effective basic reproductive number.

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}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\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)]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-{\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)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+{\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)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}-\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)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\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)

where

${M}_{5}=4{k}_{3}{k}_{4}{k}_{5}{k}_{6}{\alpha}_{2}{x}_{1}^{*}\left({\alpha}_{7}{b}_{m}{x}_{1}^{*}+\left(1-\epsilon \right){x}_{2}^{*}\right)$ . (59)

The Jacobian of (57) evaluated at ${E}_{0}$ with ${\alpha}_{2}{b}_{m}={P}^{*}$ , denoted by ${J}^{*}$ is given

${J}^{*}=\left[\begin{array}{ccccccccc}-{k}_{1}& 0& 0& 0& 0& 0& 0& 0& -{Q}_{5}\\ \nu & -{k}_{2}& 0& 0& 0& 0& 0& 0& -{Q}_{6}\\ 0& 0& -{k}_{3}& 0& 0& 0& 0& 0& 0\\ 0& 0& {\sigma}_{1}& -{k}_{4}& 0& 0& 0& 0& 0\\ 0& 0& 1-{\sigma}_{1}& 0& -{k}_{4}& 0& 0& 0& 0\\ 0& 0& 0& \gamma & \gamma & -{\mu}_{1}& 0& 0& 0\\ 0& 0& 0& {Q}_{7}& {Q}_{7}& 0& {\mu}_{2}& 0& -\beta {\Lambda}_{M}\\ 0& 0& 0& {Q}_{7}& {Q}_{7}& 0& 0& -{k}_{5}& 0\\ 0& 0& 0& 0& 0& 0& 0& {\sigma}_{2}& -{k}_{6}\end{array}\right]$ (60)

where,

$\begin{array}{l}{Q}_{5}=\frac{{\alpha}_{1}{b}_{m}{x}_{1}^{*}}{{x}_{1}+{x}_{2}+{x}_{3}+{x}_{4}+{x}_{5}+{x}_{6}},\\ {Q}_{6}=\frac{{\alpha}_{1}{b}_{m}\left(1-\epsilon \right){x}_{2}^{*}}{{x}_{1}+{x}_{2}+{x}_{3}+{x}_{4}+{x}_{5}+{x}_{6}},\\ {Q}_{7}=\frac{{\alpha}_{2}{b}_{m}{x}_{7}^{*}}{{x}_{1}+{x}_{2}+{x}_{3}+{x}_{4}+{x}_{5}+{x}_{6}}.\end{array}$ (61)

It follows that (60) has a right eigenvector denoted by $v={v}_{1},{v}_{2},{v}_{3},{v}_{4},{v}_{5},{v}_{6},{v}_{7},{v}_{8},{v}_{91}$ , where

$\begin{array}{l}{v}_{1}=\frac{-{Q}_{5}\left({k}_{1}{k}_{2}+\nu \psi \right)+\psi \left({Q}_{6}{k}_{1}+\nu {Q}_{5}\right){v}_{9}}{\left({k}_{1}{k}_{2}-\nu \psi \right){k}_{1}},\\ {v}_{2}=\frac{-\left({Q}_{6}{k}_{1}+\nu {Q}_{5}\right){v}_{9}}{{k}_{1}{k}_{2}-\nu \psi},\\ {v}_{3}=\frac{{Q}_{1}{v}_{9}}{{k}_{3}},\\ {v}_{4}=\frac{{\sigma}_{1}{Q}_{1}{v}_{9}}{{k}_{3}{k}_{4}},\\ {v}_{5}=\frac{\left(1-{\sigma}_{1}\right){Q}_{1}{v}_{9}}{{k}_{3}{k}_{4}},\end{array}$

$\begin{array}{l}{v}_{6}=\frac{\gamma {Q}_{1}{k}_{4}\left({Q}_{1}+\left(1-{\sigma}_{1}\right)\right){v}_{9}}{{k}_{3}{k}_{4}^{2}},\\ {v}_{7}=\frac{\left(\beta {\Lambda}_{M}{k}_{4}^{2}{k}_{3}+{Q}_{7}{Q}_{1}{k}_{1}\right){v}_{9}}{{k}_{3}{k}_{4}^{2}{\mu}_{2}},\\ {v}_{8}=\frac{{Q}_{7}{Q}_{1}{k}_{4}\left({\sigma}_{1}+\left(1-{\sigma}_{1}\right)\right){v}_{9}}{{k}_{3}{k}_{4}^{2}{k}_{5}},\\ {v}_{9}={v}_{9}.\end{array}$ (62)

And a left eigenvector given by $w={w}_{1},{w}_{2},{w}_{3},{w}_{4},{w}_{5},{w}_{6},{w}_{7},{w}_{8},{w}_{91}$ , where

$\begin{array}{l}{w}_{1}=\frac{\nu {w}_{2}}{{k}_{1}},\\ {w}_{2}={w}_{2},\\ {w}_{3}=\frac{{Q}_{7}{w}_{7}+\gamma {w}_{6}}{{k}_{3}{k}_{4}},\\ {w}_{4}={w}_{5}=\frac{\gamma {w}_{6}-{Q}_{7}{w}_{7}}{{k}_{3}{k}_{4}},\end{array}$

$\begin{array}{l}{w}_{6}={w}_{6},\\ {w}_{7}={w}_{7},\\ {w}_{8}={w}_{9}=0.\end{array}$ (63)

Computation of a

$\frac{{\partial}^{2}{f}_{1}}{\partial {x}_{1}\partial {x}_{9}}=\frac{{\alpha}_{1}{b}_{m}{x}_{1}^{*}}{{\left({x}_{1}^{*}+{x}_{2}^{*}\right)}^{2}}-\frac{{\alpha}_{1}{b}_{m}}{{x}_{1}^{*}+{x}_{2}^{*}},$

$\begin{array}{l}\frac{{\partial}^{2}{f}_{1}}{\partial {x}_{2}\partial {x}_{9}}=\frac{{\partial}^{2}{f}_{1}}{\partial {x}_{3}\partial {x}_{9}}=\frac{{\partial}^{2}{f}_{1}}{\partial {x}_{4}\partial {x}_{9}}=\frac{{\partial}^{2}{f}_{1}}{\partial {x}_{6}\partial {x}_{9}}=\frac{{\alpha}_{1}{b}_{m}{x}_{1}^{*}}{{\left({x}_{1}^{*}+{x}_{2}^{*}\right)}^{2}},\\ \frac{{\partial}^{2}{f}_{1}}{\partial {x}_{1}\partial {x}_{2}}=\frac{{k}_{1}-\psi}{{x}_{1}^{*}+{x}_{2}^{*}}+\frac{2\left(\psi {x}_{2}^{*}-{k}_{3}{x}_{1}^{*}\right)}{{\left({x}_{1}^{*}+{x}_{2}^{*}\right)}^{3}},\\ \frac{{\partial}^{2}{f}_{7}}{\partial {x}_{4}\partial {P}^{*}}=\frac{{\partial}^{2}{f}_{7}}{\partial {x}_{5}\partial {P}^{*}}=\frac{{\partial}^{2}{f}_{7}}{\partial {x}_{7}\partial {P}^{*}}=\frac{{x}_{7}}{{x}_{1}^{*}+{x}_{2}^{*}},\\ \frac{{\partial}^{2}{f}_{7}}{\partial {x}_{4}\partial {x}_{7}}=\frac{{\partial}^{2}{f}_{7}}{\partial {x}_{5}\partial {x}_{7}}=\frac{{\mu}_{2}}{{\left({x}_{1}^{*}+{x}_{2}^{*}\right)}^{2}}-\frac{{P}^{*}}{{x}_{1}^{*}+{x}_{2}^{*}}.\end{array}$ (64)

$\begin{array}{c}a={\displaystyle \underset{k,i,j=1}{\overset{n}{\sum}}{v}_{k}{w}_{i}{w}_{j}\frac{{\partial}^{2}{f}_{k}}{\partial {x}_{i}\partial {x}_{j}}\left(0,0\right)}\\ ={v}_{1}{w}_{9}\frac{{\alpha}_{1}{b}_{m}{x}_{1}^{*}}{{\left({x}_{1}^{*}+{x}_{2}^{*}\right)}^{2}}\left({w}_{1}+{w}_{2}+{w}_{3}+{w}_{4}+{w}_{5}+{w}_{6}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-{v}_{1}{w}_{2}\left({w}_{9}\frac{{\alpha}_{1}{b}_{m}}{{x}_{1}^{*}+{x}_{2}^{*}}+{w}_{2}\left(\frac{{k}_{3}+\psi}{{x}_{1}^{*}+{x}_{2}^{*}}+\frac{2\psi {x}_{2}^{*}+{k}_{3}{x}_{1}^{*}}{{\left({x}_{1}^{*}+{x}_{2}^{*}\right)}^{3}}\right)\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-{v}_{7}{w}_{7}\left({w}_{4}+{w}_{5}\right)\left(\frac{{\mu}_{2}}{{\left({x}_{1}^{*}+{x}_{2}^{*}\right)}^{2}}-\frac{{P}^{*}}{{x}_{1}^{*}+{x}_{2}^{*}}\right).\end{array}$ (65)

$b={\displaystyle \underset{k,i,j=1}{\overset{n}{\sum}}{v}_{k}{w}_{i}}\frac{{\partial}^{2}{f}_{k}}{\partial {x}_{i}\partial {P}^{*}}={v}_{7}\left({w}_{4}+{w}_{5}+{w}_{7}\right)\left(\frac{{x}_{7}^{*}}{{x}_{1}^{*}+{x}_{2}^{*}}\right).$ (66)

3.4. Vaccine Impact Analysis

Vaccine was believed to confer life-long immunity until 1990s. This was the norm as it was approximately correct for most available vaccine for infectious children diseases. But most vaccines used for combating adult infectious diseases today are defective and thus immunity conferred on the recipients wane with time. It is expected that the future Chikv vaccine will also be defective and hence the need to assess its effectiveness in ${R}_{C}$ a community. In this paper, the vaccine impact analysis is done by differentiating effective reproductive number with respect to the proportion p of susceptible individuals vaccinated at

equilibrium, according to [32] , $\left(p=\frac{{V}_{H}^{*}}{{N}_{H}^{*}}\right)$ i.e., $\frac{\partial {R}_{C}}{\partial p}=\frac{{R}_{C}\epsilon}{1\left(1-p\epsilon \right)}$ , i.e. since $0<\epsilon <1$ we have that $\frac{\partial {R}_{C}}{\partial p}<0$ , hence ${R}_{C}$ is a decreasing function of p. This means that a vaccination program with $p>0$ and $\epsilon >0$ at equilibrium, the future vaccine will have a positive impact. Besides, there exist a ${p}_{C}$ such that ${R}_{C}\left({p}_{C}\right)=1$ given by $\frac{1}{\epsilon}\left(1-\frac{1}{{R}_{C}}\right)$ and for vaccination of proportion of susceptible $p>{p}_{C}$ the number of new-cases reduces to zero faster than when $p<{p}_{C}$ .

4. Numerical Simulation

To further verify the analytical results in the model, the ode 45 code embedded in matlab was used to simulate some parameters of the model. Table 1 provided values of the parameters while initial values of the state variables were chosen arbitrarily. Figures 2(A)-(D) and Figures 3(A)-(D) are simulation of the various model compartments with time. Figure 4 is the simulation of some compartments

Figure 2. Plot of the various populations with parameters as in Table 1. (A) is the simulation of susceptible human against time, the plot shows that the susceptible human decreases with time due to the proportion that gets infected but slows down after some days, perhaps due to the vaccination and other control measures. (B) is the simulation of the vaccinated compartment. The plot shows a steady increase initially, but began to slope down after few days, this could be due to the fact that a proportion of the class are infectious as the vaccine is imperfect. (C) is the simulation of the exposed compartment with time, the plot shows a steady decline as members become infectious and progress to either the symptomatic or asymptomatic compartment. Finally (D) is the simulation of the symptomatic compartment with time. The plot shows a steady decline and tends to zero after about 20 days. This could be attributed to recovery from the infection.

Figure 3. Plot of the various populations with parameters as in Table 1. (A) is the simulation of the asymptomatic infected compartment with time, it shows a sharp increase at the onset of the epidemic, followed by a decline. (B) is the simulation of the recovered compartment with time, it shows a steady increase at the initial time, got to a peak and then remains a constant as time progresses. (C) is the simulation of susceptible mosquito compartment with time. It maintains a steady increase until perhaps due to short life cycle. (D) is the simulation of exposed mosquito compartment with time. The plot shows a steady decline with time as proportion progresses to infected compartment.

Figure 4. Simulation of the Human populations with varying values of β. (A) is the effect of the vertical transmission (β) on the susceptible compartment, while (B), (C) and (D) is the effect on same on the vaccinated, exposed, and symptomatic infected human compartment respectively. It is obvious from the plots that β has negligible effect in all the compartments and hence on the transmission of Chikungunya virus according to the model analysis and simulation.

with various values of the vertical transmission rate $\left(\beta \right)$ . Figure 5 is a contour plot of the effective basic reproduction number as a function of recruitment rate of susceptible mosquito $\left({\Lambda}_{M}\right)$ and vertical transmission rate $\left(\beta \right)$ while Figure 6 is the contour plot of effective basic reproductive number with varying values of vaccine efficacy $\left(\epsilon \right)$ and vaccinated proportion. Finally, Figure 7 is a simulation of the new cases of Chikungunya with different values of vaccine efficacy $\left(\epsilon \right)$ and vaccination rate $\left(\upsilon \right)$ . The figures and detailed caption are presented below.

Figure 5. Simulation of the chikv model displaying a contour graph of R_{C} as a function of recruitment rate of susceptible mosquito; and recruitment rate of infected mosquito (β) with parameter values as listed in Table 1.

Figure 6. Simulation of the chikv model displaying a contour graph of ( ${R}_{c}$ ) as a function of vaccinated human population and vaccine efficacy ( $\epsilon $ ); with parameter values as listed in Table 1.

Figure 7. Plot of the new-cases of Chikv model with varying values of vaccine efficacy and, proportion of vaccinated susceptible human population. For (A) $\nu =0.6,\epsilon =0.68,{R}_{C}=1.0076,{P}_{C}=0.63156$ For (B) $\nu =0.63,\epsilon =0.78,{R}_{C}=0.8727,{P}_{C}=0.63156$ .

5. Conclusion

A deterministic mathematical model for Chikungunya virus dynamics was developed using the standard incidence approach. The model assumed that the offspring of infected mosquito is infected at birth (vertical transmission) and also through blood meal from symptomatically and as-symptomatically infected human (horizontal transmission). For the subhuman population, only horizontal transmission was considered and the virus infection in human is assumed fatal, though with a very low rate. The disease free and endemic equilibrium was obtained and analyzed for both local and global asymptotically stability. The analysis shows that the model undergoes backward bifurcation when the effective basic reproductive number ${R}_{C}\le 1$ . Numerical simulation of the model shows that the effect of vertical transmission of the mosquito sub-population in the dynamics of the virus is negligible, even when the rate is high as shown in Figures 4(A)-(D). Further, the contour plot of the effective basic reproductive number ${R}_{C}$ with respect to the vaccine efficacy $\epsilon $ and the proportion of susceptible vaccinated (Figure 6) gave the rates at which the ${R}_{C}$ is above, below and equal to unity, this confirms that the use of imperfect vaccine will be effective. Figure 6 also reveals a linear relationship between the effective basic reproductive number and the two parameters in question unlike Figure 5. Also the graph of Chikungunya new case (Figure 7) shows a decrease in new cases with high vaccine efficacy $\epsilon $ and proportion of vaccinated susceptible $\nu $ . Hence buttressing the point made in Figure 6.

Appendix A

Castilo-Chaevz and Song [3]

Consider the following general system of ordinary differential equations with a parameter $\varphi $ .

$\frac{\text{d}x}{\text{d}t}=f\left(x,\varphi \right):{R}^{n}\times R\to {R}^{n}$ and $f\in {C}^{2}\left({R}^{n}\times R\right)$

where 0 is an equilibrium point of the system (that is, $f\left(0,\varphi \right)=0$ for all $\varphi $ ) and

(A1) $A={D}_{x}f\left(0,0\right)=\left(\frac{\partial {f}_{i}}{d{x}_{j}}\left(0,0\right)\right)$ is the linearization matrix of the system 2.10 around the equilibrium 0 with $\varphi $ evaluated at 0;

(A2) Zero is a simple eigenvalues of A and other eigenvalues of A have negative real parts;

(A3) Matrix A has a right eigenvector w and left eigenvector v (each corresponding to zero eigenvalues).

Let ${f}_{k}$ be the kth component of f and

To do this we need the values of a and b given below:

$a={\displaystyle \underset{k,i,j=1}{\overset{n}{\sum}}{v}_{k}{w}_{i}{w}_{j}\frac{{\partial}^{2}{f}_{k}}{\partial {x}_{i}\partial {x}_{j}}\left(0,0\right)}$

$b={\displaystyle \underset{k,i,=1}{\overset{n}{\sum}}{v}_{k}{w}_{i}\frac{{\partial}^{2}{f}_{k}}{\partial {x}_{i}\partial {x}_{\varphi}}\left(0,0\right)}$

then, the local dynamics of the system around equilibrium point 0 is totally determined by the signs of a and b, particularly,

1) $a>0,b>0$ , when $\varphi <0$ with $\left|\varphi \right|\ll 1$ , 0 is locally asymptotically stable and there exists a positive unstable equilibrium; when $0<\varphi \ll 1$ , 0 is unstable and there exists a negative, locally asymptotically stable equilibrium;

2) $a<0,b<0$ , when $\varphi <0$ with $\left|\varphi \right|\ll 1$ , 0 is unstable; when $0<\varphi \ll 1$ , 0 is locally asymptotically stable equilibrium and there exists a positive unstable equilibrium;

3) $a<0,b>0$ , when $\varphi $ changes from negative to positive, 0 changes its stability from stable to unstable. Correspondingly a negative unstable equilibrium becomes positive and locally asymptotically stable.

References

[1] Robinson, M.C. (1955) An Epidemic of Virus Disease in Southern Province, Tanganyika Territory, in 1952-1953. I. Clinical Features. Transactions of the Royal Society of Tropical Medicine and Hygiene, 49, 28-32.

https://doi.org/10.1016/0035-9203(55)90080-8

[2] Vazeille, M., Moutailler, S., Coudrier, D., Rousseaux, C., Khun, H., et al. (2007) Two Chikungunya Isolates from the Outbreak of La Reunion (Indian Ocean) Exhibit Different Patterns of Infection in the Mosquito, Aedes albopictus. PLoS ONE, 2, e1168.

https://doi.org/10.1371/journal.pone.0001168

[3] Moore, D.L., Readdy, S., Akinkugbe, F.M., Lee, V.H., David-West, T.S., Causey, O.R. and Carez, D.E. (1974) An Epidemic of Chikungunya Fever at Ibadan, Nigeria, 1969. Annals of Tropical Medicine & Parasitology, 68, 59-68.

https://doi.org/10.1080/00034983.1974.11686925

[4] World Health Organisation (2017).

http://www.who.int/mediacentre/factsheets/fs327/en/

[5] Chandrakant, L., and Pradhan, S.K. (2006) Emergence of Chikungunya Virus in India Subcontinent after 32 Years, a Review. Journal of Vector Borne Diseases, 43, 151-160.

[6] Pialoux, G., Gauze‘re, B.A., Jaureguiberry, S. and Strobel, M. (2007) Chikungunya, an Epidemic Arbovirosis. The Lancet, 7, 319-327.

https://doi.org/10.1016/S1473-3099(07)70107-X

[7] Peterson, L. R., and Powers, A.M. (2016) Chikungunya: Epidemiology. F1000 Research, 5, 82-89.

[8] Reiter, P. (2007) Oviposition, Dispersed and Survival in Aedes aegypti: Implication for the Efficiency of Control Strategies. Vector-Borne Zoonetic Disease, 7, 261-273.

https://doi.org/10.1089/vbz.2006.0630

[9] Bonizzoni, M., Gezperi, G., Chen, X. and James, A.A. (2013) The Invasive Mosquito Species Aedes albopictus: Current Knowledge and Future Perspective. Trends in Parasitology, 29, 460-468.

https://doi.org/10.1016/j.pt.2013.07.003

[10] Rao, G., Khan, Y.Z. and Chitnis D.S. (2008) Chikungunya Infection in Neonates. Indian Pediatics, 45, 240-244.

[11] Anderson, R.M. and May, R.M. (1991) Infectious Diseases of Humans: Dynamics and Control. Oxford University Press, Oxford.

[12] Hethcote, H.W. (2000) The Mathematics of Infectious Diseases. SIAM Review, 42, 599-653.

https://doi.org/10.1137/S0036144500371907

[13] Diekman, O. and Heesterbeck, J.A.P. (2000) Mathematical Epidemiology of Infectious Diseases. Wiley-Blackwell, New York.

[14] Dumont, Y., Chiroleu, F. and Domerg, C. (2008) On a Temporal Model for the Chikungunya Disease: Modeling, Theory and Numerics. Mathematical Biosciences, 213, 80-91.

https://doi.org/10.1016/j.mbs.2008.02.008

[15] Moulay, D., Aziz-Alaoui, M.A. and Cadivel, M. (2011) The Chikungunya Disease: Modeling, Vector and Transmission Global Dynamics. Mathematical Biosciences, 229, 50-63.

https://doi.org/10.1016/j.mbs.2010.10.008

[16] Ruiz-Moreno, D., Vargas, I.S., Olson, K.E. and Harrington, L.C. (2012) Modeling Dynamic Introduction of Chikungunya Virus in the United States. PLoS Neglected Tropical Disease, 6, e1918.

https://doi.org/10.1371/journal.pntd.0001918

[17] Pongsumpun, P. and Sangsawangl, S. (2013) Local Stability Analysis for Age Structural Model of Chikungunya Disease. Journal of Basic and Applied Scientific Research, 3, 302-312.

[18] Yakob, L. and Clements, A.C.A. (2013) A Mathematical Model of Chikungunya Dynamics and Control: The Major Epidemic on Réunion Island. PLoS ONE, 8, e57448.

https://doi.org/10.1371/journal.pone.0057448

[19] Manore, C.A., Hickmann, K.S., Xu, S., Wearing, H.J. and Hyman, J.M. (2014) Comparing Dengue and Chikungunya Emergence and Endemic Transmission in A. aegypti and A. albopictus. Journal of Theoretical Biology, 356, 174-191.

https://doi.org/10.1016/j.jtbi.2014.04.033

[20] Chitnis, N., Hyman, J.M. and Cushing, J.M. (2008) Determining Important Parameters in the Spread of Malaria through the Sensitivity Analysis of a Mathematical Model. Bulletin of Mathematical Biology, 70, 1272-1296.

https://doi.org/10.1007/s11538-008-9299-0

[21] Dumont, Y. and Chiroleu, F. (2010) Vector Control for the Chikungunya Disease. Mathematical Biosciences and Engineering, 7, 315-348.

https://doi.org/10.3934/mbe.2010.7.313

[22] Nur Aida, H., Abu Hassan, A., Nurita, A.T., Che Salmah, M.R. and Norasmah, B. (2008) Population Analysis of Aedes albopictus (Skuse) (Diptera: Culicidae) under Uncontrolled Laboratory Conditions. Tropical Biomedicine, 25, 117-125.

[23] Lahariya, C. and Pradhan, S.K. (2006) Emergence of Chikungunya Virus in Indian Subcontinent after 32 Years: A Review. Journal of Vector Borne Diseases, 43, 151-160.

[24] Delatte, H., Gimonneau, G., Triboire, A. and Fontenille, D. (2009) Influence of Temperature on Immature Development, Survival, Longevity, Fecundity, and Gonotrophic Cycles of Aedes albopictus, Vector of Chikungunya and Dengue in the Indian Ocean. Journal of Medical Entomology, 46, 33-41.

https://doi.org/10.1603/033.046.0105

[25] Hethcote, H.W. (1978) An Immunization Model for a Heterogeneous Population. Theoretical Population Biology, 14, 338-349.

https://doi.org/10.1016/0040-5809(78)90011-4

[26] Institute of Epidemiology, Disease Control and Research (2017) Chikungunya News Letter.

http://www.iedcr.gov.bd/index.php/chikungunya/chikungunya-newsletter

[27] Lakshmikantham, V., Leela, S. and Martynyuk, A.A. (1989) Stability Analysis of Non-Linear Systems. Marcel Dekker, New York.

[28] Castillo-Chavez, C., Blower, S., van den Driessche, P., Kirschner, D. and Yakubu, A.A. (2002) Mathematical Approaches for Emerging and Reemerging Infectious Diseases. Springer-Verlag, New York.

[29] Powell, D.R., Faie, J., Leclaire, R.J., Moore, L.M. and Thompson, D. (2005) Sensitivity Analysis of an Infectious Disease Model. International System Dynamic Conference, Boston, MA.

[30] Castillo-Chavez, C. and Song, B. (2004) Dynamical Models of Tuberculosis and Their Applications. Mathematical Biosciences & Engineering, 1, 361-404.

https://doi.org/10.3934/mbe.2004.1.361

[31] Agosto, F.B., Easley, S., Freeman, K. and Thomas, M. (2016) Mathematical Model of Three Age Structured Transmission Dynamics of Chikungunya Virus. Computational and Mathematical Methods in Medicine, 2016, Article ID: 4320514.

https://doi.org/10.1155/2016/4320514

[32] Garba, M.S., Gumel, A.B. and Bakar, A.M.R. (2008) Backward Bifurcation in Dengue Transmission Dynamics. Mathematical Bioscience, 215, 11-25.

https://doi.org/10.1016/j.mbs.2008.05.002