Estimating chikungunya virus transmission parameters and vector control effectiveness highlights key factors to mitigate arboviral disease outbreaks

Background Viruses transmitted by Aedes mosquitoes have greatly expanded their geographic range in recent decades. They are considered emerging public health threats throughout the world, including Europe. Therefore, public health authorities must be prepared by quantifying the potential magnitude of virus transmission and the effectiveness of interventions. Methodology We developed a mathematical model with a vector-host structure for chikungunya virus transmission and estimated model parameters from epidemiological data of the two main autochthonous chikungunya virus transmission events that occurred in Southern France, in Montpellier (2014) and in Le Cannet-des-Maures (2017). We then performed simulations of the model using these estimates to forecast the magnitude of the foci of transmission as a function of the response delay and the moment of virus introduction. Conclusions The results of the different simulations underline the relative importance of each variable and can be useful to stakeholders when designing context-based intervention strategies. The findings emphasize the importance of, and advocate for early detection of imported cases and timely biological confirmation of autochthonous cases to ensure timely vector control measures, supporting the implementation and the maintenance of sustainable surveillance systems.


Introduction
Worldwide, arboviral diseases such as dengue, Zika and chikungunya constitute a major proportion of infectious diseases emergence [1]. In recent decades, their incidence has increased dramatically, and they have substantially extended their geographic range. While tropical countries bear the heaviest burden, some temperate areas are increasingly exposed to this threat due to the presence of Aedes albopictus, an efficient vector for their transmission. Various factors are likely to drive transmission in areas where Ae. albopictus is established and which have suitable environmental conditions [2].
In mainland France, various vector-borne transmission events have been observed over the past ten years [2,3]. Dengue fever events are the most frequent in the country [2,3], reflecting the global epidemiology of the disease [4]. However, transmission events of chikungunya virus (CHIKV) in France results in outbreaks with the greatest number of cases. This can be partially explained by the fact that Ae. Albopictus is considered as an efficient vector of CHIKV, especially East-Central-South African (ECSA) genotypes [5], whereas the species is not considered currently as a primary vector of dengue [6]. Moreover, a series of adaptive mutations of CHIKV are selected by Ae. albopictus for even better transmission, as highlighted during the major CHIKV outbreak in the Indian Ocean in 2005-2006 [7,8]. The main mutation-an alanine to valine substitution at position 226 of the E1 glycoprotein-was also present in the viral genotypes at the origin of the two main episodes of CHIKV transmission that occurred in mainland France in 2014 and 2017 [9,10]. For each of the two events, an imported case (primary case) was identified and in both situations, the imported case was returning from Cameroon [3].
This increase in CHIKV fitness in Ae. albopictus, combined with the strong capacity of the latter to invade new areas, including areas with temperate climates, explains the potential for major epidemics of CHIKV to occur throughout Europe. Accordingly, health authorities, especially in France, have proactively strengthen strategies for preparedness and response to arboviral risks [2]. To prevent local transmission a better understanding of the effectiveness of current surveillance and control strategies is needed, together with a better knowledge of the factors that influence it.
The potential for CHIKV transmission in temperate areas as well as the effectiveness of vector control remain poorly understood. The acquisition of such knowledge is essential to inform evidence-based decision-making regarding surveillance and control strategies. In this perspective, we developed a compartmental epidemiological model with a vector-host structure to estimate key model parameters of both CHIKV transmission and current vector control interventions, based on observations made during previous CHIKV transmission events. We then assessed the effectiveness of these strategies according to different virus transmission scenarios during the vector activity season.

