Mathematical Modeling and Characterization of the Spread of Chikungunya in Colombia

The Chikungunya virus is the cause of an emerging disease in Asia and Africa, and also in America, where the virus was first detected in 2006. In this paper, we present a mathematical model of the Chikungunya epidemic at the population level that incorporates the transmission vector. The epidemic threshold parameter R0 for the extinction of disease is computed using the method of the next generation matrix, which allows for insights about what are the most relevant model parameters. Using Lyapunov function theory, some sufficient conditions for global stability of the the disease-free equilibrium are obtained. The proposed mathematical model of the Chikungunya epidemic is used to investigate and understand the importance of some specific model parameters and to give some explanation and understanding about the real infected cases with Chikungunya virus in Colombia for data belonging to the year 2015. In this study, we were able to estimate the value of the basic reproduction numberR0. We use bootstrapping and Markov chain Monte Carlo techniques in order to study parameters’ identifiability. Finally, important policies and insights are provided that could help government health institutions in reducing the number of cases of Chikungunya in Colombia.


Introduction
The Chikungunya virus is a type of arbovirus, so it is only transmitted by hematopoietic arthropods that become infected after biting some vertebrates.Later, the arthropods can transmit the virus to a susceptible vertebrate through a bite [1].The arthropod vectors are usually Aedes aegypti and Ae.albopictus [2][3][4].The virus has had mutations, allowing Aedes albopictus to transmit the disease [5].The Chikungunya disease originated in Sub-Saharan Africa, and it has become endemic in Africa, where there is a natural transmission between mosquitoes and primates to humans [3].The Chikungunya disease is an emerging disease in Asia.In America, it was detected in 2006, and there is an imminent risk of the virus spreading throughout South America.In the last few years, several outbreaks have occurred on the island of Reunion, in Cambodia, Comoros, Mayotte, Madagascar, Mauritius, Italy, Seychelles, and the Maldives [3,[6][7][8].The first outbreak in Europe was in the warm northeast region of Italy in July 2007.Probably, the virus came from Kerala (India), where the disease was at its highest peak [7].The Chikungunya disease causes arthritis, fever, and pain of the joints.Symptoms of chikungunya are generally resolved within 7-10 days, but some patients are plagued with chronic arthralgia, which could persist for months or years [9][10][11].There is no vaccine for the Chikungunya virus at this moment that could be used to restrict and control the transmission of the disease [6,12].Moreover, no effective drug is available for human use for any alphavirus, although analgesics and non-steroidal anti-inflammatory drugs can provide relief from symptoms [12].The incubation of the virus lasts between 5-12 days and the infectious stage between 5-15 days, in both primates and humans [13,14].
As mentioned above, the Chikungunya virus is spreading around the world, and due to world climate change, it is expected that more regions are going to be reached by it.For instance, in [15], the authors presented results regarding the risk of Zika and Chikungunya virus transmission in human population centers of the eastern United States.In this way, in order to better understand the dynamics of how the Chikungunya virus is transmitted, we propose and analyze a mathematical model given by a system of nonlinear differential equations where the populations of hosts and mosquitoes are homogeneous.This mathematical modeling approach is the standard way to study the dynamical behavior of diseases in populations from an epidemiological point of view [16][17][18][19].In particular, there is a variety of models for vector-borne diseases.For instance, in [20], the authors proposed a mathematical model for vector-borne disease with delay to consider the incubation period.In addition, in [21], the authors used optimal control to minimize the number of infections.Recently, some works extended these types of models using versatile fractional derivatives [22][23][24].
In [25], the authors proposed a deterministic mathematical model for Chikungunya infection considering that there is transmission of the virus between humans and mosquitoes.The authors used two infected human subpopulations designated as symptomatic and asymptomatic to classify the humans responsible for transmitting the virus.Additionally, they considered the subpopulation of humans carrying the virus, but had no possibility of spreading it.In this paper, the authors demonstrated the influence of humans on the infection of the latency period.They also remarked about the necessity of fitting the model to real data so that it will be useful in controlling the spread of the virus.Another interesting work was presented in [26], where the authors proposed a mathematical model of three age-structured transmissions of Chikungunya virus.The authors divided the human population into juvenile, adult, and senior subpopulations.
In [6], the authors proposed a stochastic mathematical model for a rural region in Cambodia, considering the fact that a stochastic model fits data better in small populations.These authors also used a subgroup of latent human and mosquito individuals and additionally included a larvae subgroup in the mosquito population.An innovation in that work is the introduction of three subclasses of infected humans.Other mathematical models for the Chikungunya virus propagation that include the latency stage were given in [27][28][29].
In the present work, our aim is to understand and explain some dynamics regarding the prevalence of Chikungunya infection in Colombia.The mathematical model is important since in practice, there are few economic resources to tackle or fight the Chikungunya virus' spread.Thus, the model allows one to decide on the best or most convenient health policies.In the proposed model, we consider a chronic subpopulation that has not been considered previously.The people in this particular class cannot transmit the disease, but have some types of chronic rheumatic symptoms.The main reason to consider this class is that health institutions are interested in the number of chronic cases of Chikungunya and its evolution.In addition, the latency stage for the human and mosquito populations is also introduced in the mathematical model, since this stage is observed in the real world.Although the virus may also spread to other vertebrate populations such as primates [30][31][32][33], reservoirs different from humans are not included in the analysis, since in Colombia, known cases are all human.However, future work might include reservoirs different from humans and would increase the complexity of the mathematical model and of the fitting process to the real data.
Numerical simulations are performed to support the theoretical results.In addition, the proposed mathematical model of Chikungunya is used to explain and understand the infected cases with the Chikungunya virus in Colombia.In particular, we use the proposed Chikungunya mathematical model to perform a fitting process to real data of Colombia.Additionally, we used bootstrapping and Markov chain Monte Carlo techniques in order to do analysis of the parameters' identifiability [34][35][36][37][38][39][40].Finally, important policies and insights are provided that could help government health institutions in lowering the infected cases with the Chikungunya virus in Colombia.
The paper is organized as follows.In Section 2, the mathematical model of Chikungunya is presented together with a set of definitions and basic underlying hypotheses.The proposed Chikungunya mathematical model is analyzed in Section 3. The fitting process and numerical simulations are presented in Section 4. Next, the parameter estimation of the model for a Colombia case using Markov chain Monte Carlo, as well as bootstrapping is performed in Section 5.The last part, Section 6, is devoted to the discussion and conclusions.

Mathematical Model
In this section, we set out a continuous mathematical model for the transmission and evolution of the Chikungunya infection in the human and mosquito populations.In the proposed mathematical model, vertical transmission is not included because the number of cases is small when compared to the total number of infected cases.However, if the virus prevalence increases in the female population, it may be necessary to include it in the proposed model.Furthermore, we do not consider a vaccinated class since there are no vaccines on the market [12,41].
The proposed model of the Chikungunya virus transmission dynamics incorporates a cross-transmission between the human and vector populations.In particular, it is assumed that the Chikungunya virus spreads by the effective contact between a mosquito infected with a human susceptible, and vice versa.This contact depends on different environmental factors.Some factors to consider are: the weather, the temperature, the altitude, and the mosquito bite rate.Varying these values will generate different degrees of the probability of the transmission of the disease.Here, we will not include seasonal effects or variability in the populations.Using a population-based approach of an epidemiological type, the population of humans is divided into five groups: susceptibles, latents, infected, recovered, and the ones with chronic rheumatic symptoms.It is important to remark that the people with chronic rheumatic symptoms do not have the Chikungunya virus, just some symptoms related to it, since they were infected previously.In addition, the population of the vector is divided into three groups: susceptible vector, latent vector, and infected vector.The resulting mathematical model is a nonlinear system of eight ordinary differential equations, which is analyzed to find the equilibrium points and their stability, including the well-known epidemic threshold parameter known as the basic reproduction number, R 0 .We estimate some of the unknown key epidemiological parameters of the model from real data, which allow us to compute R 0 , defined as the average number of secondary cases generated by a typical infectious individual in a fully-susceptible population.Moreover, we are able to compute approximately how many individuals are infected during an outbreak.All these estimates can help assist with outbreak planning, assessment of health strategies, and the design of future research regarding Chikungunya infection transmission.
Thus, following the basic ideas and structure of mathematical modeling in epidemiology, the Chikungunya model will be developed under the following basic hypotheses [16,17]:

•
The total population of humans N h (t) is divided into five subpopulations: humans who may become infected (susceptible S h (t)), humans exposed, but still not infected due to the existence of an incubation period of the virus (latent E h (t)), humans infected by the Chikungunya virus and that develop the disease (infected I h (t)), humans who have recovered from the Chikungunya infection (recovered R h (t)), and humans who have the disease chronically (chronic C h (t)).

•
The parameter µ h is the birth rate of humans.The birth rate µ h is assumed equal to natural death d h .

•
The mortality rate increase due to the disease is a real fact.However, since this rate is small in comparison with other rates and is not going to affect the dynamics, we assume that = 0.

•
The total population of mosquitoes N v (t) is divided into three subpopulations: mosquitoes who may become infected (susceptible S v (t)), mosquitoes in a latent stage (latent E v (t)), and mosquitoes currently infected or spreading the Chikungunya virus (infected I v (t)).

•
The parameter µ v is the birth rate of the mosquitoes, and it is assumed equal to the death rate d v .

•
A susceptible human can transit to the latent subpopulation E h (t) because of an effective transmission due to a bite of an infected mosquito at a rate of β 1 .

•
A susceptible mosquito can be infected if there exists an effective transmission when it bites an infected human, at a rate β 2 .

•
A fraction α of the latent humans passes to infected by the virus.

•
A fraction γ of the infected humans recovers, i.e., they do not have the disease anymore.

•
A fraction ρ of the recovered humans moves to the chronic class.

•
A fraction φ of the latent mosquitoes goes through to infected mosquitoes.

•
Homogeneous mixing is assumed, i.e., all susceptible humans have the same probability of being infected and all susceptible mosquitoes have the same probability of being infected.
It is important to notice that the parameter β 1 depends on two different parameter: the bite rate of an infected mosquito on susceptible people and the probability per bite to transmit the virus from the mosquito to the human.On the other hand, β 2 depends on the bite rate of a susceptible mosquito on an infected human and the probability per bite to transmit the virus from the human to the vector.In reality, these parameters are very difficult to estimate due to different environment conditions.Thus, the values of these parameters are most likely to vary from one region to another.In fact, the values might be very different for rural and urban areas.We assume that the natural birth rate is equal to the natural death rate for the human population since we are interested in fitting the model to the real data of the year 2015 and study the identifiability of the parameters.This is considered a short period where the births can be balanced with deaths.For instance, in [25], the authors did not even consider the births and deaths, even though these might affect the dynamics.For mosquitoes, there are no reliable data; thus, we adopted a conservative approach similar to the one used in [25].The model could be expanded to consider the variation of mosquito populations.This can be done adding a seasonal forcing function [42], but will introduce additional parameters to the model.Furthermore, it will compromise the identifiability more, which currently is an important issue.
Under the above hypotheses, the following diagram illustrates the interactions of the Chikungunya infection in human and mosquito populations; see Figure 1.The corresponding mathematical model is given by the following system of ordinary differential equations: Adding the first five equations, one gets: and therefore, the human population is constant.Analogously, by adding the last three equations, we have that: and likewise, the mosquito populations is constant.Thus, we will denote N h (t) = N h and N v (t) = N v , without dependency on time.
All parameters in this model are non-negative.It is easy to prove that the system (1) is well-posed, in the sense that if the initial data ( Ṡh (0), Ėh (0), İh (0), Ṙh (0), Ċh (0), Ṡv (0), Ėv (0), İv (0) ) are in the region R 8 + , then the solutions will be defined for all time t ≥ 0 and remain in this region.Normalizing the human and mosquito populations, and: Using the assumptions µ h = d h and µ v = d v , one can obtain the following system that describes the dynamics of Chikungunya in each class: where Now, we can define a more specific region, 8 }, and the solutions will be defined for all time t ≥ 0 and remain in this region Ω.

