# Dynamics of a Fractional-Order Chikungunya Model with Asymptomatic Infectious Class

Authors: Mlyashimbi Helikumi; Gideon Eustace; Steady Mushayabasa
Journal: Computational and Mathematical Methods in Medicine
 (2022)

---

## Abstract

In this paper, a nonlinear fractional-order chikungunya disease model that incorporates asymptomatic infectious individuals is proposed and analyzed. The main interest of this work is to investigate the role of memory effects on the dynamics of chikungunya. Qualitative analysis of the model’s equilibria showed that there exists a threshold quantity which governs persistence and extinction of the disease. Model parameters were estimated based on the 2015 weekly reported cases in Colombia. The Adams-Bashforth-Moulton method was used to numerically solve the proposed model. We investigated the role of asymptomatic infectious patients on short- and long-term dynamics of the diseases. We also determined threshold levels for the efficacy of preventative strategies that results in effective management of the disease. We believe that our model can provide invaluable insights for public health authorities to predict the effect of chikungunya transmission and analyze its underlying factors and to guide new control efforts.


---

## Body

## 1. Introduction

Emerging and reemerging diseases of vector-borne infections such as Zika virus and chikungunya pose a considerable public health problem worldwide [1]. Previous studies reported that vector-borne diseases are emerging at a growing rate and bearing a disproportionate segment of all new infectious diseases, the vast majority of them being viruses [1, 2]. These diseases are currently associated with high mortality and morbidity, which triggers immeasurable loss in many societies. For this reason, there is growing interest among researchers to utilize different techniques to study the dynamics of vector-borne infections with a goal to provide useful means for policymakers to evaluate potential ways of effectively managing these diseases [3, 4].Although various techniques can be used to understand the short- and long-term dynamics of vector-borne infections, mathematical modeling is now regarded as a standard and indispensable tool for understanding the mechanisms of interaction between host and vectors [5, 6]. Models are used to approach questions which are too complex, inaccessible, numerous, diverse, mutable, unique, dangerous, expensive, big, small, slow, or fast to approach by other means [7]. In this study, we seek to develop and analyze a mathematical framework for understanding the dynamics of chikungunya virus.Chikungunya is a viral infection whose causative agent is an alphavirus that infects humans through bites of Aedes mosquitoes [8]. The initial chikungunya cases were in Tanzania in the mid-1900s [9]. The clinical indications of chikungunya infections resemble those of dengue fever and include rashes and high fever [9, 10]. The onset of the disease in humans is often rapid and takes 5-7 days, and several case fatalities have been reported [11, 12]. Over the past 50 years, numerous cases of chikungunya reemergence have occurred in different parts of Africa and Asia [13, 14]. The infection has presently been established in nearly 40 countries but is endemic in Africa where its transmission is maintained between mosquitoes, primates, and humans [15]. The most notable outbreak occurred in India in 2005-2006 when the WHO approximated 1.3 million cases [16, 17].In the last two decades, a plethora of mathematical models have been proposed to explain and predict as well as quantify the effectiveness ways of managing chikungunya virus (see, for example, [8, 18–26] and references therein). Yakob and Clements [8] constructed and analyzed a deterministic mathematical model for chikungunya virus that incorporated two infected human subpopulations designated as symptomatic and asymptomatic. Among several outcomes, their study demonstrated the strong influence that both the latent period of infection in humans and the prepatent period have on the dynamics of the disease. In [20], a spatial stochastic model was utilized to show that perifocal vector control is capable of limiting the spread of chikungunya. In [23], a stochastic model that incorporated climate-based mosquito population was proposed to identify temporal windows that have epidemic risk in the United States (US). Among several outcomes, their work strongly suggested that, in the event of an introduction and establishment of chikungunya in the US, endemic and epidemic regions would emerge initially, primarily defined by environmental factors controlling annual mosquito population cycles. In [24], a system of ordinary differential equations was employed to investigate and understand the importance of different model parameters on the dynamics of chikungunya. In [26], a mathematical framework was proposed to model and analyze virus mutation dynamics of chikungunya outbreaks. Results from their study suggested that the virus mutation dynamics could play an important role in the transmission of chikungunya. Hence, there is a need to better understand the mutation mechanism.Despite these efforts, however, not much has been done to investigate the role of memory effects on chikungunya virus dynamics. Previous studies suggest that memory effects are inherent in many biological phenomena and models based on integer-order differentiation do not adequately capture the memory effects [27, 28]. Based on this assertion, fractional-order calculus has been widely and extensively used by several researchers recently, to investigate the role of memory effects in different diseases [29–33]. The fractional derivative has several definitions such as those derived from Riemann-Liouville, Caputo, Liouville, Weyl, Riesz, Grünwald-Litnikov, Marchaud and Hifler, Caputo-Fabrizio-Caputo, Atanga-Baleanu, Atanga-beta derivative, M-fractional derivative, conformable derivative, Atangana-Koca, Atanga-Gomez, and variable-order and fractal-fractional idea [34, 35].Motivated by the above-mentioned works, we derive a fractional-order model for chikungunya virus based on the Caputo derivative. The choice of using the Caputo derivative is also aided by the fact that the Caputo derivative for a given function which is constant is zero. Thus, the Caputo operator computes an ordinary differential equation, followed by a fractional integral to obtain the desired order of fractional derivative [36]. Most importantly, the Caputo fractional derivative allows the use of local initial conditions to be included in the derivation of the model [33, 36].The remaining part of this paper is organized as follows: Section2 contains some few basic definitions on fractional calculus. These definitions will be utilized to establish several important results. In Section 3, we propose and qualitatively analyze a fractional-order chikungunya virus model. In particular, we computed the reproduction number and demonstrated that it is an important threshold quantity for disease persistence and extinction. Section 4 contains numerical illustrations and discussions. The paper is concluded in Section 5


## 2. Preliminaries on Fractional Calculus

Throughout the paper, we use the Caputo fractional-order derivatives because the initial conditions of fractional differential equations with Caputo derivative take on the same form as that of the integer-order ones, which have received more attention in modeling and analysis of many real-world phenomena despite their inability to capture memory effects, which are inherent in these phenomena. In this section, we will present some essential definitions and lemmas that will be utilized to determine the dynamical behavior of the proposed model.Definition 1 (see [37]).
The Riemann-Liouville fractional integral operator of orderθof a continuous functionf:t0,+∞⟶ℝis defined as(1)Iθft=1Γθ∫0tt−sθ−1fsds,where Γ· is the gamma function and is defined by Γz=∫0∞e−ttz−1dt and θ>0.Definition 2 (see [37]).
The Caputo fractional derivative of orderθfor a functionf∈Cnt0,+∞,ℝis defined by(2)t0cDtθft=1Γn−θ∫t0tfnst−sθ−n+1ds,where t≥t0 and n is a positive integer such that n−1<θ<n∈ℕ. In particular, if 0<θ<1, then the Caputo fractional derivative takes the form
(3)bcDtθft=1Γ1−θ∫t0tf′st−sθds.Lemma 3 (see [38]).
The solution to the Cauchy problem(4)Dtθt0cxt=λxt+ft,xa=bb∈ℝwith 0<θ<1 and λ∈ℝ has the form
(5)xt=bEθλt−aθ+∫att−sθ−1Eθ,θλt−sθfsds,while the solution to the problem
(6)Dtθt0cxt=λxt,xa=bb∈ℝis given by
(7)xt=bEθλt−aθ,where Eθ· is the Mittag-Leffler function and defined by
(8)Eθz=∑k=0∞zkΓθk+1,θ>0,z∈ℂ.Lemma 4 (see [39]).
Letutbe a continuous function ont0,+∞and satisfying(9)Dtθt0cut=−λut+μ,ut0=ut0,with 0<θ<1, λ,μ∈ℝ2, and λ≠0, and t0≥0 is the initial time. Then,
(10)ut≤ut0−μλEθ−λt−t0θ+μλ.Lemma 5 (see [40]).
Letxtbe a continuous and differentiable function withxt∈ℝ+. Then, for any time instantt≥t0, one has(11)t0cDtθxt−x∗−x∗lnxtx∗≤1−x∗xtt0cDtθxt,x∗∈ℝ+,∀θ∈0,1.