Model
We developed a mathematical model based on a Susceptible-Exposed-Infectious-Recovered host and Susceptible-Exposed-Infectious vector framework ( Fig 1A). Vector demography was incorporated into the model in order to integrate the seasonal activity of Ae. albopictus. To this extent, the birth rate of mosquitoes was modulated to fit with the dynamics observed by entomological surveillance in Montpellier (Fig 1B), in Le Cannet-des-Maures ( Fig 1C) and with the dynamics of an artificial standard mosquito population (Fig 1D).
The rate of recruitment (birth and immigration) and the mortality rate of human population were ignored in the model because of the short time scale of each transmission event. Neither was disease-induced mortality considered given the limited size of events. Humans start in a susceptible state (S h ), and are then infected (E h ) through bites of infectious mosquitoes (I m ) at rate a, with the probability b. They become infectious (I h ) at a rate ω h , finally recover (R h ) at a rate σ, and then become immune. Similarly, susceptible mosquitoes (S m ) biting infectious humans (I h ) at rate a, can become infected (E m ) with the probability c and then infectious (I m ) at the rate ω m , the inverse of the extrinsic incubation period. The choice was made not to consider the presence of asymptomatic forms in order to have a model that was as simple as possible, while still meeting the objectives of decision support. This choice is justified by the significant effort of epidemiological investigation implemented in the area of virus circulation and was also reinforced by an analysis of the impact of the presence of asymptomatic forms on the estimates made (S1 Text).
We first used a deterministic model which is expressed by a set of ordinary differential equations. For the human compartments: For the vector compartments: The following expression was used to compute the basic reproduction number, R 0 [11]: Parameter values were taken from the literature (Table 1) with preference given to a selection of data adapted to Ae. albopictus and to the viral genotypes involved in the two transmission events [9,10].
To account for the stochasticity of both mosquito and human infections, we used, as a second step, a continuous-time stochastic version of the model, whose transitions and rates were specified by analogy with the deterministic equations (S1 Table). Numerical simulations were performed using Gillespie simulation algorithm.
Vector control strategy is mainly based on the use of adulticides with Ultra-Low Volume spraying of pyrethroids within a 250 m radius area around the residences of the cases, covering the entire area of viral circulation. Vector control was simulated by reducing the entire vector population of the area according to the effectiveness of the control measure (Eff, i.e. the proportion of the adult vector population reduced thanks to control measures) at the precise moment mosquito adulticide was sprayed. Such dynamics in vector population reduction thus reflect the immediate effect and the absence of persistence of insecticides used for urban vector control. Afterwards, the vector population returns after a few days to the level that would have been observed in the absence of control [20,21]. This dynamic covers not only the newly emerged mosquito adults, but also the recolonization by imagos present in the adjacent patches. It can be considered conservative, as it limits the effectiveness of vector control on disease incidence [21].
The model aimed to estimate two parameters: the efficacy of control measures, Eff and human host susceptibility to infection, b (i.e. the probability that a human host gets infected after being bitten by an infectious vector). To do this, the model was fitted using a Markov chain Monte Carlo (MCMC) method with the Metropolis-Hastings algorithm. The particle MCMC Metropolis-Hastings algorithm was used to evaluate the likelihood of stochastic models [22]. MCMC was run for five different chains with different initial values for both parameters in order to avoid converging towards local minima. Model fitting and simulations were performed with R [23] and, specifically the fitR package [24].

Data
Epidemiological data were collected by the French national public health agency, during two autochthonous transmission events of CHIKV (S2 Table). These events occurred in 2014 in Montpellier [9] and in 2017 in Le-Cannet-des-Maures, both on the Mediterranean coast [10]. The two events led to an equivalent number of cases (12 and 11 cases in Montpellier and Le Cannet-des-Maures, respectively). However, the transmission event started rather late in the vector activity season in Montpellier (August 30) while it started in mid-season in Le Cannetdes-Maures (July 10) [25]. In both events, the primary (imported) case was identified and the CHIKV strain belonged to the ECSA genotype with the A226V mutation of the E1 protein.
Following epidemiological investigations, the estimated total affected area of transmission was 7.1 and 8.4 ha, in Montpellier and Le-Cannet-des-Maures, respectively. Population density was estimated at 70 and 45 hab/ha, respectively. Active door-to-door investigations were conducted in both areas of virus circulation to identify symptomatic cases who had not consulted a physician [9,10]. Vector population dynamics in Le-Cannet-des-Maures were derived from the results of routine entomological surveillance conducted in 2017 in Nice (located 80 km from Le-Cannet-des-Maures) by the local public mosquito control agency (EID-Méditerranée) which used a network of 50 ovitraps. For Montpellier, we used surveillance data from a study conducted in the area in 2014 to describe mosquito population dynamics [26]. From previous studies dedicated to mosquito population density estimation [27] and nuisance perception during entomologic field investigations, we assumed a maximum mosquito population of 600 females/ha in Le Cannet-des-Maures and 730 females/ha in Montpellier. A standard vector population dynamic (for both events) was also derived from data collected in Nice by the public mosquito control agency between 2008 and 2017 in order to simulate the dynamics of vector populations usually observed in the South of France. For these standard population dynamics, we considered three different mosquito population densities with a maximum population size of 400, 800 and 1 200 females/ha in order to describe the different possible urban environments. Mosquito control measures were implemented at different time points during both local transmission events (S2 Table). Initial conditions for our model simulation were established on the basis of these data and are reported in S3 Table.