Equilibrium Points and Local Stability of the Chikungunya Mathematical Model
Setting the right-hand side of Equation ( 2) equal to zero, we can find the equilibrium points of the model.
The first point is denoted as the disease-free equilibrium (DFE), and the other is called the endemic equilibrium (EE).Thus, we obtain the equilibrium point DFE = (1, 0, 0, 0, 0, 1, 0, 0).On the other hand, the mathematical expression of the other equilibrium , is extremely long.Therefore, we will present this endemic equilibrium (EE) in terms of I * h .It is important to mention that our main interest is to obtain the conditions under which the population is free of disease.We will discuss later the endemic equilibria and under which conditions they exist.
For infectious diseases, the basic reproduction number is denoted as R 0 and is one of the most useful threshold parameters.It can help determine whether or not an infectious disease will spread throughout a population [43,44].R 0 also has been used for social epidemics [18,45].
In order to compute the basic reproduction number R 0 for the mathematical model ( 2), we apply the next generation technique as proposed in [43,44,46].The infectious classes in this model are e h (t), i h (t), e v (t), and i v (t).Thus, at the DFE, one gets, , and: where the matrix F is related to the rate of increase of secondary infections and V to the rate of the disease progression, death, and recovery.The next generation matrix, K = FV −1 , is non-negative, and therefore, it has a non-negative eigenvalue R 0 = ρ(FV −1 ) and a non-negative eigenvector ω associated with R 0 [43,46].There are no other eigenvalues of the matrix K with modulus greater than R 0 .Therefore, applying the next generation matrix, we obtain: or in a more compact form: Thus, we obtain the following theorem [43,44]: Theorem 1.The disease-free equilibrium DFE is locally asymptotically stable for R 0 < 1 and unstable for R 0 > 1.
This analysis shows that when the threshold parameter R 0 is less than unity, the disease disappears, while the opposite results show an endemic equilibrium in the human population.
Notice that the basic reproduction number R 0 has in the numerator all the factors that contribute to having an epidemic outbreak.For instance, we have the Chikungunya transmissibility β 1 from an infected vector to a susceptible human and the Chikungunya transmissibility β 2 from an infected human to a susceptible vector.In addition, we have the incubation rates of the virus for humans and mosquitoes α and φ, respectively.Notice that both chronic and recovered classes have a limited role in the Chikungunya dynamics, since their equations are in some sense uncoupled from the rest of the model and do not play a role in the threshold parameter R 0 .

