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.

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, 𝕄-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: 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

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 aswhere Γ(·) is the gamma function and is defined by Γ(z) = ∫0∞e−tdt and θ > 0.

Definition 2 (see [37]).

The Caputo fractional derivative of orderθfor a functionf ∈ C([t0, +∞), ℝ)is defined bywhere 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

Lemma 1 (see [38]).

The solution to the Cauchy problem with 0 < θ < 1 and λ ∈ ℝ has the form while the solution to the problem is given by where E(·) is the Mittag-Leffler function and defined by

Lemma 2 (see [39]).

Letu(t)be a continuous function on[t0, +∞)and satisfyingwith 0 < θ < 1, (λ, μ) ∈ ℝ2, and λ ≠ 0, and t0 ≥ 0 is the initial time. Then,

Lemma 3 (see [40]).

Letx(t)be a continuous and differentiable function withx(t) ∈ ℝ+. Then, for any time instantt ≥ t0, one has

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:

The total populations of humans N(t) are subdivided into classes of susceptible S(t), exposed E(t), asymptomatic infectious individuals A(t), symptomatic infectious individuals I(t), and humans who have recovered from infection R(t). The total vector population Nv(t) is subdivided into compartments of susceptible S(t), exposed E(t), and infectious I(t). 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

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 μ with j = h, v. Further, new recruits are assumed susceptible and their recruitment rate is proportional to the population and is given by μN. Meanwhile, we assume that whenever there is effective contact between a susceptible individual and an infectious vector, disease transmission will occur at rate β. Similarly, let β 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 − η) to model the reduction of infectivity of the asymptomatic individuals

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 − ε model the effects of preventative strategies utilized by humans to minimize chances of contracting the infection (naturally or through treatment)

Exposed humans are assumed to incubate the infection for 1/α days after which a proportion f become asymptomatic infectious patients and the remainder, (1 − f), become symptomatic infectious patients. The average infectious period of asymptomatic infectious individuals is modeled by 1/γ, and symptomatic infectious individuals are assumed to be infectious for an average period of 1/σ days

Exposed vectors incubate the disease for 1/α 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: 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

Figure 1

Flow chart for chikungunya disease transmission.

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(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:

3.2. Positivity and Boundedness of Model Solutions

Theorem 1 .

The fractional order (18) has a unique solution, which remains inℝ7, and the closed setis 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)

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 n(t) = s(t) + e(t) + a(t) + i(t). By adding all the equations for the host population in (6), one gets

It follows from Lemma 3 that (38) has a solution of the form

Since E(−μt) ≥ 0, so that when n(0) ≤ 1, we have n(t) ≤ 1. Similarly, by adding all the equations for the vector population (with, n = s + e + i), one gets

Again, by Lemma 3, we have

By similar arguments as before, we have n(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 disease-free equilibrium (DFE) point (denoted by E0) and the endemic equilibrium (EE) point (denoted by E∗) point, which are, respectively, given by E0 : (s0, e0, a0, i0, s0, e0, i0) = (1, 0, 0, 0, 1, 0, 0). The endemic equilibrium (EE) point E∗ = (S∗, E∗, I∗, S∗, E∗, I∗) with with

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 2 .

IfR 0 < 1, then the DFE of system ((18)) is globally asymptotically stable inΩ, otherwise it is unstable.

By considering only the infected compartments from (18), we can write where x = (e, a, i, e, i), with F and V defined as follows:

One can verify by direct calculations that V−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.,

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  DL(t) = 0 is the singleton {ℰ0}. Therefore, by LaSalle's invariance principle [43], ℰ0 is globally asymptotically stable in Γ when R0 ≤ 1.☐

Theorem 3 .

IfR 0 > 1, then the endemic equilibriumE∗of system ((18)) is globally asymptotically stable inΩ.

To prove Theorem 8, we consider the following Lyapunov functional:

Let g(i, a) = i + ηa. 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  DU(t) ≤ 0 for all (s, e, a, i, s, e, i) ≥ 0 since s∗, e∗, a∗, i∗s∗, e∗, and i∗ 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, 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 = nh, and n = 0, 1, 2, ⋯, N ∈ ℤ+, equation (39) can be discretized as follows: where

and the predicted value ψ(t) is determined by with

The error estimate is with k ∈ ℕ and p = min(2, 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]: where

Further, we have

In addition, the quantities are derived from system (53), at the points t, 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), C(t), are given by the solution of the following equation:

In order to compute the best fitting, we implemented the function where θ, η, ϕ, ε, β are variables such that

for a given (θ, η, ϕ, ε, β), solve numerically the system of differential equation (18) and equation (55) to obtain a solution , which is an approximation of the real-world data Y(t)

Set t0 = 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 i(t); that is., , …..

Compute the root mean square (RMSE) of the difference between and observed cases. This function 𝔽 returns the root-mean-square error (RMSE) where

(4) Determine a global minimum for the RMSE using Nelder-Mead algorithm

The function 𝔽 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 e(0) = 0.0002, a = i(0) = 0.0001, s(0) = 1 − e(0) − a(0) − i(0), e(0) = 0.0002, i(0) = 0.0001, and s(0) = 1 − e(0) − i(0) and the baseline values for the model parameters are in Table 1.

Table 1

Description of parameters used in the model system (18).

Symbol	Description	Value	Units
μ h	Natural mortality rate of humans	0.00004 day −1	[24]
1/αh	Incubation period of in humans	12 (5–12) days −1	[20, 24]
c	Biting rate of human by female mosquito	0.19 (0.19–0.39) days −1	[20, 24]
f h	Proportion of asymptomatic patients	0.1 (0.1–0.3)	[20]
1/γh	Infectious period for asymptomatic patients	10 (7–10) days −1	[20]
1/σh	Infectious period for symptomatic patients	7 (7–10) days −1	[20, 24]
ϕ	The number of mosquitoes per individuals	3	Fitting
η h	Reduction of infectivity from asymptomatic patients to susceptible mosquito	0.2 (0–1.0)	Fitting
ε h	Efficacy of preventative strategies	0.3 [0–1]	Fitting
β h	Probability of infection to be transmitted from humans to mosquito per bite	0.9 (0.8–1)	Fitting
β v	Probability of infection to be transmitted from an infectious to a susceptible human per bite	0.244 (0.06–0.244)	[24, 46]
1/μv	Life expectancy of mosquitoes	30 (4–30) day −1	[19, 20]
1/αv	Incubation of mosquitoes	3 (3–7) day −1	[24, 46]
R 0	Basic reproduction number	1.5–7	[20]

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

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].

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 Theorems 7 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 number R0 with respect to model parameters. Baseline values used are in Table 1.

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 R0. In particular, if the mosquitoes increase their feeding rate by 10%, then the magnitude of R0 will also increase by 10%. Parameters β (transmission rate of the disease from humans to mosquito per bite), β (transmission rate of the disease from mosquito to human per bite), ϕ (the number of mosquitoes per individual), and μ (mortality rate of mosquitoes) have a similar elasticity index in absolute value. However, an increase in β, β, and ϕ will increase the size of R0 while an increase in μ reduces the size of R0. It is also worth noting that the efficacy of preventative strategies ε 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 Figure 6(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.

Numerical illustrations in Figure 7(a) show the effects of varying mosquito mortality rate, μ on R0. From the results, we note that if the life span of mosquitoes is less than 14 days (μ = 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, μ on R0, and (b) the number of mosquitoes per individuals, ϕ on R0.

Numerical simulation in Figure 8(a) shows the effects of varying average infectious period for symptomatic patients, 1/σ 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, β 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/σ on R0, and (b) the probability of infection to be transmitted from an infectious to a susceptible human per bite, β on R0.

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

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.

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.

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

Figure 11

Numerical results for different values of f (proportion of individuals who become asymptomatic patients upon the completion of the incubation period) over time with θ = 0.64.

To explore the impact of the preventative strategies on long-term dynamics, we simulated model (18) with different values of ε with a fixed derivative order (θ = 0.64) and the results are in Figure 12. Overall, the results show that as ε approaches unity, the disease dies out in the community. In particular, the closer ε is to unity, the shorter the time it takes for the disease to become extinct.

Figure 12

Numerical results showing the effects of varying ε on long-term dynamics of the disease.

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.