## 3. Model Formulation and Analytical Results

### 3.1. Model Formulation

Owing to the merits of fractional derivatives (discussed in Introduction) in modeling real-world problems in comparison to integer ordinary differential equations, in this section, we propose a fractional-order model for chikungunya infection. The proposed model (12) is based on the Caputo fractional operator which is known to be the most suitable for modeling biological and physical phenomena [36]. The proposed model is governed by the following assumptions:
(i)
The total populations of humansNht are subdivided into classes of susceptible Sht, exposed Eht, asymptomatic infectious individuals Aht, symptomatic infectious individuals Iht, and humans who have recovered from infection Rht. The total vector population Nvt is subdivided into compartments of susceptible Svt, exposed Evt, and infectious Ivt. Once infected, vectors are assumed to remain infectious for their entire life span. Throughout this paper, we will use the subscripts h and v to represent the human and vector population, respectively(ii)
Disease transmission is assumed to occur solely when there is interaction between the host and vector. Thus, all new recruits in the both the host and vector are assumed to be susceptible to infection. In most cases, chikungunya disease does not lead to death, but its symptoms can be severe and disabling [41]. Based on this assertion, we assumed a constant size population in both the host and vector with a recruitment and non-disease-related mortality rate modeled by μj with j=h,v. Further, new recruits are assumed susceptible and their recruitment rate is proportional to the population and is given by μjNj. Meanwhile, we assume that whenever there is effective contact between a susceptible individual and an infectious vector, disease transmission will occur at rate βv. Similarly, let βh be the transmission rate of the disease from a infectious human to a susceptible vector whenever there is effective contact. Asymptomatic infectious individuals are assumed to have less parasite load compared to symptomatic infectious individuals. To account for this aspect, in our model formulation, we introduce a factor 1−ηh to model the reduction of infectivity of the asymptomatic individuals(iii)
Although there is no vaccine to prevent chikungunya virus infection, humans can minimize chances of contracting the infections through a number of strategies such as the use of insect repellent, wearing of long-sleeved shirts and pants, and treating clothing and gear. Let1−εh model the effects of preventative strategies utilized by humans to minimize chances of contracting the infection (naturally or through treatment)(iv)
Exposed humans are assumed to incubate the infection for1/αh days after which a proportion fh become asymptomatic infectious patients and the remainder, 1−fh, become symptomatic infectious patients. The average infectious period of asymptomatic infectious individuals is modeled by 1/γh, and symptomatic infectious individuals are assumed to be infectious for an average period of 1/σh days(v)
Exposed vectors incubate the disease for1/αv days after which they become infectious. Once infected, vectors are assumed to remain infectious for their entire life spanBased on the above assumptions, the proposed model (12) is summarized in Figure 1 and akes the following form:
(12)bcDtθSht=μhNht−1−εhβvcIvtNhtSht−μhSht,bcDtθEht=1−εhβvcIvtNhtSht−αh+μhEht,bcDtθAht=fhαhEht−μh+γhAht,bcDtθIht=1−fhαhEht−μh+σhIht,bcDθRht=γhAht+σhIht−μhRht,bcDtθSvt=μvNvt−1−εhβhcIht+1−ηhAhtNhtSvt−μvSvt,bcDtθEvt=1−εhβhcIht+1−ηhAhtNhtSvt−αv+μvEvt,bcDtθIvt=αvEvt−μvIvt,where θ∈0,1 is the order of the fractional derivative. Model (12) is subject to the initial conditions:
(13)Sj0=Sj0,Ej0=Ej0,Ah0=Ah0,Ij0=Ij0,Rh0=Rh0,for j=h,v. By adding all the equations that govern the dynamics of the disease in the host, one can observed that the human population is assumed to be constant. Similarly, by adding all the equations for the vector population, one can also observe that the vector population is assumed to be constant. Based on this, we can consider normalized populations for both the vector and host. Let
(14)sjt=SjtNjt,ejt=EjtNjt,ijt=IjtNjt,aht=AhtNht,rht=RhtNht,ϕ=NvtNht.Using the definitions in (18), the chikungunya model with normalized populations is given as follows:
(15)bcDtθsht=μh−1−εhβvϕcivsh−μhsh,bcDtθeht=1−εhβvϕcivsh−m1eh,bcDtθaht=fhαheh−m2ah,bcDtθiht=1−fhαheh−m3ih,bcDθrht=γhah+σhih−μhrh,bcDtθsvt=μv−1−εhβhcih+1−ηhahsv−μvsv,bcDtθevt=1−εβhcih+1−ηhahsv−m4ev,bcDtθivt=αvev−μviv,with
(16)m1=αh+μh,m2=μh+γh,m3=μh+σh,m4=αv+μv,subject to the initial conditions
(17)sj0=sj0,ej0=ej0,ah0=aj0,ij0=ij0,rh0=rh0,for j=h,v. Since the variable rht does not appear in any other equations of system (12), it suffices to analyze the dynamics of chikungunya virus infection from the following reduced system:
(18)bcDtθsht=μh−1−εhβvϕcivsh−μhsh,bcDtθeht=1−εhβvϕcivsh−m1eh,bcDtθaht=fhαheh−m2ah,bcDtθiht=1−fhαheh−m3ih,bcDtθsvt=μv−1−εβhcih+1−ηhahsv−μvsv,bcDtθev=1−εβhcih+1−ηhahsv−m4ev,bcDtθivt=αvev−μviv.


### 3.2. Positivity and Boundedness of Model Solutions