Endemic Equilibria
In order to find the endemic equilibrium point, we set the right-hand side of Equation ( 2) equal to zero, and at least one of the infected components of the model ( 2) is different than zero.
represent any arbitrary endemic equilibrium point.Solving the equations for the steady state, we obtain: .
Since we assume that I * h = 0, then we can substitute S * h and I * v in the first equation of Model (2) at the steady state.Thus, after some calculations, we get, Then, solving for I * h , one gets: Using the threshold parameter R 0 , we can rewrite the value of I * h as: Clearly, if R 0 < 1, then I * h is negative.Thus, in order to have a positive realistic endemic equilibrium point for I * h , the threshold parameter R 0 must be greater than one.
Theorem 2. The disease-free equilibrium DFE is globally asymptotically stable on Ω for R 0 < 1.
Proof.In order to establish the global stability of the disease-free equilibrium point DFE = (1, 0, 0, 0, 0, 1, 0, 0), we are going to construct a Lyapunov function V(t) [47,48].The main idea is to find some positive constant weights W i such that dV(t) dt < 0 if R 0 < 1.Let us define the following Lyapunov function: Taking the derivative with respect to time along the solutions of the model ( 2), one gets: Rearranging and using A = (α + d h ), B = (d h + γ), and C = (d v + φ), one gets: Then, we have several options in order to obtain , we can set W 3 = 1, and then, the coefficient of i h will be negative if R 0 < 1.Then, we can set the other weights in order to cancel the other terms.Thus, we can set Therefore, we have only: 2) is uniformly persistent.
Notice that the uniform persistence and the positive invariance of the compact set Ω imply the existence of the endemic equilibrium EE of (2).
In order to establish the global stability of the endemic equilibrium , we can construct a Lyapunov function.For instance, a general form of Lyapunov functions used in the literature of mathematical biology is ), originally from the first integral of a Lotka-Volterra system [48].
Theorem 4. If R 0 > 1, then the EE point is globally asymptotically stable on Ω.
Proof.The main idea is to find some positive constant weights W i such that Let us define the following Lyapunov function: Taking the derivative with respect to time along the solutions of the model ( 2), one gets: Using the equations of Model (2), we obtain: , and and using the information regarding the EE point, one gets: Doing some rearrangement, we obtain: Since the arithmetic mean is greater than or equal to the geometric mean, we have: Thus, V(t) ≤ 0 for all (s h , e h , i h , r h , c h , s v , e v , i v ) ∈ Ω, and the strict equality V(t) = 0 holds only for (s Then, the EE point is globally asymptotically stable whenever R 0 > 1.

Numerical Simulation
In this section, a set of numerical simulations using the mathematical model of Chikungunya ( 2) is performed in order to support the presented theoretical results and validate the importance of the threshold parameter R 0 .Those simulations also help us to better understand the relationships among the different groups of the human and vector populations, the parameters, and the dynamics of Chikungunya at the population level.The numerical simulations are presented using the proportions of the subpopulations in Model (2) in order to observe easily the dynamics of Chikungunya.In addition, the simulation results are presented using the parameter values presented in Table 1, which corresponds approximately to the current Colombian scenario.All simulations were done using an adaptive Runge-Kutta-Fehlberg method of order four [49].2) when the threshold parameter R 0 < 1 in order to corroborate that the infected populations vanish and to observe the dynamics of the susceptible recovered and chronic populations.We run the simulations for a time long enough so we can observe both the transient dynamics and the steady states of the system.
For this first numerical simulation, the parameter values used are those presented in Table 1, but β 1 = 1/720.In this particular scenario, the numerical value of β 1 is chosen such that R 0 < 1.
As can deducted directly from the Chikungunya model ( 2), the units of all these parameters are in days −1 , except β i , which measures the effective contacts per day, i.e., the total number of contacts, effective or not, per unit day, multiplied by the risk of infection with Chikungunya virus.
In Figure 2, it can be observed that for R 0 ≈ 0.93, the infected population I(t) dies out, and the susceptible population approaches the disease-free steady state value S d = 1.Moreover, the steady state is as expected the disease-free equilibrium DFE = (1, 0, 0, 0, 0, 1, 0, 0).In this way, if health institutions have to reduce the Chikungunya prevalence in the population, then they need to introduce changes in the health policies that affect the corresponding parameters related to R 0 in order to reduce its value.

