CUBO, A Mathematical Journal Vol.22, N◦02, (177–201). August 2020 http://dx.doi.org/10.4067/S0719-06462020000200177 Received: 28 February, 2020 | Accepted: 15 June, 2020 Mathematical Modeling of Chikungunya Dynamics: Stability and Simulation Ruchi Arora, Dharmendra Kumar, Ishita Jhamb and Avina Kaur Narang Department of Mathematics, SGTB Khalsa College, University of Delhi, Delhi-110007, India ruchi@sgtbkhalsa.du.ac.in, dharmendrakumar@sgtbkhalsa.du.ac.in, ishita.jhamb@gmail.com, kaur.avina45@gmail.com ABSTRACT Infection due to Chikungunya virus (CHIKV) has a substantially prolonged recupera- tion period that is a long period between the stage of infection and recovery. However, so far in the existing models (SIR and SEIR), this period has not been given due atten- tion. Hence for this disease, we have modified the existing SEIR model by introducing a new section of human population which is in the recuperation stage or in other words the human population that is no more showing acute symptoms but is yet to attain complete recovery. A mathematical model is formulated and studied by means of exis- tence and stability of its disease free equilibrium (DFE) and endemic equilibrium (EE) points in terms of the associated basic reproduction number (R0). RESUMEN La infección debida al virus Chikungunya (CHIKV) tiene un peŕıodo de recuperación sustancialmente prolongado, que es un peŕıodo largo entre la etapa de infección y recuperación. Sin embargo, hasta ahora en los modelos existentes (SIR y SEIR), este peŕıodo no ha recibido suficiente atención. Por tanto, para esta enfermedad, hemos modificado el modelo SEIR existente introduciendo una nueva sección de población humana que está en la etapa de recuperación o, en otras palabras, la población humana que ya no muestra śıntomas agudos pero todav́ıa no se recupera completamente. Se formula y estudia un modelo matemático a través de la existencia y estabilidad de su equilibrio libre de enfermedad (DFE) y puntos de equilibrio endémico (EE) en términos del número de reproducción básico asociado (R0). Keywords and Phrases: Equilibrium point, disease free equilibrium, endemic equilibrium, re- production number, local stability, global stability. 2020 AMS Mathematics Subject Classification: 92B05, 93A30, 93C15 c©2020 by the author. This open access article is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License. http://dx.doi.org/10.4067/S0719-06462020000200177 178 Ruchi Arora, Dharmendra Kumar, Ishita Jhamb and Avina Kaur CUBO 22, 2 (2020) 1 Introduction In recent past, the study of vector borne diseases has gained considerable attention and mathemat- ics have become a useful tool for such studies. Several temporal deterministic models have been proposed for diseases like dengue, malaria, chikungunya etc. Chikungunya is a disease caused by the chikungunya virus, an RNA genome which is a member of the Alphavirus genus in the family of Togaviridae. It is a mosquito borne viral disease which is transmitted to humans through Aedes aegypti mosquito bite [1]. In 1952, chikungunya was first confirmed as the cause of an epidemic of dengue like illness on the Comoros islands located on the eastern coast of northern Mozam- bique [2]. Since its discovery, numerous CHIKV outbreaks with irregular intervals of 2-20 years have affected Asian, African, European and American countries. In Thailand, the first report of chikungunya infection occurred in Bangkok in 1958 [3]. In India, the virus emerged in parts of Vellore, Calcutta and Maharashtra in the early 1960’s [4]. The virus continued to spread in Sri Lanka in 1969 and many countries of Southeast Asia such as Myanmar, Indonesia and Vietnam [4]. Later, some irregular cases of chikungunya fever were also seen in many provinces of Thailand in the period from 1976 to 1995 [3]. From 1999 to 2000, the reemergence of chikungunya occurred in Democratic Republic of Congo [2], 13,500 cases were reported in Lamu, Kenya in 2004 [5]. In the years 2005-2007, there occurred an outbreak in Reunion islands in the Indian Ocean. In 2007, 197 cases were reported in Europe due to chikungunya [1]. The outbreak mutated to facilitate the disease transmission by Aedes albopictus from the tiger mosquito family. It was a mutation in one of the viral envelope genes which allowed the virus to be present in the mosquito saliva only two days after the infection and seven days in Aedes aegypti mosquitoes. The results indicated that the areas where the tiger mosquitoes are present could have a greater risk of outbreak. After an effective bite from a mosquito infected with CHIKV, the incubation period (i.e., the time elapsed between exposure to pathogenic organism and when symptoms and signs are first apparent) usually lasts for 3-7 days with fever as the most prominent symptom. The symptoms of chikungunya fever differ from the normal fever as they are accompanied with acute joint pains. Other common symptoms are nausea, rashes, headache and fatigue. Some cases may result in neu- rological, retinal and carpological complications as well, which makes it difficult for older people to recover as against young people. In some instances, people live with joint pains for years which indicates that the recuperation period can last for a long time. The symptoms of chikungunya are generally mild and the disease may sometimes be misdiagnosed with Zika and Dengue due to similarity in symptoms. There have been very few cases where chikungunya resulted in death and mostly infected individuals are expected to make full recovery with lifelong immunity. As such, there is no preventive vaccine or cure for chikungunya. One can only manage the symptoms by taking medications for temporary relief. To prevent the spread of disease, breeding sites for the mosquitoes should be checked. Using mosquito repellents and wearing long sleeve clothes and full CUBO 22, 2 (2020) Mathematical Modeling of Chikungunya Dynamics . . . 179 pants can help in preventing mosquito bite. For more such information one may refer to [1]. Increasing globalization and factors contributing to climate change brought about a sudden expansion of mosquito breeding sites. This makes it necessary to improve the vector control tech- niques and to identify the indexes that monitor thresholds for such programs. Through the 20th century, mathematical modeling has been extensively used to study epidemic diseases. Futher- more, this branch of mathematics is also being used to devise optimal control strategies for various infectious diseases. Like M. Barro et al. [6] introduced an optimal control for a SIR model governed by an ODE system with time delay. And, O. K. Oare [7] considered and analyzed a deterministic multipatch hepatitis C virus model for it. In context of infection due to chikungunya virus, Y. Dumont et al. [8] proposed a model associated with the time course of the first epidemic of chikungunya in several cities of Reunion Island. A model describing the mosquito population dynamics and the virus transmission to human population was discussed by D. Moulay et al. [9]. Although simplistic, L. Yacob et al. [10] gave a model which provided a close approximation of the peak incidence of the outbreak and the final epidemic size. S. Naowarat and I. M. Tang [11] studied the model taking into consideration the presence of two species of Aedes mosquito (Aedes aegypti and Aedes albopictus). D. H. Palacio and J. Ospina [12] derived measures of disease control, by means of three scenarios, namely a single vector, two vectors, and two vectors and human and non-human reservoirs. It also showed the need to periodically evaluate the effectiveness of vector control measures. F. B. Agusto et al. [13] described the chikungunya model of three age structured transmission dynamics by considering juvenile, adult and senior population, where the dynamics of shift in individuals from one stage to another was studied. In this paper, we introduce a deterministic model to study the dynamics and transmission of chikungunya virus by considering a very significant section from the class of infected individuals. Usually, the existing models focus on the SIR or the SEIR human population model and SEI mosquito population model. Since the period from the infected stage to the complete recovery stage is quite long for this disease, so it becomes significant to study that particular class of human population which has recovered from acute symptoms of the disease but is yet to attain full recovery. Though the class no longer shows the immediate symptoms like fever, rashes, nausea etc. but at the same time they are bearing the latent and the passive effects of the disease like joint pains, fatigue, headache etc. Generally such ailments continue for a prolonged period which may vary from individual to individual. But as long as the patient is suffering from these ailments, he or she cannot be declared as fully recovered [14]. Focussing on this category of patients, we introduce a new compartment between compartments of the infected and the recovered human population within the existing SEIR model. We refer to it as the recuperation compartment and denote it by R′. So, in this paper our aim is to study, analyse and investigate in detail the model showing the interaction between the human population divided into five compartments resulting 180 Ruchi Arora, Dharmendra Kumar, Ishita Jhamb and Avina Kaur CUBO 22, 2 (2020) into a SEIR′R model and the mosquito population into the traditional three compartments which we denote by XYZ model. The paper is divided as follows: Section 2 deals with the formulation of the model, section 3 analyses its feasibility, section 4 determines the disease free equilibrium (DFE) and establishes its local and global stability , section 5 deals with the existence of endemic equilibrium (EE) and its local stability. Also by means of simulation of the formulated model, we provide a visualization to the dynamics of this disease, in section 6. Finally related to our model, some conclusions are stated. 2 Model Formulation In this section, an epidemic model is formulated for chikungunya disease. Let NH represent the total human population which is further subdivided into five categories; susceptibles (S), humans exposed to infection (E), infected humans (I), population in recuperation phase (R′) and finally the population that has attained complete recovery (R). So, the traditional SEIR epidemic model has been modified to a more relevant and practically applicable SEIR′R model. Hence in this case, at any time t NH(t) = S(t) + E(t) + I(t) + R ′(t) + R(t). (2.1) Let NM represent the total mosquito population which is further subdivided into 3 parts; suscep- tible mosquitoes (X), mosquitoes exposed to infection (Y) and infectious mosquitoes (Z). So the total mosquito population is NM (t) = X(t) + Y (t) + Z(t). For human population, let µ be the constant birth rate and ζ be the natural death rate. Then the rate of change of susceptible human population is given by dS dt = µ − λHS − ζS, (2.2) where λH = βBHZ NH . BH is the transmission probability per contact for susceptible humans (S) and β is the mosquito biting rate for transfer of infection from infectious mosquito class (Z) to susceptible human population (S). As only the susceptible human population out of the whole population is prone to get infection, thereby we divide the expression by NH. The rate of change of exposed human population is given by dE dt = λHS − αE − ζE, (2.3) where α is the rate of progression from exposed (E) to infected (I) human population. Here the inflow rate is λH and outflow rate is α + ζ. Similarly, the rate of change of infected human population is CUBO 22, 2 (2020) Mathematical Modeling of Chikungunya Dynamics . . . 181 dI dt = αE − γI − (ζ + ζ1)I, (2.4) where ζ1 is death rate due to infection and γ is progression rate of infected (I) to recuperated (R′) human population. Now, rate of change of human population in recuperation phase is dR′ dt = γI − λR′ − (ζ + ζ2)R ′, (2.5) where ζ2 is the death rate of humans in recuperated phase due to virus and λ is the rate of progression from recuperation (R′) to the recovery phase (R). Finally, rate of change of recovered human population is, dR dt = λR′ − ζR. (2.6) Again for the mosquito population, let ρ be the constant birth rate and κ be the natural death rate, then the rate of susceptible mosquito population is given by dX dt = ρ − λMX − κX, (2.7) where λM = νBM (I + R ′) NH . BM is the transmission probability per contact for susceptible mosquito population (X) and ν is the mosquito biting rate for transfer of infection from infected (I) or recuperated (R′) human population to susceptible mosquito population (X). Again there occurs division by NH because infection can be transfered to mosquitoes only by a certain fraction of human population. Now, the rate of change of exposed mosquito population is given by dY dt = λMX − ψY − κY, (2.8) where ψ is the progression rate from exposed (Y) to infectious (Z) mosquito population. Here the inflow rate is λM and outflow rate is ψ + κ. Similarly, the rate of change of mosquito population carrying infection is dZ dt = ψY − κZ. (2.9) Compiling the above discussion, we get the eight dimensional system of nonlinear ordinary differ- ential equations that forms our Chikungunya Model (CM). The parameters and the variables used in the model (CM) are described in Table 1. To get a clear view of the inter relationships between various compartments in discussion, one may refer to Figure 1 which shows the schematic flow 182 Ruchi Arora, Dharmendra Kumar, Ishita Jhamb and Avina Kaur CUBO 22, 2 (2020) diagram of the model. The model (CM) is as follows: (CM) dS dt = µ − βBHZS NH − ζS, dE dt = βBHZS NH − αE − ζE, dI dt = αE − γI − (ζ + ζ1)I, dR′ dt = γI − λR′ − (ζ + ζ2)R ′, dR dt = λR′ − ζR, dX dt = ρ − νBM (I + R ′)X NH − κX, dY dt = νBM (I + R ′)X NH − ψY − κY, dZ dt = ψY − κZ. Figure 1: Schematic diagram of Chikungunya Model (CM) CUBO 22, 2 (2020) Mathematical Modeling of Chikungunya Dynamics . . . 183 Table 1: Description of variables and parameters used in model (CM) Variables Description S Susceptible human population. E Exposed human population. (Population that is infected but yet to show symptoms). I Infected Human population showing symptoms. R′ Human population in recuperation phase. R Fully recovered human population. X Susceptible mosquito population. Y Exposed mosquito population. (carrying infection but not yet capable to spread it). Z Infectious mosquito population spreading the disease. Parameters Description µ Human birth rate. β Mosquito biting rate for transfer of infection from infectious mosquito class (Z) to susceptible human population (S). α Progression rate of exposed to infected human population. γ Progression rate of infected to recuperated human population. λ Progression rate of recuperated to fully recovered human population. ρ Mosquito birth rate. ν Mosquito biting rate for transfer of infection from infected human population(I) or population under recuperation phase (R′) to susceptible mosquito population (X). ψ Progression rate from exposed to infectious mosquito population. ζ Natural death rate for human population. ζ1 Human death rate in infected stage due to viral infection. ζ2 Human death rate due to infection under recovery phase. κ Natural death rate for mosquito population. BH Transmission probability per contact in susceptible humans. BM Transmission probability per contact in susceptible mosquitoes. NH Total human population, i.e. S+E+I+R ′+R. 184 Ruchi Arora, Dharmendra Kumar, Ishita Jhamb and Avina Kaur CUBO 22, 2 (2020) Table 2: Range of Parameters for the model (CM) Parameters Range References µ 400 × 1 15 × 365 - 400 × 1 12 × 365 [15, 16] β 0.19 - 0.39 [15, 17] α 1 4 - 1 2 [4, 15, 18, 19, 20, 21] γ 1 4 - 1 2 Estimated [14] λ 1 8 - 1 4 Estimated [14] ρ 500 × 0.015 - 500 × 0.33 [15, 16, 22, 23] ν 0.19 - 0.39 [15, 17] ψ 1 6 - 1 2 [9, 18, 20, 24] ζ 1 60 × 365 - 1 18 × 365 [13] ζ1 1 105 - 1 104 [25] ζ2 1 106 - 1 105 [25] κ 1 42 - 1 14 [9, 18, 19, 20, 21] BH 0.001 - 0.54 [8, 15, 26, 18, 27] BM 0.005 - 0.35 [8, 26, 27, 28, 29] CUBO 22, 2 (2020) Mathematical Modeling of Chikungunya Dynamics . . . 185 Table 3: Values of Parameters for Simulation Parameters R0 < 1 R0 > 1 µ 400 × 1 15 × 365 400 × 1 15 × 365 β 0.25 0.30 α 1 3 1 4 γ 1 3 1 4 λ 1 7 1 8 ρ 500 × 0.1675 500 × 0.2 ν 0.25 0.30 ψ 1 3.5 1 4 ζ 1 40 × 365 1 30 × 365 ζ1 1 104 1 105 ζ2 1 105 1 106 κ 1 14 1 30 BH 0.24 0.30 BM 0.24 0.30 186 Ruchi Arora, Dharmendra Kumar, Ishita Jhamb and Avina Kaur CUBO 22, 2 (2020) 3 Preliminary Results 3.1 Positivity of Solutions In order to establish the epidemiological meaningfullness [13], we prove the non negativity of the state variables for the formulated model at all t > 0. Theorem 3.1: The solution M(t) = (S,E,I,R′,R,X,Y,Z) of model (CM) with M(0) ≥ 0, is non negative for all t > 0. Moreover, lim t→∞ sup NH(t) = µ ζ and lim t→∞ sup NM(t) = ρ κ where NH(t) = S(t) + E(t) + I(t) + R ′(t) + R(t) and NM (t) = X(t) + Y (t) + Z(t). Proof: Let t1 = sup {t > 0 : M(t) > 0}. Clearly t1 > 0. Consider the first equation of the model (CM), dS dt = µ − βBHSZ NH − ζS. Solving the differential equation we have, d dt { S(t) exp [( ∫ t1 0 βBHZ(τ) NH(τ) dτ + ζt )]} = µexp [( ∫ t1 0 βBHZ(τ) NH(τ) dτ + ζt )] =⇒ S(t1) exp [( ∫ t1 0 βBHZ(τ) NH(τ) dτ + ζt1 )] − S(0) = ∫ t1 0 µexp [( ∫ u 0 βBHZ(τ) NH(τ) dτ + ζu )] du. Furthermore, S(t1) =S(0) exp [( − ∫ t1 0 βBHZ(τ) NH(τ) dτ + ζt1 )] + exp [( − ∫ t1 0 βBHZ(τ) NH(τ) dτ + ζt1 )] ∫ t1 0 µexp [( ∫ u 0 βBHZ(τ) NH(τ) dτ + ζu )] du > 0. Similarly, the non negativity can be shown for all the state variables, i.e., M(t1) > 0 and therefore M(t) > 0 for all t > 0. In fact, we now have, 0 < S(t) ≤ NH(t), 0 < E(t) ≤ NH(t), 0 < I(t) ≤ NH(t), 0 < R ′(t) ≤ NH(t), 0 < R(t) ≤ NH(t); 0 < X(t) ≤ NM(t), 0 < Y (t) ≤ NM(t), 0 < Z(t) ≤ NM (t). As the total human population is given by NH(t) = S(t) + E(t) + I(t) + R ′(t) + R(t), the rate of change of human population with respect to time is given by dNH dt = µ − ζ(S + E + I + R′ + R) − ζ1I − ζ2R ′ = µ − ζNH − ζ1I − ζ2R ′ ≤ µ − ζNH. (3.1) CUBO 22, 2 (2020) Mathematical Modeling of Chikungunya Dynamics . . . 187 Now for NM(t) = X(t) + Y (t) + Z(t), dNM dt ≤ ρ − κNM. Let N = µ ζ . As t → ∞, the disease will disappear. Therefore, lim t→∞ sup I(t) = 0 and lim t→∞ sup R′(t) = 0. Now, dNH dt = µ − ζNH this implies NH(t) = µ ζ + ( NH(0) − µ ζ ) e−ζt, which further implies lim t→∞ NH(t) = µ ζ = N. This follows that 0 < lim t→∞ supNH(t) ≤ N = µ ζ if lim t→∞ sup I(t) = 0 and lim t→∞ sup R′(t) = 0. And if NH > N = µ ζ then from (3.1), dNH dt < 0. Similarly, it can be seen that 0 < lim t→∞ supNM(t) ≤ ρ κ . 3.2 Invariant Region Consider ℜ = ℜH × ℜM ⊂ R 5 + × R 3 +, where ℜH = { S,E,I,R′,R : NH(t) ≤ µ ζ } , ℜM = { X,Y,Z : NM(t) ≤ ρ κ } . Now, we establish the positive invariance [13], of the region ℜ associated to the model (CM). That is, we show that solutions in ℜ remain in ℜ for all t > 0 . Theorem 3.2: The region ℜ ⊂ R8+ is positively invariant for the model (CM), with non-negative initial conditions in R8+. Proof : As seen in Theorem 3.1, dNH dt ≤ µ − ζNH and dNM dt ≤ ρ − κNM. By using standard comparison theorem [30], it can be seen that, NH(t) ≤ µ ζ = N. So, clearly every solution in ℜH remains in ℜH for all t > 0. Similar is the case for every solution of ℜM . Hence, the region ℜ is positively invariant and contains all solutions of R8+ for model (CM). In the following sections, we show the existence and stability of the disease free equilibrium (DFE) and endemic equilibrium (EE) for the model (CM). 4 Disease Free Equilibrium (DFE) In this section, we find a unique disease free equilibrium (DFE) for the model (CM) and then analyse its stability. 188 Ruchi Arora, Dharmendra Kumar, Ishita Jhamb and Avina Kaur CUBO 22, 2 (2020) 4.1 Existence of Equilibrium To determine the disease free equilibrium (DFE) of the model, we consider the sections of pop- ulations that are free from disease and put their time derivatives equal to zero. Let DFE be denoted by Ed = (S ∗,E∗,I∗,R ′ ∗,R∗,X∗,Y ∗,Z∗). As sections of susceptible and recovered hu- mans as well as susceptible mosquitoes are the only sections free from disease therefore Ed = (S∗,0,0,0,R∗,X∗,0,0). Solving the differential equations of the model (CM), DFE is obtained as Ed = ( µ ζ ,0,0,0,0, ρ κ ,0,0 ) . 4.2 Reproduction Number Let the basic reproduction number be denoted by R0, which is defined as the expected number of secondary cases produced by a single (typical) infection in a population that is completely disease free. To find the threshold quantity R0 [31, 32], we consider the next generation matrix G, which comprises of two matrices F and V −1, where F = dFi(x0) dxj and V = dVi(x0) dxj for 1 ≤ i,j ≤ 5. Here, Fi represents the new infection, whereas Vi corresponds to the transfers of infection from one compartment to another. Let x0 be the disease free equilibrium state. Hence, the reproduction number is the largest eigen value of the next generation matrix G (defined as the product of matrices F and V −1), that is the largest eigen value of the matrix, G = FV −1. Corresponding to the model (CM), F =       βBHSZ NH 0 0 νBM (I+R ′)X NH 0       and V =       αE + ζE −αE + γI + (ζ + ζ1)I −γI + λR′ + (ζ + ζ2)R ′ ψY + κY −ψY + κZ       . Next, we find the Jacobian F and V of the matrices F and V respectively and the eigen values of the matrix G = FV −1, gives the reproduction number as R0 = √ ρνψζαβBHBM (λ + γ + ζ + ζ2) κ √ µ(ψ + κ)(ζ + α)(ζ + γ + ζ1)(ζ + λ + ζ2) . CUBO 22, 2 (2020) Mathematical Modeling of Chikungunya Dynamics . . . 189 4.3 Local Stability Theorem 4.1 : The DFE of the chikungunya model (CM) is locally asymptotically stable, if R0 < 1 and unstable if R0 > 1, where R0 is the associated reproduction number. Proof : We consider the system of non linear differential equations, corresponding to the model (CM) to evaluate its Jacobian matrix. Let JD denote the Jacobian of the system at DFE that is, JD =             −ζ 0 0 0 0 0 0 −BHβ 0 −α − ζ 0 0 0 0 0 BHβ 0 α −γ − ζ − ζ1 0 0 0 0 0 0 0 γ −λ − ζ − ζ2 0 0 0 0 0 0 0 λ −ζ 0 0 0 0 0 −νBM ρζ κµ −νBM ρζ κµ 0 −κ 0 0 0 0 νBM ρζ κµ νBM ρζ κµ 0 0 −ψ − κ 0 0 0 0 0 0 0 ψ −κ             Clearly, the trace of the matrix JD is negative and determinant of matrix JD [33, 34], is given by det(JD) = −ζ 2[κ2µ(ψ + κ)(ζ(ζ + α + γ) + αγ + ζζ1 + ζ1α)(−ζ − λ − ζ2)] + ρνζψαβBHBM (ζ + λ + γ + ζ2) µ . For R0 < 1, we have √ ρνψζαβBHBM (ζ + γ + λ + ζ2) < κ √ µ(ψ + κ)(ζ + α)(ζ + λ + ζ2)(ζ + γ + ζ1). Therefore, κ 2 µ(ψ + κ)(ζ + λ + ζ2)(ζ(ζ + α + γ) + αγ + ζζ1 + ζ1α) − ψ[ρνζαβBHBM(ζ + λ + γ + ζ2)] > 0 or det(JD) > 0. Hence, DFE is locally asymptotically stable if R0 < 1 . 4.4 Global Stability Consider the feasible region ℜ1 = {D ∈ ℜ : S ≤ S ∗,X ≤ X∗} where D = (S,E,I,R′,R,X,Y,Z), S∗ and X∗ are the components of DFE (Ed). Lemma 4.1: The region ℜ1 is positively invariant for the model (CM). 190 Ruchi Arora, Dharmendra Kumar, Ishita Jhamb and Avina Kaur CUBO 22, 2 (2020) Proof: From the first equation of the model (CM), dS dt = µ − βBHZS NH − ζS ≤ µ − ζS ≤ ζ ( µ ζ − S ) ≤ ζ(S∗ − S) S ≤ S∗ + (S(0) − S∗)e−ζt Thus, if S∗ = µ ζ for all t ≥ 0 and S(0) ≤ S∗ , then S ≤ S∗ for all t ≥ 0. Similarly, for dX dt = ρ − νBM (I + R ′)X NH − κX ≤ ρ − κX ≤ κ(X∗ − X) X ≤ X∗ + (X(0) − X∗)e−κt Thus, if X∗ = ρ κ for all t ≥ 0 and X(0) ≤ X∗, then X ≤ X∗ for all t ≥ 0. Hence, it has been shown that the region ℜ1 is positively invariant and attracts all solutions in ℜ 8 + for the model (CM). Now in order to establish the global asymptotic stability of DFE [35], we rewrite the model (CM) as [ dTU dt = F(TU,TI) dTI dt = G(TU,TI), G(TU,0) = 0 ] (RM) where TU = (S,R,X) ∈ R 3 and TI = (E,I,R ′,Y,Z) ∈ R5. Let E∗D = (T ∗ U,0) be DFE of (RM) where T ∗ U = ( µ ζ ,0, ρ κ ) . We now state the following two condi- tions which must be satisfied to guarantee global asymptotic stability: (H1) For dTU dt = F(TU,0), T ∗ U is globally asymptotically stable. (H2) G(TU,TI) = ATI − Ĝ(TU,TI), Ĝ(TU,TI) ≥ 0, (TU,TI) ∈ ℜ where A = ∂G(T ∗U,0) ∂TI is an M-matrix which by definition has the off diagonal elements non-negative. Theorem 4.2: The fixed point E∗D = (T ∗ U,0) is globally asymptotic stable (g.a.s) equilibrium of (RM) provided that R0 < 1 and that assumptions (H1) and (H2) are satisfied. Proof: For the system (RM), dTU dt = F(TU,0) =    µ − ζS 0 ρ − κX    CUBO 22, 2 (2020) Mathematical Modeling of Chikungunya Dynamics . . . 191 We solve the above linear differential system to get the S(t) = µ ζ + S∗(0)e−µt, R(t) = 0 and X(t) = ρ κ + X∗(0)e−κt which implies S(t) → µ ζ , R(t) → 0 and X(t) → ρ κ as t → ∞. Therefore, disease free point T ∗U is a globally asymptotic stable (g.a.s) equilibrium of dTU dt = F(TU,0). Hence (H1) holds. Clearly it can be seen that G(TU,TI) =       βBHZS NH − αE − ζE αE − γI − (ζ + ζ1)I γI − λR′ − (ζ + ζ2)R ′ νBM (I+R ′)X NH − ψY − κY ψY − κZ       Also from (H2) G(TU,TI) = ATI − Ĝ(TU,TI), where A = ∂G(T ∗U,0) ∂TI =       −α − ζ 0 0 0 βBH α −γ − ζ − ζ1 0 0 0 0 γ −λ − ζ − ζ2 0 0 0 νBM ρζ κµ νBM ρζ κµ −ψ − κ 0 0 0 0 ψ −κ       . Therefore, ∂G(T ∗U,0) ∂TI TI =       −αE − ζE − βBHZ αE − γI − (ζ + ζ1)I γI − λR′ − (ζ + ζ2)R ′ νBM aζ κµ (I + R′) − (ψ + κ)Y ψY − κZ       . In view of (H2), Ĝ(TU,TI) = ∂G(T ∗ U ,0) ∂TI TI − G(TU,TI) which gives Ĝ(TU,TI) =        βBHZ(1 − S NH ) 0 0 νBM (I + R ′) [ ρζ µκ − X NH ] 0        . Clearly, βBHZ(1 − S NH ) ≥ 0 as S NH < 1. Also, κX ρ ≤ ζNH µ or ρζ µκ ≥ X NH ⇒ X∗ S∗ ≥ X NH and from Lemma 4.1 we know X∗ ≥ X and N∗H = S ∗ ≥ NH, which implies νBM (I + R ′) [ ρζ µκ − X NH ] ≥ 0. Therefore, (H2) holds true. Hence, E∗D = (T ∗ U,0) is globally asymptotically stable in the region ℜ whenever R0 ≤ 1. 5 Endemic Equilibrium In this section, we first determine the endemic equilibrium points for the model (CM), establish its existence and then analyse its stability. 192 Ruchi Arora, Dharmendra Kumar, Ishita Jhamb and Avina Kaur CUBO 22, 2 (2020) 5.1 Endemic Equilibrium Points Let endemic equilibrium points be denoted by Ee = (S ∗∗,E∗∗,I∗∗,R′∗∗,R∗∗,X∗∗,Y ∗∗,Z∗∗). The components of Ee are obtained by imposing constant solutions in the model (CM) and solving the algebraic equations. By computations, we have S∗∗ = µNH ζNH + Z∗∗βBH , E ∗∗ = βBHZ ∗∗µ (α + ζ)(βBHZ∗∗ + ζNH) , I∗∗ = αβBHZ ∗∗µ (γ + ζ + ζ1)(α + ζ)(βBHZ∗∗ + ζNH) , R′∗∗ = αγβBHZ ∗∗µ (λ + ζ + ζ2)(γ + ζ + ζ1)(α + ζ)(βBHZ∗∗ + ζNH) , R∗∗ = λαβBHZ ∗∗µγ ζ(λ + ζ + ζ2)(γ + ζ + ζ1)(α + ζ)(βBHZ∗∗ + ζNH) , X∗∗ = ρ λM + κ , Y ∗∗ = ρλM (λM + κ)(ψ + κ) , Z∗∗ = ρψλM κ(λM + κ)(ψ + κ) . 5.2 Existence and Uniqueness of Endemic Equilibrium(E e ) Theorem 5.1 : Chikungunya Model (CM) has a unique endemic equilibrium if R0 > 1. As seen in section 2, λM = νBM (I ∗∗ + R ′ ∗∗) NH = νBMζαβµBHZ ∗∗(ζ + ζ2 + λ + γ) µ(βBHZ∗∗ + µ)(α + ζ)(ζ + ζ1 + γ)(ζ + ζ2 + λ) = R20µZ ∗∗κ2(ψ + κ) ρψ(βBHZ∗∗ + µ) Also, λH = βBHZ ∗∗ NH = βBHρψλM κNH(λM + κ)(ψ + κ) , or equivalently λM = λHµκ 2(ψ + κ) βBHρψζ − µκ(ψ + κ)λH . CUBO 22, 2 (2020) Mathematical Modeling of Chikungunya Dynamics . . . 193 Equating both values of λM , we get the following linear equation in terms of λH: λH(ρψβBH + R 2 0µκ(ψ + κ)) = (R 2 0 − 1)βBHρψζ. The unique solution to this equation exists and is given by λH = (R20 − 1)βBHρψζ ρψβBH + R 2 0µκ(ψ + κ) , which is positive if R20 > 1. This implies Z ∗∗ > 0, for R0 > 1. Hence, unique endemic equilibrium exists for R0 > 1. 5.3 Local Stability Theorem 5.2: The endemic equilibrium of the chikungunya model (CM) is locally asymptotically stable if R0 > 1. Proof: We evaluate the Jacobian matrix for the system of nonlinear differential equations corre- sponding to the model (CM). Let Je denote the Jacobian of the system at Ee (which exists for R0 > 1). Clearly, JE = (J1,J2,J3,J4,J5,J6,J7,J8) T where J1 = ( −βBHZ ∗∗ NH + βBHZ ∗∗ S ∗∗ (NH ) 2 − ζ, βBH Z ∗∗ S ∗∗ (NH ) 2 , βBH Z ∗∗ S ∗∗ (NH ) 2 , βBH Z ∗∗ S ∗∗ (NH ) 2 , βBH Z ∗∗ S ∗∗ (NH ) 2 ,0, 0, −βBH S ∗∗ NH ) , J2 = ( βBHZ ∗∗ NH − βBHZ ∗∗S∗∗ (NH ) 2 , −βBH Z ∗∗S∗∗ (NH ) 2 − α − ζ, −βBH Z ∗∗S∗∗ (NH ) 2 , −βBH Z ∗∗S∗∗ (NH ) 2 , −βBH Z ∗∗S∗∗ (NH ) 2 ,0, 0, βBH S ∗∗ NH ) , J3 = (0,α,−γ − ζ − ζ1,0, 0,0, 0,0) , J4 = (0,0,γ,−λ − ζ − ζ2,0, 0,0, 0) , J5 = (0,0, 0,λ,−ζ,0,0, 0) , J6 = ( νBM (I ∗∗+R ′ ∗∗)X∗∗ (NH ) 2 , νBM (I ∗∗+R ′ ∗∗)X∗∗ (NH ) 2 , νBM (I ∗∗+R ′ ∗∗)X∗∗ (NH ) 2 − νBM X ∗∗ NH , νBM (I ∗∗+R ′ ∗∗)X∗∗ (NH ) 2 − νBM X ∗∗ NH , 0, − νBM (I ∗∗+R ′ ∗∗) NH − κ,0, 0 ) , J7 = ( −νBM (I ∗∗+R ′ ∗∗)X∗∗ (NH ) 2 , −νBM (I ∗∗+R ′ ∗∗)X∗∗ (NH ) 2 , −νBM (I ∗∗+R ′ ∗∗)X∗∗ (NH ) 2 + νBM X ∗∗ NH , −νBM (I ∗∗+R ′ ∗∗)X∗∗ (NH ) 2 + νBM X ∗∗ NH , 0, νBM (I ∗∗+R ′ ∗∗) NH ,−κ − ψ,0 ) , J8 = (0,0, 0,0, 0,0,ψ,−κ) Further, we reduce JE to the following upper triangular matrix (UE). UE = (U1,U2,U3,U4,U5,U6,U7,U8) T where U1 = ( −βBHZ ∗∗ NH + βBHZ ∗∗ S ∗∗ (NH)2 − ζ, βBHZ ∗∗ S ∗∗ (NH)2 , βBHZ ∗∗ S ∗∗ (NH)2 , βBHZ ∗∗ S ∗∗ (NH)2 , βBH Z ∗∗ S ∗∗ (NH)2 ,0,0, −βBHS ∗∗ NH ) , U2 = ( 0, −βBHZ ∗∗ S ∗∗ (NH)2 − α − ζ, −βBHZ ∗∗ S ∗∗ (NH)2 , −βBHZ ∗∗ S ∗∗ (NH)2 , −βBHZ ∗∗ S ∗∗ (NH)2 ,0,0, βBHS ∗∗ NH ) , U3 = (0,0,−γ − ζ − ζ1,0,0,0,0,0), U4 = (0,0,0,−λ − ζ − ζ2,0,0,0,0) , U5 = (0,0,0,0,−ζ,0,0,0), U6 = ( 0,0,0,0,0,− νBM(I ∗∗+R ′ ∗∗) NH − κ,0,0 ) , U7 = (0,0,0,0,0,0,−κ − ψ,0) , U8 = (0,0,0,0,0,0,0,−κ) 194 Ruchi Arora, Dharmendra Kumar, Ishita Jhamb and Avina Kaur CUBO 22, 2 (2020) Attached are the eigen values of UE: ( −ζ,−κ,−ψ − κ,−γ − ζ − ζ1,−λ − ζ − ζ2,− νBM (I ∗∗+R ′ ∗∗) NH − κ, Z ∗∗ βBH(S ∗∗ −NH) (NH)2 − ζ, −βBHZ ∗∗ S ∗∗ (NH)2 − α − ζ ) each of which is negative and by the criterion given in [36], the endemic equilibrium point (Ee) is locally asymptotically stable if R0 > 1. 6 Numerical Simulation The values of parameters that would be used for simulation of the model (CM) are listed in Table 3. The values used for simulation are taken with reference to their ranges, as stated in Table 2. Fig. 2a and Fig. 2b are visualizations of the existence and stability of equilibria for the cases, R0 < 1 and R0 > 1, respectively. Also, it illustrates that for R0 < 1, the infection dies out over a period of time as it is the case of DFE. However, in the same time period, it can been seen that the infection continues to persist in the population when R0 > 1 as it is the case of EE. (a) (b) Figure 2: Total number of Infected Humans (I) with respect to time. In Fig. 3a, it is clear that the recuperated population ultimately falls down to zero for the case when R0 < 1, where finally the disease dies out and ultimately the entire population will shift to the recovered section with no more inflow into the recuperated part. In contrast, for the same time period, if R0 > 1 (Fig. 3b), the disease persists in the population. Therefore, we can see a substantial proportion of population which is still in the recuperated phase. Fig. 4 and Fig. 5, both show the time duration around which the number of infected popula- tion comes to a fall which is actually the same for recuperated population to reach the peak. In Fig. 6a, again for R0 < 1, as the disease dies out so it is evidently a situation when the population of the infectious mosquitoes dies out. In contrast to it, for R0 > 1 (Fig. 6b), the number CUBO 22, 2 (2020) Mathematical Modeling of Chikungunya Dynamics . . . 195 (a) (b) Figure 3: Total number of Recuperated Humans (R′) with respect to time. Figure 4: Total number of Infected (I) and Recuperated (R′) Humans when R0 < 1. of infectious mosquitoes continue to persist in population as it is the case of endemic equilibrium (EE). Fig. 7 shows the change in the number of infected, recuperated and recovered population with respect to time in accordance with model (CM) whereas Fig. 8 is a simulation of the model (CM) without recuperated section of population. The curve representing the recovered population in Fig. 8, is an increasing curve showing a rapid increase in the number of people attaining full recovery. But this does not fit in accordance to the case of Chikungunya infection. However, in Fig. 7, we can see the convexity of the curve repre- senting recovered population for a substantial period of time and this is because of the presence of recuperation factor which has been considered in our model. During this period, the recuperation curve is rising higher which is practically more relevant and well in consensus with the nature of this particular disease. 196 Ruchi Arora, Dharmendra Kumar, Ishita Jhamb and Avina Kaur CUBO 22, 2 (2020) Figure 5: Total number of Infected (I) and Recuperated (R′) Humans when R0 > 1. (a) (b) Figure 6: Total number of Infectious Mosquitoes (Z) with respect to time. 7 Conclusion In this paper, a new deterministic model is formulated to study the transmission dynamics of Chikungunya virus (CHIKV). Making a considerable refinement to the existing models present in the literature, a so far neglected section of human population is introduced, namely the population in the recuperation phase. The study shows that the disease free equilibrium (DFE) of the model is locally as well as globally asymptotically stable whenever existence of an associated reproduction number R0, is less than 1 and unstable otherwise. Also, an endemic equilibrium (EE) exists whenever R0 is greater than 1 and is locally asymptotically stable too. Simulations of the model make it evident that introduction of the said compartment is well justified, as this model provides a more realistic illustration for Chikungunya infection wherein the quantitative behaviour of disease has given a better visualisation. Moreover, the qualitative behaviour of the disease as studied by CUBO 22, 2 (2020) Mathematical Modeling of Chikungunya Dynamics . . . 197 Figure 7: Variation of Infected, Recuperated and Recovered Human Population with time for model (CM). Figure 8: Variation of Infected and Recovered Human Population with time for model (CM) without recuperation section. various researchers in [14] is very well taken into consideration through our model. If we do not consider the recuperation section in model (CM), then the following model becomes a special case of our model. It is clearly seen that our model (CM) gives a better illustration to the dynamics of the Chikungunya virus and hence, the proposed model is indeed 198 Ruchi Arora, Dharmendra Kumar, Ishita Jhamb and Avina Kaur CUBO 22, 2 (2020) more realistic and practical. dS dt = µ − βBHZS NH − ζS, dE dt = βBHZS NH − αE − ζE, dI dt = αE − γI − (ζ + ζ1)I, dR dt = γI − ζR, dX dt = ρ − νBMIX NH − κX, dY dt = νBMIX NH − ψY − κY, dZ dt = ψY − κZ, where NH(t) = S(t) + E(t) + I(t) + R(t). Comparison of the above model with our model (CM) is done in section 6 with the help of the graphs shown in Fig. 7 and Fig. 8. Acknowledgement The authors are thankful to Dr. Sukhanta Dutta, Dr. Harsha Kharbanda and Tanvi for their valuable suggestions and inputs. This research was supported and funded by the Science Centre, SGTB Khalsa College (Project Code: SGTBKC/SC/SP/2017/08/518). CUBO 22, 2 (2020) Mathematical Modeling of Chikungunya Dynamics . . . 199 References [1] WHO. Chikungunya. http://www.who.int/denguecontrol/arbo-viral/other aborvial chikun gunya/en/, 2014. [2] WHO. Chikungunya. http://www.who.int/mediacentre/factsheets/fs327/en/, 2014. [3] National Center for Biotechnology information. Chikungunya outbreak. http://www.ncbl.nlm.nih.gov/pmc/atricles/PMC4111161, 2014. [4] G. Pialoux, B. Gäuzère, S. Jauréguiberry, and M. Strobel, Chikungunya, an epidemic arbovi- rosis. The Lancet infectious diseases, 7:319–27, 06 2007. [5] K. Sergon, C. Njuguna, R. Kalani, V. Ofula, C. Onyango, L. Konongoi, S. Bedno, H. Burke, A. M. Dumilla, J. Konde, M. Kariuki Njenga, R. Sang, and R. Breiman. Seroprevalence of chikungunya virus (chikv) infection on lamu island, kenya, october 2004. The American journal of tropical medicine and hygiene, 78:333–7, 03 2008. [6] M. Barro, A. Guiro and D. Ouedraogo, Optimal control of a SIR epidemic model with general incidence function and a time delays, Cubo 20 (2018), no. 2, 53–66. [7] O.K. Oare, Impact and optimal control of movement on a multipatch hepatitis C virus model, TWMS J. Pure Appl. Math. 5 (2014), no. 1, 80–95. [8] Y. Dumont, F. Chiroleu and C. Domerg, On a temporal model for the Chikungunya disease: modeling, theory and numerics, Math. Biosci. 213 (2008), no. 1, 80–91. [9] D. Moulay, M. A. Aziz-Alaoui and M. Cadivel, The chikungunya disease: modeling, vector and transmission global dynamics, Math. Biosci. 229 (2011), no. 1, 50–63. [10] L. Yakob and A. C. A. Clements. A mathematical model of chikungunya dynamics and control: The major epidemic on runion island. PloS one, 8:e57448, 03 2013. [11] S. Naowarat and I. M. Tang. Transmission model of chikungunya fever in the presence of two species of aedes mosquitoes. American Journal of Applied Sciences, 10:449–459, 05 2013. [12] D. Hincapie-Palacio and J. Ospina. Mathematical modeling of chikungunya fever control. page 94870Z, 05 2015. [13] F. B. Agusto et al., Mathematical model of three age-structured transmission dynamics of Chikungunya virus, Comput. Math. Methods Med. 2016, Art. ID 4320514, 31 pp. [14] A. Mohan, D. H. N. Kiran, I. Chiranjeevi Manohar, and D. Prabath Kumar. Epidemiology, clinical manifestations, and diagnosis of chikungunya fever: Lessons learned from the re- emerging epidemic. Indian journal of dermatology, 55:54–63, 01 2010. 200 Ruchi Arora, Dharmendra Kumar, Ishita Jhamb and Avina Kaur CUBO 22, 2 (2020) [15] C. Manore, K. Hickmann, S. Xu, H. Wearing, and J. Hyman. Comparing dengue and chikun- gunya emergence and endemic transmission in a. aegypti and a. albopictus. Journal of Theo- retical Biology, 356:174191, 09 2014. [16] N. Chitnis, J. M. Hyman and J. M. Cushing, Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model, Bull. Math. Biol. 70 (2008), no. 5, 1272–1296. [17] H. Delatte, G. Gimonneau, A. Triboire, and D. Fontenille. Influence of temperature on imma- ture 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, 02 2009. [18] Y. Dumont and F. Chiroleu, Vector control for the Chikungunya disease, Math. Biosci. Eng. 7 (2010), no. 2, 313–345. [19] C. Lahariya and S. Pradhan. Emergence of chikungunya virus in indian subcontinent after 32 years: A review. Journal of vector borne diseases, 43:151–60, 01 2007. [20] M. R. Sebastian, R. Lodha, and S. Kabra. Chikungunya infection in children. The Indian Journal of Pediatrics, 76:185–89, 03 2009. [21] O. Schwartz and M. Albert. Biology and pathogenesis of chikungunya virus. nat rev microbiol 8:491-500. Nature reviews. Microbiology, 8:491–500, 07 2010. [22] K. Costanzo, K. Mormann, and S. Juliano. Asymmetrical competition and patterns of abun- dance of aedes albopictus and culex pipiens (diptera: Culicidae). Journal of medical entomol- ogy, 42:559–70, 08 2005. [23] N. A. Hashim, A. Hassan, O. N. Abu Tahir, M. Salmah, and N. Basari. Population analysis of aedes albopictus (skuse) (diptera:culicidae) under uncontrolled laboratory conditions. Tropical biomedicine, 25:117–25, 09 2008. [24] M. Dubrulle, L. Mousson, S. Moutailler, M. Vazeille, and A. B. Failloux. Chikungunya virus and aedes mosquitoes: Saliva is infectious as soon as two days after oral infection. PloS one, 4:e5895, 02 2009. [25] D. Mavalankar, P. Shastri, T. Bandyopadhyay, J. Parmar, and K. Ramani. Increased mortality rate associated with chikungunya epidemic, ahmedabad, india. Emerging infectious diseases, 14:412–5, 04 2008. [26] P. Poletti, G. Messeri, M. Ajelli, R. Vallorani, C. Rizzo, and S. Merler. Transmission potential of chikungunya virus and control measures: The case of italy. PloS one, 6:e18860, 05 2011. CUBO 22, 2 (2020) Mathematical Modeling of Chikungunya Dynamics . . . 201 [27] M. Turell, J. R. Beaman, and R. F. Tammariello. Susceptibility of selected strains of aedes aegypti and aedes albopictus (diptera: Culicidae) to chikungunya virus. Journal of medical entomology, 29:49–53, 02 1992. [28] E. Massad, S. Ma, M. Burattini, Y. Tun, F. Coutinho, and L. Ang. The risk of chikungunya fever in a dengueendemic area. Journal of travel medicine, 15:147–55, 05 2008. [29] K. Pesko, C. J. Westbrook, C. Mores, L. Philip Lounibos, and M. Reiskind. Effects of infec- tious virus dose and bloodmeal delivery method on susceptibility of aedes aegypti and aedes albopictus to chikungunya virus. Journal of medical entomology, 46:395–9, 04 2009. [30] V. Lakshmikantham, S. Leela, and A.A. Martynyuk. Stability Analysis of Nonlinear Systems. Systems & Control: Foundations & Applications. Springer International Publishing, 2015. [31] P. van den Driessche and J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci. 180 (2002), 29–48. [32] J. H. Jones. Notes on r0. Department of Anthropological Sciences Stanford University, 2007. [33] G.F. Simmons. Differential Equations with Applications and Historical Notes. Textbooks in Mathematics. CRC Press, 2016. [34] C. Bhunu, W. Garira, and Z. Mukandavire. Modeling hiv/aids and tuberculosis coinfection. Bulletin of mathematical biology, 71:1745–80, 06 2009. [35] C. Castillo-Chavez, Z. Feng and W. Huang, “On the Computation of R0 and Its Role on Global Stability,” In: Mathematical Approaches for Emerging and Reemerging Infectious Diseases: An Introduction, Springer-Verlag, New York, 2002, pp. 229-250. [36] EX. DeJesus and C. Kaufman. Routh-hurwitz criterion in the examination of eigenvalues of a system of nonlinear ordinary differential equations. Physical review. A, 35:5288–5290, 07 1987. Introduction Model Formulation Preliminary Results Positivity of Solutions Invariant Region Disease Free Equilibrium (DFE) Existence of Equilibrium Reproduction Number Local Stability Global Stability Endemic Equilibrium Endemic Equilibrium Points Existence and Uniqueness of Endemic Equilibrium(Ee) Local Stability Numerical Simulation Conclusion