Definition of scenarios
We tested the impact of different scenarios on the magnitude of outbreaks, using the estimated parameters. First, we considered the context of the outbreaks that occurred in both areas. To do this, we initially forecasted what the size of the outbreaks in Montpellier and Le-Cannetdes-Maures would have been in the absence of any vector control action. We then compared the hypothetical sizes of each outbreak according to different timings of vector control action (i.e. scenarios with earlier or later interventions with respect to the actual timing of interventions). More specifically, we simulated response delays in vector control implementation between 0 and 90 days after primary case introduction while keeping the number of treatments and the spacing between each treatment similar to what actually occurred during the transmission events in 2014 and 2017. Finally, for each of the two outbreaks, different vector control sequences were tested by varying the number and frequency (i.e., the sequence) of control measures.
Second, with a view to generalizing the model predictions to assess the impact of response delay and time of virus introduction during a standard vector activity period, we based our model simulations on a standard mosquito population dynamic (Fig 1D). Different simulations were run for three mosquito densities corresponding to low (maximum of 400 females/ ha), medium (800 females/ha) and high density (1200 females/ha). In this context, reduced effectiveness of vector control was also considered since the value of the effectiveness of the vector control measures is not absolute. More specifically, effectiveness can be reduced because of various constraints and difficulties including limited access to the areas to be treated, absence of residents, physical barriers created by buildings, public opposition to insecticides, and unfavourable weather conditions. As part of this modelling framework, we also considered the impact of the sequences of vector control measures. In particular, we simulated different time points of primary case introduction throughout the activity season of Ae. albopictus in the South of France, from May 1 to November 30. For each outbreak, we simulated (i) different response delays between 20 and 90 days, (ii) 1 to 6 treatments, and (iii) intervals between treatments of 5, 7 and 10 days.

Results
Infectious disease transmission is mainly a stochastic process, especially in the early stage of an epidemic. Hence, particular emphasis was given thereafter to the results of the stochastic model. The deterministic model provides a good fit for the data, especially for the event in Le Cannet-des-Maures. However, we were not able to fit the stochastic version of the model with particle filtering [22], probably due to the limited number of cases. Nevertheless, the simulations of the stochastic model using the parameters estimated by the deterministic approach provided a good fit of the data for both events (Fig 2). Based on these results, R 0 was estimated at 1.86 (95% CI: 1.83-1.88) in Montpellier and 1.78 (95% CI: 1.72-1.84) in Le-Cannet-des-Maures.

Estimates of the size of the outbreaks according to different scenarios
Absence of vector control. Considering the previous estimates for both foci, we quantified the potential epidemiological impact of these two outbreaks in the absence of vector control measures by simulating 500 replicates. Accordingly, in Montpellier, the chikungunya  (Table 3 and S2 Fig). The simulations of these two different events illustrate the influence of the starting date of the outbreak on its potential size. Forecasts exhibit a large variation in the size of the outbreak due to the random nature of the stochastic equations and the high rate of extinction, as illustrated by the median values. The deterministic simulations display a higher mean number of cases with less variation in the results (Table 3 and S3 Fig).
Impact of response delay on outbreak size. The sizes of the outbreaks were simulated for different delays of vector control implementation. Response times between 0 and 90 days after primary case introduction-as described in the scenarios introduced in the Methods sectionwere considered for both events (Fig 3). For instance, the implementation of measures 10 days  Table 3). The comparison of the two events allowed us to assess both the impact of the response delay after primary case introduction and the impact of the date of primary case introduction on the size of the resulting outbreak. The delay in implementing control measures after the introduction of the primary case was longer in Montpellier (51 days later) than in Le Cannet-des-Maures (32 days). Both events reached a similar final size (12 and 11 cases in Montpellier and Le Cannet-des-Maures, respectively). However, the transmission event started at the end of the vector activity season in Montpellier (August 30) whereas it started in mid-season in Le Cannet-des-Maures (July 10) [25].   Table 3. They illustrate the epidemiological benefit of successive treatments. Unlike the time interval between treatments, the number of treatments has a strong impact. The importance of the latter parameter is related to the occurrence of the event during the period with the highest risk as illustrated by the simulated scenarios at Le-Cannet-des-Maures, while the impact of the number of treatments is not as large for scenarios simulated for Montpellier as the event there took place at the end of the vector activity period.

