Modelling and Analyzing Virus Mutation Dynamics of Chikungunya Outbreaks

Chikungunya fever, caused by chikungunya virus (CHIKV) and transmitted to humans by infected Aedes mosquitoes, has posed a global threat in several countries in 2015. Recent outbreaks in La Réunion, Italy and China are related with a new variant of CHIKV with shorter extrinsic incubation period in contaminated mosquitoes, but the role of this new variant on the spread of chikungunya fever is unclear. We develop a mathematical model that incorporates the virus mutation dynamics in the transmission of CHIKV among mosquitoes and humans. Our numerical simulations show that a substantial virus mutation rate combined with high virus transmission probabilities from mosquito to human, could result in sustainable chikungunya fever outbreaks. Further, we apply Markov Chain Monte Carlo sampling method to fit our model to the 2007 chikungunya fever outbreak data in North-Eastern Italy where the mutant strain was detected. We conclude that the basic reproduction number might be underestimated without considering the mutation dynamics, and our estimation shows that the basic reproduction number of the 2007 Italy outbreak was \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\pmb{\mathscr{R}}}_{{\bf{0}}}$$\end{document}R0 = 2.035[95%Cl: 1.9424 - 2.1366]. Sensitivity analysis shows that the transmission rate of the mutant strain from mosquitoes to human is more influential on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\pmb{\mathscr{R}}}_{{\bf{0}}}$$\end{document}R0 than the shortened extrinsic incubation period. We conclude that the virus mutation dynamics could play an important role in the transmission of CHIKV, and there is a crucial need to better understand the mutation mechanism.

. CHIKV transmission flow chart. Subindices 1 and 2 correspond to the non-mutant and mutant strains, human populations are stratified into exposed, asymptomatically/symptomatically infected and recovered compartments. The mutation dynamics is characterized by the flow from the mosquitoes infected with the non-mutant strain to the mosquitoes with the mutant strain. Solid arrows represent the movements of population among compartments. Compartments responsible for the transmission of CHIKV are colored correspondingly to the transmission rates.

Mathematical model
We stratify the infected mosquitoes and human populations in terms of non-mutant and mutant strains to describe the mutation dynamics. For mosquitoes, the total mosquito population at time t is denoted by N t ( ) M , which includes the following epidemiological compartments: susceptible S ( ) M , exposed E ( M1 and E ) M2 , infectious I ( M1 and I ) M2 , with sub-index 1 for non-mutant and 2 for mutant strain. Susceptible mosquitoes move to the latent compartment = E i ( 1,2) Mi , upon successful exposure to the respective strain. The mosquito recruitment and death rates are denoted by λ M and µ M . Exposed mosquitoes become infectious after the latent period γ 1/ M1 or γ 1/ M2 . We assume that the virus could mutate in contaminated mosquitoes all through their life time at a rate of δ. For humans, we take into account of the asymptomatic and symptomatic infections. The total human population N H includes the following epidemiological compartments: susceptible S ( ), H exposed E ( H1 and E ), H2 asymptomatically infected I I ( , ) H N H N 1 2 , symptomatically infected I ( H1 , I ) H2 , and recovered R ( H1 , R ) H2 . We assume that a proportion φ of exposed individuals will become symptomatic after the latent period γ 1/ H , and both symptomatic and asymptomatic individuals have the same probability of transmitting the virus to mosquitoes via each bite. We assume that the virus only mutates in the mosquito population because of the selective pressure that has been demonstrated in mosquito bodies but not been found in humans 11 . It has been shown in 30 that people experienced the non-mutant CHIKV outbreak in 1975, Cambodia, obtain a lower risk of infection during the outbreak of the mutant strain in 2012. Thus we make a simplified assumption on the existence of cross-immunity between the two strains, i.e., people recovered from one strain are immune to both of the strains for life time. Following the transmission diagram shown in Fig. 1 here, b is the average biting rate of the mosquitoes. β H M 1 (β H M 2 ) is the transmission probability from human with non-mutant strain to mosquitoes (from human with mutant strain to mosquitoes), correspondingly, β M H 1 (β M H 2 ) is the transmission probability from mosquitoes with non-mutant strain to human (from mosquitoes with mutant strain to human). Symptomatic and asymptomatic individuals have the same recovery rate q. The parameter δ is the virus mutation rate (which is essentially the mutant strain adaptation rate) within the bodies of contaminated mosquitoes. CHIKV usually does not cause deaths of human, therefore, we do not consider the disease related deaths in model (1). λ H is the human recruitment rate and µ 1/ H is the average lifespan of humans. We consider model (1) with non-negative initial conditions and parameter values, and we adopt the ranges of parameters as illustrated in Table 1.
It is obvious that system (1) always has a disease free equilibrium = We follow the calculation of basic reproduction number developed in 31 where = x  E I E  I  E I I E  I I  ( ,  ,  ,  ,  , , ,  ,  ,  )   M  M  M  M  H  H  N  H  H  H  N  H  T  1  1  2  2  1  1  1  2  2   2 and The well-known result in 31 shows that the maximum real part of all eigenvalues of the matrix − F V is negative if and only if the spectral radius of the next generation matrix ρ < .
It is straightforward to obtain that  www.nature.com/scientificreports www.nature.com/scientificreports/ The eigenvalues of the matrix − FV 1 are given by the roots of the characteristic equation: The basic reproduction number 0  , defined as the average number of secondary cases arising from an average primary case in an entirely susceptible population, is the spectral radius of ρ − FV ( ) 1 31,32 ). Therefore, These explicit formula for the non-trivial components of boundary equilibria led to the following theorem.
then the mutant strain prevalent equilibrium E 2 exists.
The explicit formula for the positive equilibria with co-existence of both strains for system (1) is hard to derive algebraically. We thus seek to gain our insights from numerical simulations for the existence and stability of the positive equilibria, the boundary equilibria and possible periodic oscillations. We will mainly discuss three cases about the mutation dynamics i ( ) the infectivity from mosquitoes to humans of the mutant strain is smaller comparing with that of the non-mutant strain, i.e., β β <  Fig. 2, the non-mutant strain dominates for small mutation rates, co-existence of both strains could occur under moderate mutation rates, and mutant strain dominates for large mutation rates. Figure 3 shows that, if the mutant strain already possesses a higher transmission probability from mosquito to human, then it will dominate at any mutation rates. In the scenario of Fig. 4 where the two strains obtain same level of transmission probability, co-existence happens when there is no mutation, and mutant strain dominates as soon as the mutation rate perturbs from 0. These results indicate that the mutation dynamics in the bodies of Aedes albopictus mosquitoes could result in a progressive take-over of the non-mutant strain by the mutant strain in the human population level in the long run.

