Spatial modes for transmission of chikungunya virus during a large chikungunya outbreak in Italy: a modeling analysis

Background The spatial spread of many mosquito-borne diseases occurs by focal spread at the scale of a few hundred meters and over longer distances due to human mobility. The relative contributions of different spatial scales for transmission of chikungunya virus require definition to improve outbreak vector control recommendations. Methods We analyzed data from a large chikungunya outbreak mediated by the mosquito Aedes albopictus in the Lazio region, Italy, consisting of 414 reported human cases between June and November 2017. Using dates of symptom onset, geographic coordinates of residence, and information from epidemiological questionnaires, we reconstructed transmission chains related to that outbreak. Results Focal spread (within 1 km) accounted for 54.9% of all cases, 15.8% were transmitted at a local scale (1–15 km) and the remaining 29.3% were exported from the main areas of chikungunya circulation in Lazio to longer distances such as Rome and other geographical areas. Seventy percent of focal infections (corresponding to 38% of the total 414 cases) were transmitted within a distance of 200 m (the buffer distance adopted by the national guidelines for insecticide spraying). Two main epidemic clusters were identified, with a radius expanding at a rate of 300–600 m per month. The majority of exported cases resulted in either sporadic or no further transmission in the region. Conclusions Evidence suggest that human mobility contributes to seeding a relevant number of secondary cases and new foci of transmission over several kilometers. Reactive vector control based on current guidelines might allow a significant number of secondary clusters in untreated areas, especially if the outbreak is not detected early. Existing policies and guidelines for control during outbreaks should recommend the prioritization of preventive measures in neighboring territories with known mobility flows to the main areas of transmission.


Additional information on data
Human population density data at a resolution of 250m were obtained from the Global Human Settlement database [1] and administrative unit shapefiles were downloaded from the National Institute of Statistics [2]. The categorical location of cases was assigned on the basis of the place of residence with respect to the municipality shapefiles for Anzio and Rome; cases whose residence was outside the borders of either municipality were classified as "other". Table S1 summarizes the information from epidemiological questionnaires administered during the outbreak. These data were used for the definition of cases exported from Anzio, i.e. those with an epidemiological link to this town in the days before symptom onset. Exported cases were excluded from the transmission chain reconstruction in the main analysis.

Baseline model (main analysis)
The model assumes the following expression for the force of infection to which any case j was exposed at time t: ! ( ) = & ' + , !" ; 01 ( − " ; , ) where ! ( ) is the set of notified chikungunya cases resident in the same location of j who have been infected before time t. and are two scaling factors for the transmission rates relative to the local and focal mode of infection, respectively; , !" ; 0 is a spatial kernel with parameter , regulating the distancedependent probability of focal transmission; !" is the distance between the residences of individuals j and i; " is the time at which individual i acquired infection; and ( ; , ) is the gamma-distributed (with shape and rate ) probability distribution function of the generation time, i.e. the time elapsed between infection times in a primary host and in its secondary cases. Since only information on the dates of symptom onset are available, the set of infection times = { ! } was unknown and was considered as a set of nuisance parameters (augmented data) to be estimated together with model parameters. The generation time is assumed to be to be gamma-distributed with shape a and rate b; it reflects a time-varying level of infectiousness of cases after their initial infection and accounts for the length of incubation periods in both humans and mosquitoes (intrinsic and extrinsic incubation periods), the duration of human infectiousness and the lifespan of mosquitoes. The spatial kernel accounts for focal transmission and is assumed to follow a negative exponential distribution with parameter , , !" ; 0 = ) * +*, !" [3].

Baseline model without epidemiological information
The model remains identical as in the main analysis, except for the fact the potential focal infector is searched among the set of all individuals in the dataset who have been infected before time t, N( ), independently of the location of j (N( ) substitutes N 1 ( ) in the expression of the force of infection). In addition, there are no cases classified as exported, therefore all cases for which a focal infector cannot be identified are considered as local, and the model is computed over all cases (n=412).

Power-law kernel
Here we considered a different functional form for the spatial kernel, adopting a power-law instead of a negative exponential one: , with > 1. The kernel so defined is normalized in such a way that its integral over all dij is equal to 1.

Demographic-specific susceptibility
In this model, we allowed parameter β to vary depending on the demographics of the infected individual, i.e. whether it was a child, c, defined as an individual of any gender with age strictly lower than 16 years (32 cases in the dataset), an adult male, m (168 cases), or an adult female, f (212 cases). The force of infection becomes: where δ 1 ∈ { , , }.