Numerical Simulations for R 0 > 1
The second numerical simulation of the Chikungunya model ( 2) is when the threshold parameter R 0 > 1.For this particular scenario, the infected populations persist over time, as well as the susceptible and recovered populations.As in the previous scenario, we run the simulations for a time long enough so we can observe both the transient dynamics and the steady states of the system.
In Figure 3, it can be observed that for R 0 ≈ 1.1977, the infected subpopulations i(t) and i v (t) persist and approach an endemic steady state value.However, notice that despite the time-invariant parameter values, the infected subpopulations oscillate at the beginning and then reach a steady state.Thus, some slight changes of the parameter values due to environmental events may change the outcome dynamics.It is important to remark that if health institutions want to reduce the Chikungunya prevalence in the population, then they need to introduce changes in the health policies that affect the corresponding parameters related to R 0 in order to reduce its value.

Sensitivity Analysis of the Transmission Parameters
Here, we investigate the relevance of each transmission parameter of the Chikungunya model (2) on the dynamics of the subpopulation of humans infected by the virus.To answer this key question, we need to compare the impact of those parameters on the outcome of the infected.As has been mentioned, the threshold parameter R 0 plays a determinant role in the different subpopulations of the model.
In the previous numerical simulation, we can see in Figure 4 that varying the transmission parameter β 2 allows us to observe the effect of this parameter on the dynamics transmission of the disease in the human and mosquito population.In all these simulated scenarios, the human population is exposed to the mosquitoes infected by the Chikungunya virus.Based on the fact that the human population exposed to the Aedes aegypti in Colombia is 19 million [51] and assuming different values of the transmission parameter β 2 , we can roughly estimate the number of infected humans by Chikungunya virus, as given in Table 2.