Fitting the Italy 2007 outbreak data
Aedes albopictus was introduced in 1990 to Italy and has since been established in urban settings, and was identified as the principal vector for the 2007 CHIKV outbreak in Italy. The mutant strain ECSA-V is believed to be associated to this outbreak 33 , moreover, the co-existence of mutant and nonmutant CHIKV has been reported in this outbreak 34 . We thus use our model to fit this outbreak data and investigate the role played by the mutation dynamics. Figure 5 shows the reported case numbers between June 23rd and September 14th, 2007, where the blue bars up to day August 22nd represent the initial outbreak phase before the start of public health implementations. Since the control measures will greatly alter the model parameters, we fit the outbreak data from June 23rd to August 22nd.
Due the short epidemic time scale in comparison to the demographic time scale, we do not consider the natural growth of human population and set λ µ (1). Further, the outbreak happened mainly in two small villages Castiglione di Cervia and Castiglione di Ravenna (Ravenna province, in northeastern Italy). Therefore, we assume the human population as a constant N H . This simplifies (1) to model (S1) as shown in the Supplementary material.
For the purpose of comparison, we consider a corresponding model with no mutation dynamics, we set δ = 0 and our model (S1) reduces to model (S2), which was used to understand the 2006 Réunion Island outbreak in 24 . We would like to see if incorporating the mutation dynamics and fitting the data would lead to any substantial differences in the estimation of the basic reproductive number.
Parameters Estimation. All parameter definitions, ranges and relevant references in model (S1) are shown in Table 1. The differences between the mutant and non-mutant strains are reflected in the following parts of the parameterization.
1. Vector competence: the intrinsic ability of a vector to support viral replication so that the virus disseminates from the midgut to the salivary glands for transmission during subsequent vector blood meals. The www.nature.com/scientificreports www.nature.com/scientificreports/ ECSA-V mutant strain can enhance the replication and the dissemination of CHIKV hence led a significant increase in CHIKV infectivity from Ae.albopictus to human during transmission 11,12,35,36 . Therefore, we assume the lower bound of β M H 2 being larger than that of β M H 1 in our sampling (S1); 2. Extrinsic incubation period (EIP): the duration for a mosquito to become infectious since it has taken a viremic blood meal. With relatively short life span of mosquitoes, it is well known that longer EIPs reduce transmission efficiency simply because fewer mosquitoes can live long enough to transmit the virus. However, experimental results indicate that the new ECSA-V variant could shorten the EIP 11,12 . Therefore, we assume γ =  www.nature.com/scientificreports www.nature.com/scientificreports/ where ⁎ W q ( ) is the real value at time q, W q ( ) is its fitting value and n is the number of data used for prediction. AIC, MAPE and RMSPE are shown in Table 4, and we can conclude that our model with mutation dynamics fits better than the simple model as (S1) has a much lower AIC value.

