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

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.


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][19][20][21][22][23][24][25][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,     [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: Section 2 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

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 realworld 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 function f : ½t 0 , +∞Þ ⟶ ℝ is defined as where Γð·Þ is the gamma function and is defined by ΓðzÞ = Ð ∞ 0 e −t t z−1 dt and θ > 0: Definition 2 (see [37]). The Caputo fractional derivative of order θ for a function f ∈ C n ð½t 0 ,+∞Þ, ℝÞ is defined by where t ≥ t 0 and n is a positive integer such that n − 1 < θ     Table 1. 3 Computational and Mathematical Methods in Medicine < n ∈ ℕ: In particular, if 0 < θ < 1, then the Caputo fractional derivative takes the form Lemma 3 (see [38]). The solution to the Cauchy problem with 0 < θ < 1 and λ ∈ ℝ has the form while the solution to the problem   Computational and Mathematical Methods in Medicine is given by where E θ ð·Þ is the Mittag-Leffler function and defined by Lemma 4 (see [39]). Let uðtÞ be a continuous function on ½ t 0 , +∞Þ and satisfying with 0 < θ < 1, ðλ, μÞ ∈ ℝ 2 , and λ ≠ 0, and t 0 ≥ 0 is the initial time. Then, Lemma 5 (see [40]). Let xðtÞ be a continuous and differentiable function with xðtÞ ∈ ℝ + . Then, for any time instant t ≥ t 0 , one has

Model Formulation.
Owing to the merits of fractional derivatives (discussed in Introduction) in modeling real-world problems in comparison to integer ordinary differential  (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:    (18) with different order derivatives. The solutions were obtained upon setting ε = 0:1 giving R 0 = 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. 7 Computational and Mathematical Methods in Medicine  Computational and Mathematical Methods in Medicine (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 μ j N j : 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. Let 1 − ε 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 for 1/α h days after which a proportion f h become asymptomatic infectious patients and the remainder, ð1 − f h Þ, 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 for 1/α v days after which they become infectious. Once infected, vectors are assumed to remain infectious for their entire life span Based on the above assumptions, the proposed model (12) is summarized in Figure 1 and akes the following form:      Computational and Mathematical Methods in Medicine where θ ∈ ð0, 1 is the order of the fractional derivative.
Model (12) is subject to the initial conditions: 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 Using the definitions in (18), the chikungunya model with normalized populations is given as follows: with subject to the initial conditions for j = h, v. Since the variable r h ðtÞ 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:

Positivity and Boundedness of Model Solutions
Theorem 6. The fractional order (18) has a unique solution, which remains in ℝ 7 , and the closed set is a positive invariant set of system (18).
Proof. We begin by demonstrating that the solution of system (18) 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) It follows from Lemma 3 that (38) has a solution of the form Since E θ ð−μ h t θ Þ ≥ 0, so that when n h ð0Þ ≤ 1, we have n h ðtÞ ≤ 1. Similarly, by adding all the equations for the vector population (with, Again, by Lemma 3, we have By similar arguments as before, we have n v ðtÞ ≤ 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 diseasefree equilibrium (DFE) point (denoted by E 0 ) and the endemic equilibrium (EE) point (denoted by E * ) point, which are, respectively, given by E 0 : with In equation (25), we can note that the disease persists in the community if R 0 > 1. Precisely, the threshold quantity R 0 demonstrates the power of the disease to persist or become extinct in the host population. Thus, the expression R 0 defines the basic reproduction number of the proposed fractional-order model.

Theorem 7.
If R 0 < 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 where x = ðe h , a h , i h , e v , i v Þ T , with F and V defined as follows: One can verify by direct calculations that V −1 F is a nonnegative and irreducible matrix and ρðV −1 FÞ = R 0 . It follows from the Perron-Frobenius theorem [42] that V −1 F has 12 Computational and Mathematical Methods in Medicine positive left eigenvector w associated with R 0 , i.e., Since wV −1 is a positive vector, we propose the following Lyapunov function to study the global stability of DFE: Differentiating L along solutions of (12) leads to It can be easily verified that the largest invariant subset of Γ where c b D θ t LðtÞ = 0 is the singleton fE 0 g. Therefore, by LaSalle's invariance principle [43], E 0 is globally asymptotically stable in Γ when R 0 ≤ 1.☐ Theorem 8. If R 0 > 1, then the endemic equilibrium E * of system ( (18)) is globally asymptotically stable in Ω: Proof. To prove Theorem 8, we consider the following Lyapunov functional: Let gði h , a h Þ = i h + η h a h . At the endemic equilibrium, we have the following identities: Let After some algebraic manipulations, one gets Since the arithmetic mean is greater than or equal to the geometric mean, it follows that terms in the brackets are Hence, we conclude that c and i * v are nonnegative whenever R 0 > 1. Therefore, by Lasalle's invariance principle [43], the endemic equilibrium point is globally asymptotically stable whenever R 0 > 1.☐

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, with the following initial conditions: Now, with operating by the fractional integral operator on equation (37), we can obtain on the solution ψðtÞ by solving the following equation: Diethelm [44] used the predictor-corrector scheme based on the ABM algorithm to numerically solve (39). Setting h = T/N, t n = nh, and n = 0, 1, 2, ⋯, N ∈ ℤ + , equation (39) can be discretized as follows: 13 Computational and Mathematical Methods in Medicine and the predicted value ψ p h ðt n+1 Þ is determined by The error estimate is with k ∈ ℕ and p = min ð2, 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]: are derived from system (53), at the points t n+1 , n = 1, 2, 3, ⋯:, m:

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), CðtÞ, are given by the solution of the following equation: In order to compute the best fitting, we implemented the function where θ, η h , ϕ, ε h , β h are variables such that (1) for a given ðθ, η h , ϕ, ε h , β h Þ, solve numerically the system of differential equation (18) Table 1. Figure 2 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 fractionalorder model. Hence, we conclude that the fractional model presents better forecasts compared to the integer model. Figure 4 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].

Sensitivity of the Reproduction Number to Model
Parameters. From the results in Theorems 7 and 8, it is apparent that basic reproduction number R 0 is an integral parameter for chikungunya persistence and extinction in the community; hence, there is a need to determine the sensitivity of R 0 to each parameter. One of the techniques of determining the sensitivity of R 0 to each parameter is through the computation of an elasticity index; that is, for a parameter ν, the elasticity index is found by the for-mulaðν/R 0 Þð∂R 0 /∂νÞ. Model parameters that have a positive elasticity index will increase the magnitude of R 0 whenever they are increased while those with negative index will decrease R 0 whenever they are increased. Thus, this linearized sensitivity analysis gives an idea of parameters that are vital in reducing R 0 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.
From the results in Figure 5, one can observe that the mosquito feeding (biting) rate (parameter c) has the largest elasticity index than all the other model parameters that define R 0 . In particular, if the mosquitoes increase their feeding rate by 10%, then the magnitude of R 0 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 R 0 while an increase in μ v reduces the size of R 0 . It is also worth noting that the efficacy of preventative strategies ε h has a significant effect on reducing the size of R 0 . An increase in the efficacy of preventative strategies by 10% will decrease the size of R 0 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 R 0 by a significant margin whenever they are increased.
Simulation results in Figure 6(a) depict the effects of varying mosquito feeding (biting) rate, c on R 0 . The results show that an increase in mosquito biting rate (modeled by parameter c) increases the size of R 0 . In particular, whenever c > 0:13, then R 0 > 1 which implies that the disease persists in the community. Figure 6(b) demonstrates the effects of ε (efficacy of preventative strategies) on R 0 . 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 R 0 as a function of c and ε. We observe that when c becomes larger or when ε is reduced, R 0 increases, implying a higher disease risk.
Numerical illustrations in Figure 7(a) show the effects of varying mosquito mortality rate, μ v on R 0 . 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 R 0 . 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.
Numerical simulation in Figure 8(a) shows the effects of varying average infectious period for symptomatic patients, 1/σ h on R 0 : 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 R 0 : We can note that whenever the 16 Computational and Mathematical Methods in Medicine probability of transmission from infectious mosquito to susceptible human is greater than 10%, the disease persists in the area, otherwise it dies out.

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 Figure 9 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.
To gain insights on the role of asymptomatic patients on chikungunya dynamics, we simulated model (18) with different values of f h (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 f h 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 f h , the increase is not highly significant compared to the other populations.
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.

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.

Data Availability
Data used in the model fitting can be obtained free of charge from doi:10.3390/mca24010006. However, the data for parameters used in the model simulation is obtained from the literature, and the references are cited in the manuscript.