Estimates of the size of an outbreak for a standard vector population dynamic
Using the results outlined above for the specific two real-world events, we estimated the size of a hypothetical outbreak for a standard vector population dynamic (Fig 1D) in an attempt to generalize the previous results. We first studied the size of hypothetical transmission events as a function of the date of primary case introduction and of the response delay. We simulated three different vector densities and two values of vector control effectiveness. Our results suggest little influence of a reduction in vector control efficacy on the size of transmission events, provided that reduction was offset by an increase in the number of treatments. A contrario, vector density had a significant positive impact on outbreak size (S4 and S5 Figs). These simulations highlight the importance of rapid implementation of control measures, in particular during the peak season of vector activity, and are shown for a medium vector density with a vector efficacy of 90%, on the basis of the results of the estimates (Fig 4). Further simulations of vector control sequences (Fig 5) suggest that a single measure would have a limited impact. The optimal number of treatments in terms of vector population reduction is between 3 and 5. Our model indicates that the control of any hypothetical event of autochthonous transmission requires at least 3 treatments, or 4 or 5 if it occurs before the peak of seasonal vector activity. The impact of several successive treatments is, however, counterbalanced by the delay in response and, to a lesser extent, by the late occurrence of the virus introduction during the season. We varied the spacing between treatments to 5, 7 and 10 days, which corresponds to what is regularly implemented. A 5 days spacing is a reasonable frequency in terms of feasibility, but can sometimes be increased due to logistical or meteorological constraints. Such variation of the spacing between vector control treatments had little impact on the size of the outbreak (S6 Fig). Simulations of virus transmission under standard conditions lead to qualitatively similar conclusions and support the above results based on real-world situations. Overall, this approach allows to generalize the conclusions for different vector densities and for virus introductions throughout the vector activity period.

Discussion
In this work, we estimated transmission parameters and vector control effectiveness based on two chikungunya transmission events that occurred in mainland France using a compartmental epidemiological model. Estimated values were similar for both events with vector control effectiveness between 85 and 97%. While this result may seem quite high, it is consistent with the impact on vector population of the vector control measures actually implemented during the outbreak in Montpellier. Following the first insecticide treatment, the vector population declined drastically in the treated area, by 97% [9].