Calculation of Basic Reproduction Number.
Based on the estimated parameter values from the MCMC method, we calculate the corresponding uncertainty on the basic reproduction numbers from the sampled parameters. The median and confidence interval of the distribution of the basic reproduction numbers (see Fig. 8(a) . This means that based on our assumption of incorporating the linear mutation rate among mosquitoes, our data fitting results show that the CHIKV outbreak in Italy is majorly associated to the mutant strain, and there should be no outbreak of the non-mutant strain.
Moreover, we obtained that the median and confidence interval of the distribution of the basic reproduction numbers 1  for model (S2) as 1.862 (95% CI: [1.8181, 1.8858]). Figure 8(b) compares the basic reproduction numbers between the two models, 1  and  0 , and we observe a significant underestimation if the mutant dynamics is not considered.

Sensitivity Analysis and Intervention Priorities.
We perform sensitivity analysis for model (S1) by using the Latin Hypercube Sampling method 37 to generate 5,000 parameter combinations with each parameter uniformly distributed in its range in Table 1. Figure 9 shows the PRCC values 38 of these parameters, and we conclude that the basic reproduction number 0  is most sensitive to the mosquito biting and mortality rates; and is more sensitive to the mosquito transmission rate than the EIP; and is more sensitive to the transmission rate of the mutant strain than that of the non-mutant strain. 6   www.nature.com/scientificreports www.nature.com/scientificreports/ Figure 10(a,b) show that reducing either the λ M and β M H 2 pair, or the b and β M H 2 pair, will help bring  0 down below 1, hence prevent the outbreak from happening. On the other hand, as shown in Fig. 10(c), intervention strategies that only reducing the mosquito recruitment and biting rates are not efficient for eliminating the outbreak.

Conclusion and Discussion
The ECSA-V mutant strain was first discovered in the 2005-2006 chikungunya fever outbreak in La Réunion. In particular, the mutation was not detected during the first outbreak wave before September 2005, but was detected during the second wave from December 2005 which involved as many as ten times of the cases reported in the first wave 39 . It is believed that the mutant strain could have been spread by travelers later to Italy, as study confirmed the co-circulation of the non-mutant and mutant strains in the 2007 Chikungunya fever outbreak in Italy 34 . This motivates us in extending a SEIR-type model based on classic mathematical epidemical theory 40 to examine the impact of virus mutation. We obtained estimates of the basic reproduction number by using some parameters from published literature and our estimation shows that if the mutant strain was involved in the 2007 outbreak, then the disease burden could be underestimated.
The possible high infectivity possessed by the mutant CHIKV strains will not only increase the risk of infection and epidemic size, but could also initiate and sustain chikungunya fever outbreak. In our data fitting, the basic reproduction number caused by non-mutant strain is 0.698 (95%CI: [0.5213,0.8017]), but the estimated value by the mutant strain is 2.035 (95% CI: [1.9424, 2.1366]). Our simulations on varying some controllable model parameters show that the transmission rate of mutant strain from mosquito to human is an important parameter in the control of the disease spread. Some biological studies demonstrated that the adaptation in the Aedes albopictus mosquito vector due to the mutation may be associated with evolutionary success, which could further complicate the outbreak control. More seriously, ECSA-V can be seen as the initial step of adaptive mutation of the CHIKV, the second step of adaptive mutation has already been reported 36,41 . Thus new epidemiological     42 -where three out of four Aedes albopictus mosquitoes have their saliva samples detected with the mutant strain ten days after infections of the pre-epidemic strain. An estimation of the experimental mutation rate can be calculated easily: assume M t ( ) as the mosquitoes with their infection dominated by the non-mutant strain and δ as the linear mutation rate, then  The box plot of the basic reproduction numbers  1 for model (S2) and  0 for model (S1). Figure 9. The partial rank correlation coefficient (PRCC) of the basic reproduction number 0  in model (S1) with respect to some model parameters. For each parameter, the absolute value of its PRCC represents the sensitivity of the parameter -the larger the value is, the more sensitive  0 is to the corresponding parameter. * denotes the value of PRCC which is not zero significantly, where the significance level is 0.05. www.nature.com/scientificreports www.nature.com/scientificreports/

Limitations
Our model can be improved if more detailed biological evidence becomes available in the future. Here we would like to discuss some neglected factors in our model and future directions for this study. First of all, although the mutant strain have been discovered in many of the recent chikungunya outbreaks, there is limited knowledge on the mutation mechanism. Our model assumes a linear mutation rate happening among mosquitoes infected with the non-mutant strain, which is only one of the many possibilities. Thus future work should consider and model other possible mutation dynamics, fitting models with various hypothesis to real outbreak data could help select the reasonable ones. Secondly, we use real data to fit parameters with undetermined broad ranges while leaving the other parameters as fixed values including mosquito biting rate and the proportion of symptomatic individuals. However, the biting rate varies in Aedes species and could also depend on temperature and urbanization, and the proportion of symptomatic cases is always difficult to estimate on the population level. Lastly, our model structure should be adjusted based on further knowledges, such as the cross-immunity between the mutant and non-mutant strains after recovery, waning immunity 43 , and transmission ability difference between symptomatic and asymptomatic individuals.