Estimation of Parameters for the Colombian Scenario
The dynamics of systems are, in general, too complex to allow intuitive predictions and require the support of mathematical modeling for quantitative assessments and a reliable understanding of the system functioning [56].Moreover, one of the most difficult tasks of mathematical modeling is the estimation of model parameters.
The proposed Chikungunya mathematical model (2) establishes mathematical relationships among the eight different sub-populations and allows for the flux of individuals between sub-populations.Therefore, several parameters that regulate these relationships and fluxes appear in the model.However, the estimation of these model parameters is not a straightforward task.

Fitting Algorithm
In order to use the Chikungunya mathematical epidemiological model (2) to simulate the dynamics of Chikungunya virus in the Colombian population, it is necessary to set the parameter values of the model.However, some of the parameter values are not accurately known or belong to different regions.Here, we rely on some parameter values that are available in scientific journals or based on medical considerations.One aim here is to explain the qualitative and quantitative behavior of the Chikungunya infection dynamics at the population level in Colombia for the year 2015.
The Chikungunya disease data were collected using available cumulative data from a national health institution in Colombia.In Table 3, we can see the seroprevalence of Chikungunya infection for different weeks of 2015 in Colombia.In order to adjust the Chikungunya mathematical model ( 2) to the time-series data of Chikungunya seroprevalence in Colombia, the only parameters to be estimated by a fitting process to real data are the Chikungunya transmissibility β 1 from infected vector to susceptible human, Chikungunya transmissibility β 2 from an infected human to a susceptible vector, and the initial infected human and mosquito proportions.The parameters β depend on the number of bites per unit of time and the transmission probability per bite [57].
In regard to initial conditions for the year 2015, we have several assumptions due to the lack of real data for the vector and human populations.However, we think that the assumptions are biological plausible since we take into account some available data in conjunction with parameter values that are known from the related literature.As we have mentioned before, we used the proportions of the subpopulations in the Chikungunya mathematical model (2).It is important to point out that the data presented in Table 3 are related to individual cases, but we transform these data to proportional values related to the real demographic data of Colombia.We use a total population of 19,471,223, since this is the number of inhabitants that live within the 2200 m where the mosquitoes that transmit the Chikungunya virus live [51,58].The following initial conditions based on the proportion of infected humans i 0 (i h (0)) are assumed: e 0 = i 0 /2 based on the fact that the latent or exposed stage is half of the infected one; r 0 = i 0 /40 based on the fact that it takes 300 days for a recovered person to start suffering chronic symptoms and with a probability of approximately 5% [54,55].Finally, we set the initial condition for the proportion of chronic individuals as c 0 = 0.05 × r 0 , and the initial condition for the susceptible population is just found using the relation s 0 = 1 − e 0 − i 0 − r 0 − c 0 .For the vector population, we do not have any real data.We will assume that the proportion of vectors at the eclipse phase is a tenth of the infected vectors based on the duration of the latent and infectious stages of the mosquitoes.Here, the simulation interval period of 2015 has been chosen according to the available Chikungunya seroprevalence data in Colombia.
The fitting process to adjust the Chikungunya mathematical model ( 2) to the time-series data of Chikungunya in Colombia corresponding to the year 2015 is done minimizing the sum of squared errors (SSR).We use two algorithms to find the minimum SSR.We initially used a genetic algorithm [59], which performs a broad search of the parameter space and is less dependent on the initial guess.Once the genetic algorithm had found a good fit, these parameters were used as the initial guess for the trust-region-reflective and interior point algorithms [60,61], which search a more localized region of the parameter space.The use of several different algorithms increases the probability of finding the global minimum for the SSR.
In order to compute the best fitting, we carried out computations, and we implemented the SSR function, where β 1 , β 2 , and i v (0) are variables such that: 1.

3.
Compute the sum of square of the difference between îh (0), îh (1), îh (2),. .., îh (51), and infectious data in Table 1.This function F returns the sum of squared errors (SSR), where for the Colombia data are given by: Find a global minimum for the the sum of squared errors (SSR) using genetic, trust-region-reflective, and interior point algorithms.
The function F takes values in R 3 and returns a positive real number, the SSR that measures the closeness of the scaled infectious population (i h (t)), provided by the model, to time-series data.In order to ensure that parameter estimates are biologically realistic, we placed bounds on some of the parameters.Thus, for instance, we can obtain values for i v (0) only from the interval (0, 1) or positive values for β 1 and β 2 .