Seasonal transmissibility
In order to account for seasonal variations in the abundance of mosquitoes in the neighborhood of cases, and therefore on the probability of focal transmission, we assumed that transmissibility for focal infections may be dependent on time: we assumed a truncated sinusoidal function for d ( ): With t0 corresponding with May 15 th , so that focal transmission is assumed to occur between mid-May and mid-November, with a peak probability on mid-August.

Step-change in transmissibility
As an alternative way to account for seasonal variations in transmission, we included a model where a step change in the probability of focal transmission is allowed [4]: Where , and t0 are free model parameters.

Calibration
The calibration for all models was performed by adopting a Monte Carlo Markov Chain (MCMC) procedure. For free model parameters ( , , η, α, and one or more values for β, depending on model formulation), we considered recursive normal jumps and uninformative (uniform) priors. For nuisance parameters, i.e. the unknown times of infection ! , we initialized them by subtracting from the date of symptom onset a number of days sampled from the distribution of the intrinsic incubation period [5], i.e. a uniform distribution with a minimum of 2 and a maximum of 4. At each step of the MCMC, the times of infection ! for 10 randomly chosen individuals was updated, sampling from the same distribution. Parameter sets, together with the set of infection times ! and the ensuing imputed transmission chain, were accepted or rejected on the basis of likelihood values according to a standard Metropolis-Hastings algorithm. We run each model for 2.5 million iterations and considered the last 500,000 MCMC samples for analysis.
Average parameter values and their 95% credible interval computed from the corresponding quantiles of the posterior distributions estimated by the MCMC procedure are reported in Table S2 for all considered models.  To ensure convergence of the MCMC for the main analysis, we have run four parallel chains with different starting points and computed the Gelman-Rubin diagnostic Rc [6] for each parameter, finding the following values: for a: 1.45; for b: 1.43; for η: 1.23; for β: 1.11; for α: 1.18. The only diagnostics that were significantly non-compliant with the convergence criterion of Rc < 1.2 [6] are those for parameters a and b, determining the shape and scale of the Gamma-distributed generation time; this is likely due to slight variations in the accepted augmented infection times, Ej, that can result in corresponding accommodations of the curve of the generation time distribution and therefore in some correlation among the accepted parameter values. Nonetheless, posterior distributions for a and b and for the corresponding distribution of the generation time were substantially robust in the four MCMCs from the point of view of epidemiological interpretation (see Table S3).
As a further diagnostics we computed the effective sample size of the MCMC using the equation [7] = 6 )27 ∑ 9(-) ' where n=500,000 is the number of iterations and ( ) is the correlation among likelihood values at lag k. We obtained a value of ESS=175.0 (ESS>100 is a generally accepted requirement for convergence of the MCMC).

Branching process
The branching process model attempts to reproduce the transmission of chikungunya within the municipality of Anzio using estimates for the distribution of notified secondary cases, the generation time and the spatial kernel obtained via the transmission chain reconstruction model. Each simulated outbreak is initialized with a single index case, which is marked as not notified, as in the actual outbreak. For each symptomatic case, a number s of symptomatic secondary cases that will be notified is first assigned by sampling from a distribution . In particular, we assume to be a Poisson distribution with rate S, which varies over time and is set equal to the weekly estimated average number of secondary cases in the transmission chain reconstruction model ( Figure S1). Secondary cases are iteratively generated with the following characteristics: • the date of symptom onset is assigned as 6:; = "6,:< + , where is sampled from ( ; , ); • the position of the new case is assigned as ( 6:; , 6:; ) = ( "6,:< , "6,:< ) + ⋅ (cos( ) , sin( )) , where is sampled from the reconstructed distribution of transmission distances ( Figure S2), including both local and focal cases, and is uniformly sampled in [0, 2 ); • each individual may be assigned an asymptomatic status based on a Bernoulli sample with probability ps [8]; • symptomatic individuals may be notified based on a Bernoulli sample with probability equal to a reporting rate pd, while asymptomatic individuals will not be detected; For each case, secondary cases are generated until the number of symptomatic notified cases exceeds s (the case in excess is discarded). Outbreaks were halted when they reached the same number of symptomatic notified cases as the actual outbreak, i.e. 190; outbreaks resulting in stochastic extinctions in the first weeks were discarded and re-simulated.
All parameters and probability distributions for the branching process are either obtained from literature or estimated from the transmission chain reconstruction model (Table S4). We assume a baseline value for the reporting rate of 100% in the main analysis and show the model behavior when its value is reduced to up to 10% of all incident cases. In addition, we show the model performance when using different values of the notification rate for symptomatic cases before and after outbreak discovery on September 7. To do so, we switch the value of the notification rate when the number of notified symptomatic cases in the branching process reaches the value of 93, i.e. the number recorded in the Anzio outbreak on September 7. Results are reported in section 3.