PLOS NEGLECTED TROPICAL DISEASES
A very wide range of values is proposed in the literature for human susceptibility to CHIKV infections [28]. Our results are very close to those found by some authors [11,[29][30][31]. Basic reproduction numbers were estimated at 1.86 in Montpellier (2014) and 1.78 and in Le-Cannet-des-Maures (2017). Our estimates of R 0 are the first for CHIKV transmission in mainland France. They add important results to the very fragmentary current level of knowledge about the transmission of arboviruses in temperate climates. A recent review of the basic reproduction number for dengue, Zika and chikungunya across global climate zones, proposed an average R 0 estimate of 1.88 (range: 0.46-2.94) for CHIKV in temperate areas from 9 different studies using different methodologies [32], which is consistent with our estimates.
The value of these estimates can be discussed with regard to different hypotheses, such as in particular the choice to consider the absence of asymptomatic infections. The proportion of CHIKV asymptomatic infections is usually estimated between 3 and 25% and remains highly variable across outbreaks [33]. A meta-analysis suggest that the proportion of inapparent CHIKV infections is lineage dependent and that less inapparent infections are associated with the ECSA lineage than the Asian lineage [34], which is the lineage incriminated in the two events reported here. The decision not to consider the occurrence of asymptomatic forms was mainly supported by the significant investigation effort deployed around the two transmission foci (door-to-door case search, awareness of health professionals from the area, and communication in the local media), which strongly reduces the probability of paucisymptomatic cases that would not have been identified. Moreover, a model allowing to integrate different values for the proportion of asymptomatic forms was also considered in order to evaluate the impact on the estimation of the parameters. This sensitivity analysis shows that the impact of taking into account a proportion of asymptomatic infections remains limited (S1 Text). Based on these different elements, the choice to consider only symptomatic cases was made with a view to meeting the decision-making issues discussed here, while keeping the model as simple as possible. However, this choice could be reconsidered for other objectives or even other epidemiological situations, such as large-scale epidemics.
Our study highlights the influence of factors such as the month of introduction of the virus and the delay of vector control measures. The comparison between the two events suggests that the number of cases expected in the absence of treatment is strongly influenced by the date of occurrence, which can effectively be considered as a proxy of vector population dynamics. The results based on standard dynamics allow to assess the relative importance of the various factors of interest (response time, date of start of the epidemic). The results of the model are also consistent with the two chikungunya outbreaks that occurred in Italy. They can partially explain the larger size of the Italian episodes that arose earlier in the season in Emilia-Romagna and Lazio [35,36] and for which the vector control measures were implemented approximately 2 months after the introduction of the primary case [35,37]. Under unfavourable conditions (peak vector activity and high vector densities), the potential for viral transmission remains high, even in the event of a rapid response, within 30 days. This observation justifies the need for increased surveillance during this period and the prompt implementation of preventive measures to reduce vector populations immediately after introduction of an infected case, to prevent autochthonous transmission. However, these results cannot be considered to have a 100% predictive value, since a climatic anomaly, such as heavy rainfall, is always possible and can have an impact on vector populations, as was observed in Montpellier in 2014 [26]. Our results on vector control sequences must be considered with caution, in particular those related to the spacing between consecutive treatments. Indeed, in this case, the dynamics of the populations of vectors depend strongly on the assumption we made regarding the recolonization of the treated areas. This fact underlines the need for a better understanding of recolonization after treatment. From an operational point of view, this dearth of knowledge justifies daily entomological monitoring of adult mosquito populations during all events of autochthonous transmission, not only to improve knowledge of these particular situations, but also to trigger vector control measures if the vector population increases.
Further developments in our model may be considered for the future. These include, in particular, the inclusion of climatic and environmental data for the real-time estimation of vector populations [38], and taking into account the age structure of the vector populations, another critical driver of transmission dynamics [39][40][41]. Scaling-up of the model is also important. Our model was developed based on two real-world transmission events which occurred in single patches of less than 10 ha, and assumed homogenous mixing within each patch. The definition of a spatially explicit metapopulation model would enable the spread of the virus at larger scales [42].
Quantitative results cannot be extrapolated to other viruses transmitted by Ae. albopictus (dengue, Zika). Human and vector susceptibilities to infections, extrinsic and intrinsic incubation periods, and the duration and level of viraemia are key parameters of virus transmission and may differ substantially from one virus to another. Recent events of dengue virus circulation suggest that this vector system is less efficient than that involving CHIKV [3]. However, we can reasonably assume similar conclusions in qualitative terms. Caution is likewise required in extending these conclusions to all CHIKV genotypes, as vector competence for a specific pathogen is governed by complex interactions between vector population, virus genotype and environmental conditions [43]. Experimental studies suggest a higher risk for ECSA strains to emerge in Europe, compared to Asian CHIKV strains [44]. These ECSA strains include CHIKV genotypes harbouring an alanine at position 226 of the E1 glycoprotein (i.e. without the E1-A226V mutation known for increasing virus replication and transmission by Ae. albopictus) as highlighted by the 2017 outbreak in Italy [45].
This modelling-based work emphasizes the importance of, and advocates for early detection of imported cases and timely biological confirmation of autochthonous cases to ensure timely vector control measures. In addition, given that vector density is a critical parameter for transmission after introduction, the preventive reduction of the size of vector populations is of outmost importance. To that end, awareness and mobilization of all involved stakeholders remain key elements in the control of vector-borne diseases: travellers, patients, doctors and laboratory workers for early detection and reporting, public health authorities and vector control agencies for the timely implementation of control measures, and the general public for routine control of breeding sites. Finally, our results also justify the current strategy of implementing vector control measures around imported cases, in a preventive manner, before any autochthonous transmission.
Supporting information S1 Text. Parameter inference for Chikungunya virus outbreaks for different rates of asymptomatic infections. (DOCX) S1 Table. Transitions and rates of the stochastic model.