Numerical Simulation of the Chikungunya Mathematical Model
In this section, numerical results for the solution of the Chikungunya mathematical model (2) are presented.Since no analytic solution to the nonlinear fractional system ( 2) is available, we use a Runge-Kutta-type method to compute the solution numerically.The data from Colombia related to Chikungunya are from Weeks 1-52.It is important to remark that the data only show a small increasing interval for infected cases and the full decreasing interval.This is due to the fact that at the beginning of the Chikungunya epidemic, the health institutions diagnosed the cases using just symptoms instead of the lab tests that confirm a real infected case.Therefore, health institutions in Colombia disregarded most of the first increasing period in order to avoid inaccurate reports of Chikungunya.Notice that Dengue virus causes similar symptoms to the Chikungunya virus in human populations.Thus, lab tests are necessary to differentiate between both epidemics.
In order to fit the Chikungunya mathematical model ( 2) to the time-series data of Chikungunya cases in Colombia, we need to estimate the parameters β 1 , β 2 , and i v (0).As a first approach, we set the initial proportion of latent and infected humans based on the real data.Thus, we only need to estimate initially β 1 , β 2 , and i v (0).In other works, different assumptions regarding the initial infected and latent populations for humans and vectors have been assumed.For instance, in [62], the authors assumed that the initial numbers of latent and infectious people were equal, and analogously for the mosquito population.
The confirmed cases of Chikungunya from the year 2015 in Colombia and the best model fit along with the best fit parameter values are shown in Figure 5.It can be seen graphically that the Chikungunya mathematical model (2) produces a good adjustment to the real data and predicts that the epidemic will disappear in Colombia due to a R 0 = 0.79.In addition, the epidemic peak of the mathematical model approximates well the peak of the real data.It is important to mention that the Chikungunya mathematical model (2) generates epidemic data that fit well in terms of the SSR.The Chikungunya mathematical model (2) gives a smoother curve due to its deterministic nature, as was expected.In order to catch the natural irregularity of the real data, it would be necessary to introduce stochastic factors to the model, which allows one to obtain a more accurate fitting [63].However, introducing stochastic or temporal factors would require more detailed information regarding the dynamics of the population in Colombia, and the complexity of the model would increase.It is important to mention that the irregularity of the real data has been observed in many other studies related to other diseases, and can be explained due to different reasons such as weather, under-reporting, the stochastic nature of the virus diffusion, spatial effects, or even the heterogeneity of the human and mosquito population [19].In Figure 6, the long-term dynamics of the humans with chronic rheumatoid symptoms can be observed.It can be seen that this class population will reach a peak and then start to decrease due to the natural death rate and the decay of new infected humans over the long term.
Additionally, we considered different scenarios fixing the initial conditions for the mosquito infected population i v (0) and fitting the initial infected human population i h (0).However, the most important difference with the previous best fit scenario is that the values of the parameters β 2 and β 1 varied in such way that they compensated each other in order to be able to fit the real data level of infected humans.This compensation makes total sense since the transmission parameters β s are going to depend on the initial infected vector population i v (0).The smaller i v (0) is, the greater should be the transmission parameter from vector to humans in order to reproduce the real number of infected cases.However, we notice that the basic reproduction number R 0 stays approximately constant for all the scenarios.Thus, from the results, we obtained a robust numerical value for the basic reproduction number R 0 , regardless of the initial proportion of infected mosquitoes.
Finally, we present the numerical simulations for larger values of i v (0).We vary i v (0) from 0.2 (20%) to 0.5 (50%) in order to observe other potential realistic scenarios.The idea is to test hypotheses that consider larger values for the initial proportion of infected vectors i v (0).The best model fits can be seen in Figure 7.It can be observed graphically that the mathematical model (2) does not fit the real data well for large values of i v (0).Based on these results and the real data available, we can infer that these large value scenarios are unlikely for Colombia in 2015.Moreover, notice that for larger values of the initial proportion of infected vectors i v (0), the model fits to the real data deteriorate.Thus, we have highlighted some likely scenarios for the parameter values of β 1 , β 2 , i v (0), and i 0 that can help to analyze some potential scenarios beyond the data and with some intervention strategies.Moreover, this allows us to compare with different regions where other parameter values have been obtained and then raise some questions regarding the explanations of these differences.In the next section, we will provide further analysis regarding the identifiability of the parameters, the initial conditions, and the basic reproduction number R 0 .