Spatial mechanistic model
We developed a spatial mechanistic model of chikungunya transmission and insecticide spraying interventions to evaluate, on one hand, the effect of different intervention protocols, and on the other hand, the effect of long-distance mobility on the overall effectiveness of the currently adopted protocol.
The model represents a spatial square of side approximately 20km, divided into a 251x251 grid of cells of size . Each cell was initialized with a total number of human inhabitants and mosquitoes sampled from a Poisson distribution. For humans, the Poisson rate was given by = = = 7 , where = is the average population density of the residential area of Anzio. Similarly, for mosquitoes we consider a Poisson rate given by > = > 7 , where the density of mosquitoes > is computed by inverting the classical equation for R0 [9]: We chose R0 as the average number of secondary transmissions in absence of interventions, as estimated by the transmission chain reconstruction model; all other epidemiological parameters (see Table S5 where = and = are the respectively the number of infectious and total individuals in a given cell.
Similarly, we assumed that human infections produced by mosquitoes belonging to a given cell depend on the number of humans residing in neighboring cells with a force of infection To account for human mobility, the number of resulting new infections is distributed over the grid according to the estimated distribution of transmission distances in Anzio for focal and local cases ( Figure S2), choosing at random the direction of displacement.
Insecticide spraying interventions are modeled as follows: after outbreak detection (occurring at time in the simulation), any new infectious case triggered treatment in a square of size 2 centered on the residence of the new case, with a delay after symptom onset that was uniformly sampled over the interval [1, 2 ] days. Treatment resulted in the sudden death without replacement of a proportion of all susceptible, exposed and infectious mosquitoes in treated cells; the mosquito population recovered to pre-treatment equilibrium population sizes after a fixed number of days Θ (Table S5). , and are considered parameters of the intervention protocols and are changed one at a time in the different scenarios shown in the main analysis ( Figure 6). Duration of incubation period in mosquitos, 1/ 2.5 days [9] Duration of incubation period in humans, 1/ ! 3 days [5] Mosquito biting rate, 0.09 bites/mosquito/day [11] Probability of human to mosquito transmission per bite, # 77 % [9] Probability of mosquito to human transmission per bite, ! 70 % [9] Effectiveness of insecticide spraying in reducing the mosquito population, 86 % [12] Time to recovery of the mosquito population, Θ 14 days [12] 2. Sensitivity analyses with respect to different model formulations

Baseline model
Considering a random subsample of 10,000 reconstructed chains in the baseline model (out of a total of 500,000) rather than only the most likely presented in the main text, we found a very limited variability in the identified infectors. In particular, 60% of all non-exported cases were consistently assigned the same infector in all transmission chains; in addition, a further 36% of them were assigned the same infector in more than three quarters of transmission chains; only for the remaining 4% the candidate infector was more uncertainly identified. The proportion of focally acquired cases varied from 54.8% to 55.1%. Therefore, the conclusions of this study are robust with respect to variability in reconstructed transmission chains.   Figure S3 shows the histogram of serial intervals for reconstructed focal transmissions (maximum likelihood estimate) and the corresponding probability distribution function of the generation time estimated by the MCMC procedure (average and 95%CI computed from the joint posterior distribution of a and b).
In Table S6, we report a comparison of the main quantitative estimates across the different models and reconstructed transmission chain. In the remainder of this section we show results computed using the maximum-likelihood transmission chains for each model.

Baseline model without epidemiological information
The model still identifies two main clusters in Anzio, A1 and A2, with 78 and 57 focally transmitted cases overall, and one in Rome (with 20 secondary cases); the average speed of diffusion remained in the same range (300-600 m/month; Figure S4). The transmission distance was slightly longer compared to the main analysis, with a median of about 200m (instead of 127m), and this resulted in a slightly larger proportion of cases (68.9%, compared to 54.9% in the main analysis) being classified as focal ( Figure S5). This increase in the proportion of focal cases also resulted from the identification of plausible focal infectors for cases which were exported from Anzio. The distribution of distances for local transmission in this model shows a clear bimodality, with a first peak between 1 and 20km and a second peak, substituting cases which were classified as exported in the main analysis, between 20 and 70km.

Power-law kernel
The model with a power-law kernel does not deviate significantly from the negative-exponential kernel shown in the main analysis: parameters estimated in the two cases lead to practically identical probability distribution functions (see Figure S6). All results from the ensuing reconstruction of the transmission chains are also quite consistent with the main analyses, except for the fact that cluster A2 is split into two separate clusters (A2 and A3, Figure S7). The speed of the cluster diffusion is, again, between 250 and 400 m/month. The proportion of focally transmitted cases is estimated at about 59% and, of these, about 75% occur within 200m ( Figure S8). Figure S6. Comparison between estimated negative exponential and power law kernels.

Demographic-specific susceptibility
Conclusions from the main analysis are also robust when considering a differential susceptibility by demographic class (adult males and females, children). The main clusters and the relative size were still the same ( Figure S9) with a diffusion speed of 350-650m/month; the proportion of focal cases (68.1%) and the transmission distances (median: 175m) were similar to the model with uniform susceptibility across demographic classes shown in Section 2.2 ( Figure S10). Significant differences in the estimated values of susceptibility for different classes (Table S2) did not result in a higher proportion of transmitted cases by any demographic class, just as in the main analysis ( Figure S11).

Seasonal transmissibility
Results from the model with seasonal transmissibility were also substantially similar to all other analysis. Cluster A2 is split into two simultaneous and close clusters because of the reduced transmissibility in the early part of the outbreak; diffusion speeds were 190-450 m/month ( Figure S12). 62.6% of cases were classified as focal and 70% of these were transmitted within 200m ( Figure S13).

Model with step-change in transmissibility
When considering the possibility of a step-change in transmissibility, the model estimated an 80.1% (95%CI 72.5-87.8%) relative reduction in transmission around September 18 th (95%CI Sep 16 th -Sep 22 nd ), i.e. just a few days after systematic treatment in the area had been implemented on September 14 th . It is tempting to consider this result as an estimate of treatment effectiveness; however such an interpretation would be farfetched, since mosquito population reductions due to falling temperatures in the same period of the year may have also played a significant role. In any case, the reconstructed spatiotemporal dynamics were analogous to those identified in the main analysis: the diffusion speed of the three main clusters were 180-600m/month ( Figure S14), 67.2% of transmissions were classified as focal and 60% of them were transmitted within 200m ( Figure S15).  Figure S16 shows that the final outbreak diameter was substantially unaffected by the assumption on the value of the reporting rate, and that reporting rates pD below 50% are implausible: in such case, the outbreak would be expected to evolve over a much shorter time scale compared to observations, due to widespread background transmission.   Figure S17 shows the mean squared error between the observed number of notified symptomatic cases per week and corresponding model predictions, using different combinations of the notification rate before and after outbreak detection on September 7. The figure confirms that the best model performance occurs when the notification rate of symptomatic cases is above 70% before outbreak detection and between 90 and 100% after outbreak detection. For reference, the 2007 outbreak in Emilia Romagna, which was characterized by similar timing, intervention delay and surveillance effort [13], had an overall 78% reporting rate of symptomatic cases [8].

Impact of long-distance transmission on the effectiveness of interventions
Here, we show additional results obtained with the spatial mechanistic model under the baseline protocol of intervention ( =200m, =6 days and =7 th September) and two alternative scenarios for human mobility: one where local transmission is not allowed ("focal" scenario), and one where only focal transmission within 200m is considered ("hyperfocal" scenario). Results show that about 47% (95%CI: 36-58%) and 81% (95%CI: 71-91%) less cases after outbreak detection would be expected in the focal and hyperfocal scenarios respectively ( Figure S18). This outcome depends on the fact that secondary cases residing within a buffer of 200m from their infector will have a much lower probability of onward transmission, because the treatment performed for their infector has lowered significantly the mosquito abundance in the surrounding of their own residence: thus, the higher the proportion of cases transmitted at short distances, the better the improvement in control. This results also shows that, while current protocols for outbreak control are likely appropriate in territories with low human mobility [2], the frequent spread of infection over long distances observed for the Lazio outbreak may have significantly hurdled the ability of adopted interventions to effectively curb transmission.  Figure S18. Effect of the alternative human mobility scenarios on the potential effectiveness of the current protocol of insecticide spraying for infection control.

Statistics on average numbers of secondary cases
In Figure 4 of the main analysis we show the number of secondary cases focally transmitted by each individual, disaggregated by different covariates, i.e. demographics (adult females and males more than 15 years old, children of age below or equal 15 years), age, and population density in a 500m radius of the case.
To assess the statistical significance of these relationship, we computed Generalized Linear Models assuming a Poisson distribution of the independent variable (the number of secondary cases) with respect to the explanatory variable. For the demographic disaggregation, we found non-significant relative coefficients for either adult males (p=value: 0.164) or children (p-value: 0.340) with respect to adult females, supporting the conclusion of no difference in transmission among groups. Similarly, we found a p-value of 0.72 when age was considered as a predictor of the number of transmitted secondary cases. Conversely, the population density had a significant effect on transmission, with a p-value of 6.8 10 -7 .