Theorem 6.
The fractional order (18) has a unique solution, which remains inℝ7,and the closed set(19)Ω=sh,eh,ah,ihsv,ev,iv∈ℝ+7sh≥0,eh≥0,ah≥0,ih≥0sh+eh+ah+ih≤1sv≥0,ev≥0,iv≥0sv+ev+iv≤1is a positive invariant set of system (18).Proof.
We begin by demonstrating that the solution of system (18) is always nonnegative for all t≥0. Based on system (18), we have (for j=h,v(20)t0cDtθshtsj=0=μj≥0,t0cDtθehteh=0=1−εhβvϕcivsh≥0,t0cDtθahtah=0=fhαheh≥0,t0cDtθihtih=0=1−fhαheh≥0,t0cDtθtev=0=1−εhβhcih+1−ηhahsv≥0,t0cDtθivtiv=0=αvev≥0.
Results in equation (20) demonstrate that the vector field given by the right-hand side of (18) on each coordinate plane either is tangent to the coordinate plane or points to the interior of ℝ+7. Hence, the domain ℝ+7 is a positively invariant region. Moreover, if the initial conditions of system (18) are nonnegative, then it follows that the corresponding solutions of model (18) are nonnegative. Let nht=sht+eht+aht+iht. By adding all the equations for the host population in (6), one gets
(21)t0cDtθnht≤μh−μhnh.
It follows from Lemma3 that (38) has a solution of the form
(22)nht≤nh0−1Eθ−μhtθ+1.
SinceEθ−μhtθ≥0, so that when nh0≤1, we have nht≤1. Similarly, by adding all the equations for the vector population (with, nv=sv+ev+iv), one gets
(23)t0cDtθnvt=μv−μvnv.
Again, by Lemma3, we have
(24)nvt≤nv0−1Eθ−μvtθ+1.
By similar arguments as before, we havenvt≤1. This completes the proof.☐


### 3.3. Model Equilibria, Their Existence, and Global Stability

Direct calculations of system (18) show that the proposed model admits two equilibrium points, namely, the disease-free equilibrium (DFE) point (denoted by E0) and the endemic equilibrium (EE) point (denoted by E∗) point, which are, respectively, given by E0:sh0,eh0,ah0,ih0,sv0,ev0,iv0=1,0,0,0,1,0,0. The endemic equilibrium (EE) point E∗=Sh∗,Eh∗,Ih∗,Sv∗,Ev∗,Iv∗ with
(25)sh∗=m4βh1−εhcαhμh1−fhm2+1−ηhfhm3H1,eh∗=μvμhm2m3m4H1R02−1,ah∗=μvμhm3m4fhαhH1R02−1,ih∗=μvμhm2m41−fhαhH1R02−1,sv∗=μvm1m2m3βv1−εhϕcαv+μhm4H2,ev∗=μv2μhm1m2m3βv1−εϕcαvH2R02−1iv∗=μvμhm1m2m3βv1−εhϕcH2R02−1,with
(26)H1=βhc1−εhαhβvϕ1−εhcαv+μhm41−fhm2+1−ηfhm3,H2=1−fhm2αhβh1−εhcμh+m31−ηfhαh1−εhβhμh+m1m2μv,R0=βvβhϕc21−εh2αvαhμvm1m41−ηhfhm2+1−fhm3.In equation (25), we can note that the disease persists in the community if R0>1. Precisely, the threshold quantity R0 demonstrates the power of the disease to persist or become extinct in the host population. Thus, the expression R0 defines the basic reproduction number of the proposed fractional-order model.Theorem 7.
IfR0<1, then the DFE of system ((18)) is globally asymptotically stable inΩ, otherwise it is unstable.Proof.
By considering only the infected compartments from (18), we can write
(27)bcDtθx=F−Vx,where x=eh,ah,ih,ev,ivT, with F and V defined as follows:
(28)F=0000βv1−εθc00000000000βh1−ε1−ηcβh1−εc0000000,V=m10000−fhαhm2000−1−fhαh0m300000m40000−αvμv.
One can verify by direct calculations thatV−1F is a nonnegative and irreducible matrix and ρV−1F=R0. It follows from the Perron-Frobenius theorem [42] that V−1F has positive left eigenvector w associated with R0, i.e.,
(29)wV−1F=R0w.
SincewV−1 is a positive vector, we propose the following Lyapunov function to study the global stability of DFE:
(30)Lt=wV−1x.
DifferentiatingL along solutions of (12) leads to
(31)bcDtθLt=wV−1bcDtθx≤wV−1F−Vx=R0−1wx≤0ifR0≤1.
It can be easily verified that the largest invariant subset ofΓ where bcDtθLt=0 is the singleton E0. Therefore, by LaSalle’s invariance principle [43], E0 is globally asymptotically stable in Γ when R0≤1.☐Theorem 8.
IfR0>1,then the endemic equilibriumE∗of system ((18)) is globally asymptotically stable inΩ.Proof.
To prove Theorem8, we consider the following Lyapunov functional:
(32)Ut=sh−sh∗−sh∗lnshsh∗+eh−eh∗−eh∗lneheh∗+w1ah−ah∗−ah∗lnahah∗w2ih−ih∗−ih∗lnihih∗+w3sv−sv∗−sv∗lnsvsv∗+w3ev−ev∗−ev∗lnevev∗+w4iv−iv∗−iv∗lniviv∗.
Letgih,ah=ih+ηhah. At the endemic equilibrium, we have the following identities:
(33)μh=1−εβvϕciv∗sh∗+μhsh∗m1eh∗=1−εβvϕciv∗sh∗,m2ah∗=fhαheh∗,m3ih∗=1−fhαheh∗,μv=1−εβhcgih∗,ah∗sv∗−μvsv∗,m4ev∗=1−εβhgih∗,ah∗sv∗,μviv∗=αvev∗.
Let(34)w1=βviv∗sh∗1−fhαheh∗,w2=βviv∗sh∗fhαheh∗,w3=βviv∗sh∗βhgih∗,ah∗sv∗,w4=βviv∗sh∗αvev∗.
After some algebraic manipulations, one gets(35)t0cDUt≤μh2−shsh∗−sh∗sh+μvw32−SvSv∗−Sv∗Sv+βviv∗sh∗7−sh∗sh−sheh∗ivsh∗ehiv∗−ahah∗−ehah∗eh∗ah+eheh∗−ehih∗eh∗ih−ihih∗−sv∗sv+gih,ahgih∗,ah∗−svevgih∗,ah∗sv∗ev∗gih∗,ah∗−eviv∗ev∗iv.
Since the arithmetic mean is greater than or equal to the geometric mean, it follows that terms in the brackets are(36)2≤sjsi0+si0sij=h,v,7≤sh∗sh+sheh∗ivsh∗ehiv∗+ahah∗+ehah∗eh∗ah−eheh∗+ehih∗eh∗ih+ihih∗+sv∗sv−gih,ahgih∗,ah∗+svevgih∗,ah∗sv∗ev∗gih∗,ah∗+eviv∗ev∗iv.
Hence, we conclude thatt0cDUt≤0 for all sh,eh,ah,ih,sv,ev,iv≥0 since sh∗,eh∗,ah∗,ih∗sv∗,ev∗, and iv∗ are nonnegative whenever R0>1. Therefore, by Lasalle’s invariance principle [43], the endemic equilibrium point is globally asymptotically stable whenever R0>1.☐


## 3.1. Model Formulation

Owing to the merits of fractional derivatives (discussed in Introduction) in modeling real-world problems in comparison to integer ordinary differential equations, in this section, we propose a fractional-order model for chikungunya infection. The proposed model (12) is based on the Caputo fractional operator which is known to be the most suitable for modeling biological and physical phenomena [36]. The proposed model is governed by the following assumptions:
(i)
The total populations of humansNht are subdivided into classes of susceptible Sht, exposed Eht, asymptomatic infectious individuals Aht, symptomatic infectious individuals Iht, and humans who have recovered from infection Rht. The total vector population Nvt is subdivided into compartments of susceptible Svt, exposed Evt, and infectious Ivt. Once infected, vectors are assumed to remain infectious for their entire life span. Throughout this paper, we will use the subscripts h and v to represent the human and vector population, respectively(ii)
Disease transmission is assumed to occur solely when there is interaction between the host and vector. Thus, all new recruits in the both the host and vector are assumed to be susceptible to infection. In most cases, chikungunya disease does not lead to death, but its symptoms can be severe and disabling [41]. Based on this assertion, we assumed a constant size population in both the host and vector with a recruitment and non-disease-related mortality rate modeled by μj with j=h,v. Further, new recruits are assumed susceptible and their recruitment rate is proportional to the population and is given by μjNj. Meanwhile, we assume that whenever there is effective contact between a susceptible individual and an infectious vector, disease transmission will occur at rate βv. Similarly, let βh be the transmission rate of the disease from a infectious human to a susceptible vector whenever there is effective contact. Asymptomatic infectious individuals are assumed to have less parasite load compared to symptomatic infectious individuals. To account for this aspect, in our model formulation, we introduce a factor 1−ηh to model the reduction of infectivity of the asymptomatic individuals(iii)
Although there is no vaccine to prevent chikungunya virus infection, humans can minimize chances of contracting the infections through a number of strategies such as the use of insect repellent, wearing of long-sleeved shirts and pants, and treating clothing and gear. Let1−εh model the effects of preventative strategies utilized by humans to minimize chances of contracting the infection (naturally or through treatment)(iv)
Exposed humans are assumed to incubate the infection for1/αh days after which a proportion fh become asymptomatic infectious patients and the remainder, 1−fh, become symptomatic infectious patients. The average infectious period of asymptomatic infectious individuals is modeled by 1/γh, and symptomatic infectious individuals are assumed to be infectious for an average period of 1/σh days(v)
Exposed vectors incubate the disease for1/αv days after which they become infectious. Once infected, vectors are assumed to remain infectious for their entire life spanBased on the above assumptions, the proposed model (12) is summarized in Figure 1 and akes the following form:
(12)bcDtθSht=μhNht−1−εhβvcIvtNhtSht−μhSht,bcDtθEht=1−εhβvcIvtNhtSht−αh+μhEht,bcDtθAht=fhαhEht−μh+γhAht,bcDtθIht=1−fhαhEht−μh+σhIht,bcDθRht=γhAht+σhIht−μhRht,bcDtθSvt=μvNvt−1−εhβhcIht+1−ηhAhtNhtSvt−μvSvt,bcDtθEvt=1−εhβhcIht+1−ηhAhtNhtSvt−αv+μvEvt,bcDtθIvt=αvEvt−μvIvt,where θ∈0,1 is the order of the fractional derivative. Model (12) is subject to the initial conditions:
(13)Sj0=Sj0,Ej0=Ej0,Ah0=Ah0,Ij0=Ij0,Rh0=Rh0,for j=h,v. By adding all the equations that govern the dynamics of the disease in the host, one can observed that the human population is assumed to be constant. Similarly, by adding all the equations for the vector population, one can also observe that the vector population is assumed to be constant. Based on this, we can consider normalized populations for both the vector and host. Let
(14)sjt=SjtNjt,ejt=EjtNjt,ijt=IjtNjt,aht=AhtNht,rht=RhtNht,ϕ=NvtNht.Using the definitions in (18), the chikungunya model with normalized populations is given as follows:
(15)bcDtθsht=μh−1−εhβvϕcivsh−μhsh,bcDtθeht=1−εhβvϕcivsh−m1eh,bcDtθaht=fhαheh−m2ah,bcDtθiht=1−fhαheh−m3ih,bcDθrht=γhah+σhih−μhrh,bcDtθsvt=μv−1−εhβhcih+1−ηhahsv−μvsv,bcDtθevt=1−εβhcih+1−ηhahsv−m4ev,bcDtθivt=αvev−μviv,with
(16)m1=αh+μh,m2=μh+γh,m3=μh+σh,m4=αv+μv,subject to the initial conditions
(17)sj0=sj0,ej0=ej0,ah0=aj0,ij0=ij0,rh0=rh0,for j=h,v. Since the variable rht does not appear in any other equations of system (12), it suffices to analyze the dynamics of chikungunya virus infection from the following reduced system:
(18)bcDtθsht=μh−1−εhβvϕcivsh−μhsh,bcDtθeht=1−εhβvϕcivsh−m1eh,bcDtθaht=fhαheh−m2ah,bcDtθiht=1−fhαheh−m3ih,bcDtθsvt=μv−1−εβhcih+1−ηhahsv−μvsv,bcDtθev=1−εβhcih+1−ηhahsv−m4ev,bcDtθivt=αvev−μviv.


## 3.2. Positivity and Boundedness of Model Solutions

Theorem 6.
The fractional order (18) has a unique solution, which remains inℝ7,and the closed set(19)Ω=sh,eh,ah,ihsv,ev,iv∈ℝ+7sh≥0,eh≥0,ah≥0,ih≥0sh+eh+ah+ih≤1sv≥0,ev≥0,iv≥0sv+ev+iv≤1is a positive invariant set of system (18).Proof.
We begin by demonstrating that the solution of system (18) is always nonnegative for all t≥0. Based on system (18), we have (for j=h,v(20)t0cDtθshtsj=0=μj≥0,t0cDtθehteh=0=1−εhβvϕcivsh≥0,t0cDtθahtah=0=fhαheh≥0,t0cDtθihtih=0=1−fhαheh≥0,t0cDtθtev=0=1−εhβhcih+1−ηhahsv≥0,t0cDtθivtiv=0=αvev≥0.
Results in equation (20) demonstrate that the vector field given by the right-hand side of (18) on each coordinate plane either is tangent to the coordinate plane or points to the interior of ℝ+7. Hence, the domain ℝ+7 is a positively invariant region. Moreover, if the initial conditions of system (18) are nonnegative, then it follows that the corresponding solutions of model (18) are nonnegative. Let nht=sht+eht+aht+iht. By adding all the equations for the host population in (6), one gets
(21)t0cDtθnht≤μh−μhnh.
It follows from Lemma3 that (38) has a solution of the form
(22)nht≤nh0−1Eθ−μhtθ+1.
SinceEθ−μhtθ≥0, so that when nh0≤1, we have nht≤1. Similarly, by adding all the equations for the vector population (with, nv=sv+ev+iv), one gets
(23)t0cDtθnvt=μv−μvnv.
Again, by Lemma3, we have
(24)nvt≤nv0−1Eθ−μvtθ+1.
By similar arguments as before, we havenvt≤1. This completes the proof.☐


## 3.3. Model Equilibria, Their Existence, and Global Stability

Direct calculations of system (18) show that the proposed model admits two equilibrium points, namely, the disease-free equilibrium (DFE) point (denoted by E0) and the endemic equilibrium (EE) point (denoted by E∗) point, which are, respectively, given by E0:sh0,eh0,ah0,ih0,sv0,ev0,iv0=1,0,0,0,1,0,0. The endemic equilibrium (EE) point E∗=Sh∗,Eh∗,Ih∗,Sv∗,Ev∗,Iv∗ with
(25)sh∗=m4βh1−εhcαhμh1−fhm2+1−ηhfhm3H1,eh∗=μvμhm2m3m4H1R02−1,ah∗=μvμhm3m4fhαhH1R02−1,ih∗=μvμhm2m41−fhαhH1R02−1,sv∗=μvm1m2m3βv1−εhϕcαv+μhm4H2,ev∗=μv2μhm1m2m3βv1−εϕcαvH2R02−1iv∗=μvμhm1m2m3βv1−εhϕcH2R02−1,with
(26)H1=βhc1−εhαhβvϕ1−εhcαv+μhm41−fhm2+1−ηfhm3,H2=1−fhm2αhβh1−εhcμh+m31−ηfhαh1−εhβhμh+m1m2μv,R0=βvβhϕc21−εh2αvαhμvm1m41−ηhfhm2+1−fhm3.In equation (25), we can note that the disease persists in the community if R0>1. Precisely, the threshold quantity R0 demonstrates the power of the disease to persist or become extinct in the host population. Thus, the expression R0 defines the basic reproduction number of the proposed fractional-order model.Theorem 7.
IfR0<1, then the DFE of system ((18)) is globally asymptotically stable inΩ, otherwise it is unstable.Proof.
By considering only the infected compartments from (18), we can write
(27)bcDtθx=F−Vx,where x=eh,ah,ih,ev,ivT, with F and V defined as follows:
(28)F=0000βv1−εθc00000000000βh1−ε1−ηcβh1−εc0000000,V=m10000−fhαhm2000−1−fhαh0m300000m40000−αvμv.
One can verify by direct calculations thatV−1F is a nonnegative and irreducible matrix and ρV−1F=R0. It follows from the Perron-Frobenius theorem [42] that V−1F has positive left eigenvector w associated with R0, i.e.,
(29)wV−1F=R0w.
SincewV−1 is a positive vector, we propose the following Lyapunov function to study the global stability of DFE:
(30)Lt=wV−1x.
DifferentiatingL along solutions of (12) leads to
(31)bcDtθLt=wV−1bcDtθx≤wV−1F−Vx=R0−1wx≤0ifR0≤1.
It can be easily verified that the largest invariant subset ofΓ where bcDtθLt=0 is the singleton E0. Therefore, by LaSalle’s invariance principle [43], E0 is globally asymptotically stable in Γ when R0≤1.☐Theorem 8.
IfR0>1,then the endemic equilibriumE∗of system ((18)) is globally asymptotically stable inΩ.Proof.
To prove Theorem8, we consider the following Lyapunov functional:
(32)Ut=sh−sh∗−sh∗lnshsh∗+eh−eh∗−eh∗lneheh∗+w1ah−ah∗−ah∗lnahah∗w2ih−ih∗−ih∗lnihih∗+w3sv−sv∗−sv∗lnsvsv∗+w3ev−ev∗−ev∗lnevev∗+w4iv−iv∗−iv∗lniviv∗.
Letgih,ah=ih+ηhah. At the endemic equilibrium, we have the following identities:
(33)μh=1−εβvϕciv∗sh∗+μhsh∗m1eh∗=1−εβvϕciv∗sh∗,m2ah∗=fhαheh∗,m3ih∗=1−fhαheh∗,μv=1−εβhcgih∗,ah∗sv∗−μvsv∗,m4ev∗=1−εβhgih∗,ah∗sv∗,μviv∗=αvev∗.
Let(34)w1=βviv∗sh∗1−fhαheh∗,w2=βviv∗sh∗fhαheh∗,w3=βviv∗sh∗βhgih∗,ah∗sv∗,w4=βviv∗sh∗αvev∗.
After some algebraic manipulations, one gets(35)t0cDUt≤μh2−shsh∗−sh∗sh+μvw32−SvSv∗−Sv∗Sv+βviv∗sh∗7−sh∗sh−sheh∗ivsh∗ehiv∗−ahah∗−ehah∗eh∗ah+eheh∗−ehih∗eh∗ih−ihih∗−sv∗sv+gih,ahgih∗,ah∗−svevgih∗,ah∗sv∗ev∗gih∗,ah∗−eviv∗ev∗iv.
Since the arithmetic mean is greater than or equal to the geometric mean, it follows that terms in the brackets are(36)2≤sjsi0+si0sij=h,v,7≤sh∗sh+sheh∗ivsh∗ehiv∗+ahah∗+ehah∗eh∗ah−eheh∗+ehih∗eh∗ih+ihih∗+sv∗sv−gih,ahgih∗,ah∗+svevgih∗,ah∗sv∗ev∗gih∗,ah∗+eviv∗ev∗iv.
Hence, we conclude thatt0cDUt≤0 for all sh,eh,ah,ih,sv,ev,iv≥0 since sh∗,eh∗,ah∗,ih∗sv∗,ev∗, and iv∗ are nonnegative whenever R0>1. Therefore, by Lasalle’s invariance principle [43], the endemic equilibrium point is globally asymptotically stable whenever R0>1.☐


## 4. Numerical Results and Discussions

In this section, we will make use of MATLAB programming language to simulate model (18) so as to support analytical findings and determine the implications of time-dependent controls. To simulate model (12), we make use of the technique so-called generalized Adams-Bashforth-Moulton (ABM) method [36]. For any nonlinear fractional differential equation,
(37)bcDtθψt=ft,ψt,0≤t≤T,with the following initial conditions:
(38)ψm0=ψ0m,m=0,1,2,3,⋯,θ−1.Now, with operating by the fractional integral operator on equation (37), we can obtain on the solution ψt by solving the following equation:
(39)ψt=∑m=0θ−1ψ0mm!tm+1Γθ∫0tt−τθ−1fτ,ψτdτ.Diethelm [44] used the predictor-corrector scheme based on the ABM algorithm to numerically solve (39). Setting h=T/N, tn=nh, and n=0,1,2,⋯,N∈ℤ+, equation (39) can be discretized as follows:
(40)ψhtn+1=∑m=0∣θ∣−1ψ0mm!tn+1m+hθΓθ+2∑q=0naq,n+1ftq,ψq+hθΓθ+2ftn+1,ψn+1p,where
(41)aq,n+1=nθ+1−n−θn+θθ,q=0,n−q+2θ+1+n−qθ+1−2n−q+1θ+1,1≤q≤n,1,ifq=n+1,and the predicted valueψhptn+1 is determined by
(42)ψtn+1p=∑m=0∣θ∣−1ψ0mm!tn+1m+1Γθ∑q=0nbq,n+1ftq,ψhtq,with
(43)bq,n+1=hθθn+1−qθ−n−qθ.The error estimate is(44)max0≤q≤kψtq−ψhtq=Ohp,with k∈ℕ and p=min2,1+θ.


### 4.1. Application of the ABM Method to the Proposed Model

In this subsection, we solve numerically the nonlinear fractional model using the ABM method. In view of the generalized ABM method, the numerical scheme for the proposed model (18) is given in the following form [36]:
(45)shtn+1=sh0+hθΓθ+2fshtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1+hθΓθ+2∑q=0naq,n+1fshtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,ehtn+1=eh0+hθΓθ+2fehtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1+hθΓθ+2∑q=0naq,n+1fehtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,ahtn+1=ah0+hθΓθ+2fahtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1+hθΓθ+2∑q=0naq,n+1fahtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,ihtn+1=ih0+hθΓθ+2fihtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1+hθΓθ+2∑q=0naq,n+1fihtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,svtn+1=sv0+hθΓθ+2fsvtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1+hθΓθ+2∑q=0naq,n+1fsvtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,evtn+1=ev0+hθΓθ+2fevtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1+hθΓθ+2∑q=0naq,n+1fevtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,ivtn+1=iv0+hθΓθ+2fivtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1+hθΓθ+2∑q=0naq,n+1fsvtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,where
(46)shptn+1=sh0+1Γθ∑q=0nbq,n+1fshtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,ehptn+1=eh0+1Γθ∑q=0nbq,n+1fehtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,ahptn+1=ah0+1Γθ∑q=0nbq,n+1fahtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,ihptn+1=ih0+1Γθ∑q=0nbq,n+1fihtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,svptn+1=sv0+1Γθ∑q=0nbq,n+1fsvtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,evptn+1=ev0+1Γθ∑q=0nbq,n+1fshtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,shptn+1=iv0+1Γθ∑q=0nbq,n+1fivtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq.Further, we have(47)fshtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq=μh−1−εhβvqcivsh−μhsh,(48)fehtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq=1−εhβvqcivsh−m1eh,(49)fahtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq=fhαheh−m2ah,(50)fihtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq=1−fhαheh−m3ih,(51)fsvtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq=μv−1−εβhcihsv−μvsv−1−εhβhc1−ηhahsv,(52)fevtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq=1−εβhcihsv−m4ev+1−εβhc1−ηhahsv,(53)fivtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq=αvev−μviv.In addition, the quantities(54)fshtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1,fehtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1,fahtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1,fihtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1,fsvtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1,fevtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1,fivtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1are derived from system (53), at the points tn+1, n=1,2,3,⋯.,m.


### 4.2. Fractional-Order and Parameter Estimation Using Real Chikungunya Data

In this section, we will utilize the chikungunya data for Colombia (weekly cases observed in 2015) presented in [24] to numerically solve system (18) and determine the best fractional-order parameter θ that minimizes the deviations between the real data and model estimates. Data constitute of weekly reported chikungunya cases for Colombia in the year 2015. Despite it being a challenging task, estimation of fractional-order and model parameters is an integral part in mathematical modeling of infectious diseases. Several techniques such as the maximum likelihood iterated filtering and the nonlinear least-squares curve fitting are often used to validate proposed epidemiological models with data.The fitting process in this paper was done by making use of the least-squares method and Nelder-Mead algorithm [45]. The Nelder-Mead algorithm was used to perform a broad search of the parameter space and is less dependent on initial guesses. Once a good fit was determined, these parameters as the initial guess to search for a more localized region. We fitted the model to cumulative daily new infection data presented in [24]. Since our proposed model constitutes a normalized population, we have normalized the weekly cases by assuming that the total human population in the area where these cases were reported is assumed to be around 19,471,223 [24]. The cumulative new infections predicted by model (18), Ct, are given by the solution of the following equation:
(55)bcDtθCt=αhθ1−fheht.In order to compute the best fitting, we implemented the function(56)F:ℝθ,ηh,ϕ,εh,βh5⟶ℝθ,ηh,ϕ,εh,βh,where θ,ηh,ϕ,εh,βh are variables such that
(1)
for a givenθ,ηh,ϕ,εh,βh, solve numerically the system of differential equation (18) and equation (55) to obtain a solution Y^it=s^hi,e^hi,a^hi,i^hi,s^vi,e^vi,i^vi, which is an approximation of the real-world data Yt(2)
Sett0=1 (the fitting process starts in week 1) and for t=2,3,⋯,52, corresponding to weeks in where data are available, evaluate the computed numerical solution for iht; that is., i^h2, i^h3,…..,i^h52(3)
Compute the root mean square (RMSE) of the difference betweeni^h2,i^h3,.⋯,i^h52 and observed cases. This function F returns the root-mean-square error (RMSE) where(57)RMSE=1n∑k=152ihk−i∧hk2(4)
Determine a global minimum for the RMSE using Nelder-Mead algorithmThe functionF takes values in ℝ5 and returns a positive real number, the RMSE that measures the closeness of the model predictions to the observed data. On performing the fitting process, we assumed the following initial conditions eh0=0.0002, ah=ih0=0.0001, sh0=1−eh0−ah0−ih0, ev0=0.0002, iv0=0.0001, and sv0=1−ev0−iv0 and the baseline values for the model parameters are in Table 1.Figure 1
Flow chart for chikungunya disease transmission.Table 1
Description of parameters used in the model system (18).
SymbolDescriptionValueUnitsμhNatural mortality rate of humans0.00004 day−1[24]1/αhIncubation period of in humans12 (5–12) days−1[20, 24]cBiting rate of human by female mosquito0.19 (0.19–0.39) days−1[20, 24]fhProportion of asymptomatic patients0.1 (0.1–0.3)[20]1/γhInfectious period for asymptomatic patients10 (7–10) days−1[20]1/σhInfectious period for symptomatic patients7 (7–10) days−1[20, 24]ϕThe number of mosquitoes per individuals3FittingηhReduction of infectivity from asymptomatic patients to susceptible mosquito0.2 (0–1.0)FittingεhEfficacy of preventative strategies0.3 [0–1]FittingβhProbability of infection to be transmitted from humans to mosquito per bite0.9 (0.8–1)FittingβvProbability of infection to be transmitted from an infectious to a susceptible human per bite0.244 (0.06–0.244)[24, 46]1/μvLife expectancy of mosquitoes30 (4–30) day−1[19, 20]1/αvIncubation of mosquitoes3 (3–7) day−1[24, 46]R0Basic reproduction number1.5–7[20]Figure2 depicts the root-mean-square error (RMSE) for different derivative orders. The simulation results show that a global minimum error of estimation for the given data set occurs for θ=0.64 with RMSE=2.08×10−7. Figure 3 illustrates the model estimates for θ=0.64 and θ=1 (the classical model). From the illustration, one can be able to compare the predictability strength between the integer and fractional models. From the illustrations, we can observe that estimates for the integer model are close to the real data for the first 4 weeks; thereafter, the estimates significantly deviate from the reported cases compared to the estimates of the fractional-order model. Hence, we conclude that the fractional model presents better forecasts compared to the integer model.Figure 2
The root-mean-square error (RMSE) of the model estimation for different derivative orders. The minimum error estimation is obtained forθ=0.640, and it corresponds to RMSE=2.08×10−7.Figure 3
The estimation of the fractional-order model withθ=0.64 with RMSE=2.08×10−7 and R0=1.53.Figure4 shows the simulation of residuals against the predicted values of chikungunya cases in Colombia. It was noted that the residuals did not follow any particular path (exhibited a random pattern), suggesting that the model system (18) was a good fit well to chikungunya cases in Colombia as reported in [24].Figure 4
Simulation of residuals against predicted values of chikungunya cases in Colombia for the model system (18) as reported in [24].


### 4.3. Sensitivity of the Reproduction Number to Model Parameters

From the results in Theorems7 and 8, it is apparent that basic reproduction number R0 is an integral parameter for chikungunya persistence and extinction in the community; hence, there is a need to determine the sensitivity of R0 to each parameter. One of the techniques of determining the sensitivity of R0 to each parameter is through the computation of an elasticity index; that is, for a parameter ν, the elasticity index is found by the formulaν/R0∂R0/∂ν. Model parameters that have a positive elasticity index will increase the magnitude of R0 whenever they are increased while those with negative index will decrease R0 whenever they are increased. Thus, this linearized sensitivity analysis gives an idea of parameters that are vital in reducing R0 below unity in order for the disease to die out in the community. Results on computations on the elasticity indices are shown in Figure 5.Figure 5
Sensitivity analysis of the basic reproduction numberR0 with respect to model parameters. Baseline values used are in Table 1.From the results in Figure5, one can observe that the mosquito feeding (biting) rate (parameter c) has the largest elasticity index than all the other model parameters that define R0. In particular, if the mosquitoes increase their feeding rate by 10%, then the magnitude of R0 will also increase by 10%. Parameters βh (transmission rate of the disease from humans to mosquito per bite), βv (transmission rate of the disease from mosquito to human per bite), ϕ (the number of mosquitoes per individual), and μv (mortality rate of mosquitoes) have a similar elasticity index in absolute value. However, an increase in βh, βv, and ϕ will increase the size of R0 while an increase in μv reduces the size of R0. It is also worth noting that the efficacy of preventative strategies εh has a significant effect on reducing the size of R0. An increase in the efficacy of preventative strategies by 10% will decrease the size of R0 by 4.28%. From the results, one can conclude that the recovery of symptomatic patients, mortality rate of mosquitoes, and efficacy of preventative strategies are capable of reducing the size of R0 by a significant margin whenever they are increased.Simulation results in Figure6(a) depict the effects of varying mosquito feeding (biting) rate, c on R0. The results show that an increase in mosquito biting rate (modeled by parameter c) increases the size of R0. In particular, whenever c>0.13, then R0>1 which implies that the disease persists in the community. Figure 6(b) demonstrates the effects of ε (efficacy of preventative strategies) on R0. From the results, we can conclude that whenever the efficacy of preventative strategies is less than 54% all the time, then the disease persists in the community. However, if the reverse is true, it dies off. Figure 6(c) shows a contour plot of R0 as a function of c and ε. We observe that when c becomes larger or when ε is reduced, R0 increases, implying a higher disease risk.Figure 6
Effects of varying (a) mosquito feeding (biting) rate,c on R0; (b) the efficacy of preventative strategies, ε on R0; (c) both ε and c on R0.
(a)(b)(c)Numerical illustrations in Figure7(a) show the effects of varying mosquito mortality rate, μv on R0. From the results, we note that if the life span of mosquitoes is less than 14 days (μv=0.07 day−1), then the disease dies out, otherwise it persists. In Figure 7(b), we explore the impact of the number of mosquitoes per individuals, ϕ on R0. We can note that whenever the number of vectors per human in the area is greater than 1.3, then the infection persists in the area, otherwise it dies out.Figure 7
Effects of varying (a) mosquito mortality rate,μv on R0, and (b) the number of mosquitoes per individuals, ϕ on R0.
(a)(b)Numerical simulation in Figure8(a) shows the effects of varying average infectious period for symptomatic patients, 1/σh on R0. From the results, we note that if the average infectious period for symptomatic patients is greater than 40%, the disease dies out, otherwise it persists. In Figure 7(b), we assessed the effect of probability of infection to be transmitted from an infectious mosquito to a susceptible human per bite, βv on R0. We can note that whenever the probability of transmission from infectious mosquito to susceptible human is greater than 10%, the disease persists in the area, otherwise it dies out.Figure 8
Effects of varying (a) average infectious period for symptomatic patients,1/σh on R0, and (b) the probability of infection to be transmitted from an infectious to a susceptible human per bite, βv on R0.
(a)(b)


### 4.4. Simulation Results to Support Analytical Findings in the Study

In this section, we will carry out numerical simulations so as to support analytical findings in the study. Simulation results in Figure9 show the convergence of model solutions to the disease-free equilibrium with different derivative orders. We can note that all model solutions converge to the DFE point despite different derivative orders and this is in agreement with Theorem 7. In Figure 10, we can observe that order of the derivative has a significant effect on dynamical behavior on the infected vector and host population over time. In particular, we observe that when the derivative order θ is reduced from 1, the memory effect of the system increases, and as a result, the infection grows slowly and the number of infected vectors and host increases over time.Figure 9
Dynamical solutions of model (18) with different order derivatives. The solutions were obtained upon setting ε=0.6 giving R0=0.8722. Overall, we can note that as the order of the derivative approaches unity, then time taken by model solutions to converge to the unique DFE equilibrium increases.
(a)(b)(c)(d)(e)Figure 10
Dynamical solutions of model (18) with different order derivatives. The solutions were obtained upon setting ε=0.1 giving R0=1.9625. Overall, we observe that when the derivative order θ is reduced from 1, the memory effect of the system increases, and as a result, the infection grows slowly and the number of infected vectors and host increases over time.
(a)(b)(c)(d)(e)To gain insights on the role of asymptomatic patients on chikungunya dynamics, we simulated model (18) with different values of fh (proportion of individuals who become asymptomatic patients upon the completion of the incubation period) and the results are shown in Figure 11. The results show that as fh increases, the exposed human population, asymptomatic patients, and infectious vectors increase remarkably while the infectious human population decreases slightly. Although the exposed vector population increases with increasing fh, the increase is not highly significant compared to the other populations.Figure 11
Numerical results for different values offh (proportion of individuals who become asymptomatic patients upon the completion of the incubation period) over time with θ=0.64.
(a)(b)(c)(d)(e)To explore the impact of the preventative strategies on long-term dynamics, we simulated model (18) with different values of εh with a fixed derivative order (θ=0.64) and the results are in Figure 12. Overall, the results show that as εh approaches unity, the disease dies out in the community. In particular, the closer εh is to unity, the shorter the time it takes for the disease to become extinct.Figure 12
Numerical results showing the effects of varyingεh on long-term dynamics of the disease.
(a)(b)(c)(d)(e)


## 4.1. Application of the ABM Method to the Proposed Model

In this subsection, we solve numerically the nonlinear fractional model using the ABM method. In view of the generalized ABM method, the numerical scheme for the proposed model (18) is given in the following form [36]:
(45)shtn+1=sh0+hθΓθ+2fshtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1+hθΓθ+2∑q=0naq,n+1fshtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,ehtn+1=eh0+hθΓθ+2fehtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1+hθΓθ+2∑q=0naq,n+1fehtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,ahtn+1=ah0+hθΓθ+2fahtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1+hθΓθ+2∑q=0naq,n+1fahtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,ihtn+1=ih0+hθΓθ+2fihtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1+hθΓθ+2∑q=0naq,n+1fihtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,svtn+1=sv0+hθΓθ+2fsvtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1+hθΓθ+2∑q=0naq,n+1fsvtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,evtn+1=ev0+hθΓθ+2fevtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1+hθΓθ+2∑q=0naq,n+1fevtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,ivtn+1=iv0+hθΓθ+2fivtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1+hθΓθ+2∑q=0naq,n+1fsvtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,where
(46)shptn+1=sh0+1Γθ∑q=0nbq,n+1fshtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,ehptn+1=eh0+1Γθ∑q=0nbq,n+1fehtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,ahptn+1=ah0+1Γθ∑q=0nbq,n+1fahtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,ihptn+1=ih0+1Γθ∑q=0nbq,n+1fihtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,svptn+1=sv0+1Γθ∑q=0nbq,n+1fsvtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,evptn+1=ev0+1Γθ∑q=0nbq,n+1fshtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq,shptn+1=iv0+1Γθ∑q=0nbq,n+1fivtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq.Further, we have(47)fshtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq=μh−1−εhβvqcivsh−μhsh,(48)fehtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq=1−εhβvqcivsh−m1eh,(49)fahtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq=fhαheh−m2ah,(50)fihtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq=1−fhαheh−m3ih,(51)fsvtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq=μv−1−εβhcihsv−μvsv−1−εhβhc1−ηhahsv,(52)fevtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq=1−εβhcihsv−m4ev+1−εβhc1−ηhahsv,(53)fivtq,shtq,ehtq,ahtq,ihtq,svtq,evtq,ivtq=αvev−μviv.In addition, the quantities(54)fshtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1,fehtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1,fahtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1,fihtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1,fsvtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1,fevtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1,fivtn+1,shPtn+1,ehptn+1,ahptn+1,ihptn+1,svPtn+1,evptn+1,ivptn+1are derived from system (53), at the points tn+1, n=1,2,3,⋯.,m.


## 4.2. Fractional-Order and Parameter Estimation Using Real Chikungunya Data

In this section, we will utilize the chikungunya data for Colombia (weekly cases observed in 2015) presented in [24] to numerically solve system (18) and determine the best fractional-order parameter θ that minimizes the deviations between the real data and model estimates. Data constitute of weekly reported chikungunya cases for Colombia in the year 2015. Despite it being a challenging task, estimation of fractional-order and model parameters is an integral part in mathematical modeling of infectious diseases. Several techniques such as the maximum likelihood iterated filtering and the nonlinear least-squares curve fitting are often used to validate proposed epidemiological models with data.The fitting process in this paper was done by making use of the least-squares method and Nelder-Mead algorithm [45]. The Nelder-Mead algorithm was used to perform a broad search of the parameter space and is less dependent on initial guesses. Once a good fit was determined, these parameters as the initial guess to search for a more localized region. We fitted the model to cumulative daily new infection data presented in [24]. Since our proposed model constitutes a normalized population, we have normalized the weekly cases by assuming that the total human population in the area where these cases were reported is assumed to be around 19,471,223 [24]. The cumulative new infections predicted by model (18), Ct, are given by the solution of the following equation:
(55)bcDtθCt=αhθ1−fheht.In order to compute the best fitting, we implemented the function(56)F:ℝθ,ηh,ϕ,εh,βh5⟶ℝθ,ηh,ϕ,εh,βh,where θ,ηh,ϕ,εh,βh are variables such that
(1)
for a givenθ,ηh,ϕ,εh,βh, solve numerically the system of differential equation (18) and equation (55) to obtain a solution Y^it=s^hi,e^hi,a^hi,i^hi,s^vi,e^vi,i^vi, which is an approximation of the real-world data Yt(2)
Sett0=1 (the fitting process starts in week 1) and for t=2,3,⋯,52, corresponding to weeks in where data are available, evaluate the computed numerical solution for iht; that is., i^h2, i^h3,…..,i^h52(3)
Compute the root mean square (RMSE) of the difference betweeni^h2,i^h3,.⋯,i^h52 and observed cases. This function F returns the root-mean-square error (RMSE) where(57)RMSE=1n∑k=152ihk−i∧hk2(4)
Determine a global minimum for the RMSE using Nelder-Mead algorithmThe functionF takes values in ℝ5 and returns a positive real number, the RMSE that measures the closeness of the model predictions to the observed data. On performing the fitting process, we assumed the following initial conditions eh0=0.0002, ah=ih0=0.0001, sh0=1−eh0−ah0−ih0, ev0=0.0002, iv0=0.0001, and sv0=1−ev0−iv0 and the baseline values for the model parameters are in Table 1.Figure 1
Flow chart for chikungunya disease transmission.Table 1
Description of parameters used in the model system (18).
SymbolDescriptionValueUnitsμhNatural mortality rate of humans0.00004 day−1[24]1/αhIncubation period of in humans12 (5–12) days−1[20, 24]cBiting rate of human by female mosquito0.19 (0.19–0.39) days−1[20, 24]fhProportion of asymptomatic patients0.1 (0.1–0.3)[20]1/γhInfectious period for asymptomatic patients10 (7–10) days−1[20]1/σhInfectious period for symptomatic patients7 (7–10) days−1[20, 24]ϕThe number of mosquitoes per individuals3FittingηhReduction of infectivity from asymptomatic patients to susceptible mosquito0.2 (0–1.0)FittingεhEfficacy of preventative strategies0.3 [0–1]FittingβhProbability of infection to be transmitted from humans to mosquito per bite0.9 (0.8–1)FittingβvProbability of infection to be transmitted from an infectious to a susceptible human per bite0.244 (0.06–0.244)[24, 46]1/μvLife expectancy of mosquitoes30 (4–30) day−1[19, 20]1/αvIncubation of mosquitoes3 (3–7) day−1[24, 46]R0Basic reproduction number1.5–7[20]Figure2 depicts the root-mean-square error (RMSE) for different derivative orders. The simulation results show that a global minimum error of estimation for the given data set occurs for θ=0.64 with RMSE=2.08×10−7. Figure 3 illustrates the model estimates for θ=0.64 and θ=1 (the classical model). From the illustration, one can be able to compare the predictability strength between the integer and fractional models. From the illustrations, we can observe that estimates for the integer model are close to the real data for the first 4 weeks; thereafter, the estimates significantly deviate from the reported cases compared to the estimates of the fractional-order model. Hence, we conclude that the fractional model presents better forecasts compared to the integer model.Figure 2
The root-mean-square error (RMSE) of the model estimation for different derivative orders. The minimum error estimation is obtained forθ=0.640, and it corresponds to RMSE=2.08×10−7.Figure 3
The estimation of the fractional-order model withθ=0.64 with RMSE=2.08×10−7 and R0=1.53.Figure4 shows the simulation of residuals against the predicted values of chikungunya cases in Colombia. It was noted that the residuals did not follow any particular path (exhibited a random pattern), suggesting that the model system (18) was a good fit well to chikungunya cases in Colombia as reported in [24].Figure 4
Simulation of residuals against predicted values of chikungunya cases in Colombia for the model system (18) as reported in [24].


## 4.3. Sensitivity of the Reproduction Number to Model Parameters

From the results in Theorems7 and 8, it is apparent that basic reproduction number R0 is an integral parameter for chikungunya persistence and extinction in the community; hence, there is a need to determine the sensitivity of R0 to each parameter. One of the techniques of determining the sensitivity of R0 to each parameter is through the computation of an elasticity index; that is, for a parameter ν, the elasticity index is found by the formulaν/R0∂R0/∂ν. Model parameters that have a positive elasticity index will increase the magnitude of R0 whenever they are increased while those with negative index will decrease R0 whenever they are increased. Thus, this linearized sensitivity analysis gives an idea of parameters that are vital in reducing R0 below unity in order for the disease to die out in the community. Results on computations on the elasticity indices are shown in Figure 5.Figure 5
Sensitivity analysis of the basic reproduction numberR0 with respect to model parameters. Baseline values used are in Table 1.From the results in Figure5, one can observe that the mosquito feeding (biting) rate (parameter c) has the largest elasticity index than all the other model parameters that define R0. In particular, if the mosquitoes increase their feeding rate by 10%, then the magnitude of R0 will also increase by 10%. Parameters βh (transmission rate of the disease from humans to mosquito per bite), βv (transmission rate of the disease from mosquito to human per bite), ϕ (the number of mosquitoes per individual), and μv (mortality rate of mosquitoes) have a similar elasticity index in absolute value. However, an increase in βh, βv, and ϕ will increase the size of R0 while an increase in μv reduces the size of R0. It is also worth noting that the efficacy of preventative strategies εh has a significant effect on reducing the size of R0. An increase in the efficacy of preventative strategies by 10% will decrease the size of R0 by 4.28%. From the results, one can conclude that the recovery of symptomatic patients, mortality rate of mosquitoes, and efficacy of preventative strategies are capable of reducing the size of R0 by a significant margin whenever they are increased.Simulation results in Figure6(a) depict the effects of varying mosquito feeding (biting) rate, c on R0. The results show that an increase in mosquito biting rate (modeled by parameter c) increases the size of R0. In particular, whenever c>0.13, then R0>1 which implies that the disease persists in the community. Figure 6(b) demonstrates the effects of ε (efficacy of preventative strategies) on R0. From the results, we can conclude that whenever the efficacy of preventative strategies is less than 54% all the time, then the disease persists in the community. However, if the reverse is true, it dies off. Figure 6(c) shows a contour plot of R0 as a function of c and ε. We observe that when c becomes larger or when ε is reduced, R0 increases, implying a higher disease risk.Figure 6
Effects of varying (a) mosquito feeding (biting) rate,c on R0; (b) the efficacy of preventative strategies, ε on R0; (c) both ε and c on R0.
(a)(b)(c)Numerical illustrations in Figure7(a) show the effects of varying mosquito mortality rate, μv on R0. From the results, we note that if the life span of mosquitoes is less than 14 days (μv=0.07 day−1), then the disease dies out, otherwise it persists. In Figure 7(b), we explore the impact of the number of mosquitoes per individuals, ϕ on R0. We can note that whenever the number of vectors per human in the area is greater than 1.3, then the infection persists in the area, otherwise it dies out.Figure 7
Effects of varying (a) mosquito mortality rate,μv on R0, and (b) the number of mosquitoes per individuals, ϕ on R0.
(a)(b)Numerical simulation in Figure8(a) shows the effects of varying average infectious period for symptomatic patients, 1/σh on R0. From the results, we note that if the average infectious period for symptomatic patients is greater than 40%, the disease dies out, otherwise it persists. In Figure 7(b), we assessed the effect of probability of infection to be transmitted from an infectious mosquito to a susceptible human per bite, βv on R0. We can note that whenever the probability of transmission from infectious mosquito to susceptible human is greater than 10%, the disease persists in the area, otherwise it dies out.Figure 8
Effects of varying (a) average infectious period for symptomatic patients,1/σh on R0, and (b) the probability of infection to be transmitted from an infectious to a susceptible human per bite, βv on R0.
(a)(b)


## 4.4. Simulation Results to Support Analytical Findings in the Study

In this section, we will carry out numerical simulations so as to support analytical findings in the study. Simulation results in Figure9 show the convergence of model solutions to the disease-free equilibrium with different derivative orders. We can note that all model solutions converge to the DFE point despite different derivative orders and this is in agreement with Theorem 7. In Figure 10, we can observe that order of the derivative has a significant effect on dynamical behavior on the infected vector and host population over time. In particular, we observe that when the derivative order θ is reduced from 1, the memory effect of the system increases, and as a result, the infection grows slowly and the number of infected vectors and host increases over time.Figure 9
Dynamical solutions of model (18) with different order derivatives. The solutions were obtained upon setting ε=0.6 giving R0=0.8722. Overall, we can note that as the order of the derivative approaches unity, then time taken by model solutions to converge to the unique DFE equilibrium increases.
(a)(b)(c)(d)(e)Figure 10
Dynamical solutions of model (18) with different order derivatives. The solutions were obtained upon setting ε=0.1 giving R0=1.9625. Overall, we observe that when the derivative order θ is reduced from 1, the memory effect of the system increases, and as a result, the infection grows slowly and the number of infected vectors and host increases over time.
(a)(b)(c)(d)(e)To gain insights on the role of asymptomatic patients on chikungunya dynamics, we simulated model (18) with different values of fh (proportion of individuals who become asymptomatic patients upon the completion of the incubation period) and the results are shown in Figure 11. The results show that as fh increases, the exposed human population, asymptomatic patients, and infectious vectors increase remarkably while the infectious human population decreases slightly. Although the exposed vector population increases with increasing fh, the increase is not highly significant compared to the other populations.Figure 11
Numerical results for different values offh (proportion of individuals who become asymptomatic patients upon the completion of the incubation period) over time with θ=0.64.
(a)(b)(c)(d)(e)To explore the impact of the preventative strategies on long-term dynamics, we simulated model (18) with different values of εh with a fixed derivative order (θ=0.64) and the results are in Figure 12. Overall, the results show that as εh approaches unity, the disease dies out in the community. In particular, the closer εh is to unity, the shorter the time it takes for the disease to become extinct.Figure 12
Numerical results showing the effects of varyingεh on long-term dynamics of the disease.
(a)(b)(c)(d)(e)


## 5. Conclusions

In this paper, we have formulated a mathematical model chikungunya virus transmission that incorporates the human and mosquito populations. The proposed model is based on the Caputo fractional derivative. Fractional calculus has been widely used in mathematical modeling of evolutionary systems with memory effect on dynamics. Previous studies suggest that biological and physical phenomena are characterized by memory effects. Although there are several definitions of fractional derivatives, our choice of using Caputo derivative is as follows: (i) the Caputo derivative for a given function which is constant is zero. Thus, the Caputo operator computes an ordinary differential equation, followed by a fractional integral to obtain the desired order of fractional derivative; (ii) the Caputo fractional derivative allows the use of local initial conditions to be included in the derivation of the model. We computed the model equilibria and determined their existence and global stability. We have shown that the model has two equilibrium points: the disease-free equilibrium (DFE) and the endemic equilibrium (EE), which are both globally stable whenever the reproduction number is less than unity and greater than unity, respectively. We used the 2015 weekly reported cases in Colombia to estimate model parameters. We carried out a comprehensive sensitivity analysis of the model parameters, so as to determine the correlation between the model parameters and the reproduction number. Utilizing results on the sensitivity analysis, we determined threshold values for different model parameters that results in either extinction or persistence of the disease. In particular, this analysis was carried out for model parameters that were found to have a strong correlation with the reproduction number.The Adams-Bashforth-Moulton method was used to numerically solve the proposed model, and we observed that when the derivative order is reduced from 1, the memory effect of the system increases, and as a result, the infection grows slowly and the number of infected vectors and host increases over time. We believe that our model can provide invaluable insights for public health authorities to predict the effect of chikungunya transmission and analyze its underlying factors and to guide new control efforts.



---

*Source: 5118382-2022-02-07.xml*