Identifiability of the Parameters
As we mentioned before, the parameters β 1 and β 2 are not fully identifiable.Moreover, we do not have real data regarding initial conditions for the different subpopulations, but based on our assumptions, all of them can be set in terms of the initial conditions for the infected human and vector populations.Here, we use two approaches: for the first one, we fixed the initial condition of the infected humans to the real data, and consequently, we minimized the function F in R 3 ; the second approach was just to leave the initial condition of infected humans as a parameter, and then, we minimized the function F in R 4 .
We now extend our parameters' identifiability analysis, by incorporating two numerical techniques.These techniques are bootstrapping [34,64] and Markov chain Monte Carlo [35,36,39], which allow us to introduce further information regarding the identifiability of the parameters and corroborate the identifiability of the basic reproduction number R 0 .It is recommended to ensure that the optimum parameter values of the model can be uniquely determined by the available real data.For non-linear-based models, the issue of identifiability is not straightforward.
An important numerical tool often used to assess the uncertainty in estimated values is bootstrapping [34,64].The bootstrapping process begins with the generation of artificially-generated datasets, which are created by sampling the best fit curve and adding error such that the SSR is equal to the SSR of the original data.The mathematical model ( 2) is fitted to each of the surrogated datasets, leading to new parameter estimates.A total of 3000 bootstrap replicates were performed using the real data.In this way, we obtained estimates of the distribution of each parameter, as can be seen in Figure 8. Roughly, it can be seen from the histograms that all the parameter estimates follow a Gaussian distribution.These distributions were used to give the 95% confidence intervals, which are shown in Table 4.This information is useful to have a measure of error in the parameter estimates.There is some small skew present in the histograms of parameters T h and T M , which can be explained by their correlation.However, we are more interested in the threshold parameter R 0 , which as can be seen in Figure 9, is relatively stable around R 0 = 0.78.
The bootstrapping method allows us to find correlations between estimated parameters.In Figure 9, several two-parameter scatter plots for each of the 3000 bootstrap replicates can be observed.In order to have two uncorrelated parameters, the scatter plot should be roughly a circle, with no clear relationship between the two parameters.For the parameter i v (0), the plots show a little bit of correlation, and this was in some way expected since i v (0) is a factor on how strongly the epidemic starts.There is a clear correlation between the transmission parameters β 1 and β 2 .This correlation in some way is expected since both transmission parameters are related to the exposure of humans to mosquitoes and vice versa.However, it is important to remark that we are more interested in their product since the epidemic threshold parameter R 0 is proportional to this product.These results suggest that the threshold parameter R 0 of the mathematical model ( 2)) is identifiable.Thus, we can use the model with caution to describe the dynamics of the spread of Chikungunya in the population of Colombia in the year 2016.The second numerical technique is based on Markov chain Monte Carlo [35,36].In particular, we used a stable algorithm of the affine-invariant ensemble sampler for Markov chain Monte Carlo (MCMC) proposed in [65].An algorithm that is affine invariant performs equally well under all linear transformations; it will therefore be insensitive to covariances among parameters [66].Markov chain Monte Carlo is designed to sample from-and thereby provide sampling approximations to-the posterior PDF efficiently [35,66].In the algorithm that is used here, we draw samples from a multivariate Gaussian density.It is necessary to set up the specific values of the hyperparameters in three dimensions (β 1 , β 2 , i v (0)).In addition, we need to decide how many walkers there are, which are in some way independent paths to reach the maximum of the likelihood function [65].There are different options to initialize each of the walkers.One of the best techniques is to start in a small ball around the a priori preferred position [66].It is important to mention that the walkers will spread out and explore the whole space for the parameters.The main goal is to find the maximum of the likelihood function, which is a Gaussian where the variance is underestimated by some fractional amount.We can use Markov Chain Monte Carlo to obtain estimates and confidence intervals for each of the parameters β 1 , β 2 , and i v (0).Moreover, we can observe some possible correlations among the parameters.
In Figure 10, all the one-and two-dimensional projections of the posterior probability distributions of the parameters can be seen.This allows us to observe all of the covariances between parameters.In addition, the histograms show the marginalized distribution for a parameter or a set of parameters using the results of Markov chain Monte Carlo.Thus, Figure 10 shows the marginalized distribution for each parameter independently in the histograms along the diagonal and then the marginalized two-dimensional distributions in the other panels [66].Notice that the maximum likelihood function profile is also presented.We can say that in order to have two uncorrelated parameters, the level curves' plot should be close to circles, with no clear relationship between the two parameters.For the parameter i v (0), the plots show less correlation in comparison with the ones with respect to β 1 and β 2 , and this was in some way expected since i v (0) is a factor of how strongly the epidemic starts.There is a clear correlation between the transmission parameters β 1 and β 2 , as we also observed with the bootstrapping results.In fact, this relationship has the same pattern.This correlation in some way is expected since both transmission parameters are related to the exposure of humans to mosquitoes and vice versa.However, it is important to remark that we are more interested in their product since the epidemic threshold parameter R 0 is proportional to this product, as we have mentioned before.These results agree with previous results regarding identifiability for epidemic models with vectors [67].

Conclusions
We present a mathematical model of the spread of the Chikungunya disease at the population level that incorporates the transmission vector by including cross-transmission between the human and vector populations.The proposed model includes a chronic subpopulation, which to the best of our knowledge has not been considered in other mathematical models.We determined the epidemic threshold parameter R 0 for the extinction of disease using the method of the next generation matrix.Using Lyapunov function theory, some sufficient conditions for the global stability of the the disease-free equilibrium were obtained.Based on this parameter, we found the parameters that affect the basic reproduction number R 0 and therefore what would be the best policies to control the spread of the Chikungunya disease.We verified that when the threshold parameter R 0 is less than unity, the disease disappears, while when threshold parameter values are larger than one, the disease persists in the population.Numerical simulations were presented to support the established theoretical results.
Using the proposed mathematical model of Chikungunya diffusion, we were able to analyze the dynamics of infection during the 2015 outbreak in Colombia.In particular, we estimated the numerical values of the epidemiological parameters β 1 , β 2 , and the reproduction number R 0 .Based on numerical results of the model of Chikungunya, we were able to better explain and understand the variation of the number of infected cases with the Chikungunya virus in Colombia.We found that the transmission parameters β 1 and β 2 were highly correlated.
Our estimated values for the reproduction number R 0 ranged from 0.74-0.84across the three different scenarios proposed and also from the results of bootstrapping and Markov chain Monte Carlo methods.However, notice that rain seasons have been changing recently due to climate change, and that could affect some of the parameter values, which is one limitation of this model.The fitting of the model to the observed weekly reports during the 2015 outbreak in Colombia is relatively good, despite the natural irregularity of the real data.This irregularity of the data has been observed in many other studies related to other diseases and can be due to different reasons such as weather, under-reporting, the stochastic nature of the virus diffusion, spatial effects, or even the heterogeneity of the human and mosquito population.All these factors were not considered here explicitly and are other limitations of this work.Other limitations are the fact that the initial proportion of the different subpopulations is unknown.More complex models are necessary to study all the aforementioned factors, and future works might consider these factors or at least some of them.Another important aspect that we found in this study, based on the previous numerical simulations of the different scenarios, is that the initial infected population of mosquitoes is likely no greater than 20% for the year 2015 in Colombia.To the best of our knowledge, we do not know if there are real data regarding this particular description.
In this paper, we performed identifiability analysis in order to estimate some parameters of the model and found correlations among the parameters of the model.We used the bootstrapping technique and Markov chain Monte Carlo in order to assess the identifiability of the parameters of the model.We found out that the parameters were not fully identifiable with the prevalence of the real data that we had.However, we were able to identify that the reproduction number R 0 ranged from 0.74-0.84.It is important to notice that this result depends on some fixed parameter values such as the duration of the infectious stage and the life-span of the mosquitoes.Changing values for these parameters can change the aforementioned range for the reproduction number R 0 .However, regardless of some potential changes, we presented a study with a methodology to deal with these types of epidemic models when real data are available.
Finally, we can conclude that in order to reduce the Chikungunya disease diffusion, it is important first to decrease the transmission parameters.For instance, people could use clothes covering most of the body to avoid mosquito bits and could also use repellents.Any action that could reduce the infected mosquitoes' bites will affect the reproduction number R 0 .The aforementioned actions would reduce the transmission of the virus in both directions.Another way to reduce the value of the reproduction number R 0 is to increase the mortality of the mosquito population.Therefore, a health policy that could be implemented is to increase this mortality.Other parameters that are included in R 0 such as incubation or infected period are natural characteristics of the Chikungunya virus and are difficult to modify.

Figure 1 .
Figure 1.Dynamics of the Chikungunya virus with transmission vector.The boxes represent the subpopulation and the arrows the transition between the subpopulations.Arrows are labeled by their corresponding model parameters.

Figure 4 .
Figure 4.The dynamics of the different subpopulations are sensitive to the changes of the parameter β 2 .

Figure 6 .
Figure 6.Dynamic of the chronic infected individual humans using the fit of the Chikungunya mathematical model (2) to the time-series data of Chikungunya in Colombia corresponding to the year 2015.

Figure 7 .
Figure 7. Best fits of the Chikungunya mathematical model (2) to the time-series data of Chikungunya in Colombia corresponding to the year 2015.The initial proportion of infected vectors i v (0) varies from 0.2 (20%) to 0.5 (50%).

Figure 9 .
Figure 9. Identifiability assessment of the mathematical model (2) fit to the prevalence data of Chikungunya in the population of Colombia in the year 2016.Correlation plots generated with parameter estimates from bootstrap fits.Correlation between the parameters β 1 , β 2 , i M (0), and R 0 .

Table 2 .
Variation of the β 2 parameter and its effects on the infected human population.

Table 3 .
Data provided by the National Institute of Health-SIVIGILA, Colombia.The first row corresponds to the number of the week of the year 2015.The second row shows the number of infectious individuals for whom Chikungunya was detected in each week because those people went to see a doctor and it was reported.
Best fit of the Chikungunya mathematical model (2) to the time-series data of Chikungunya in Colombia corresponding to the year 2015.Red points give the real data, and the blue line shows the best model fit.The best fit parameter values are given in the table.SSR, sum of squared errors.

Table 4 .
(2)imated parameters using the mathematical model(2)and the prevalence data of Chikungunya in the population of Colombia in the year 2016.Identifiability assessment of the mathematical model (2) fit to the prevalence data of Chikungunya in the population of Colombia in the year 2016.Correlation plots generated with parameter estimates from bootstrap fits.
The 95% confidence intervals (bootstrap fits) are given in parentheses.