Effectiveness and economic assessment of routine larviciding for prevention of chikungunya and dengue in temperate urban settings in Europe

In the last decades, several European countries where arboviral infections are not endemic have faced outbreaks of diseases such as chikungunya and dengue, initially introduced by infectious travellers from tropical endemic areas and then spread locally via mosquito bites. To keep in check the epidemiological risk, interventions targeted to control vector abundance can be implemented by local authorities. We assessed the epidemiological effectiveness and economic costs and benefits of routine larviciding in European towns with temperate climate, using a mathematical model of Aedes albopictus populations and viral transmission, calibrated on entomological surveillance data collected from ten municipalities in Northern Italy during 2014 and 2015.We found that routine larviciding of public catch basins can limit both the risk of autochthonous transmission and the size of potential epidemics. Ideal larvicide interventions should be timed in such a way to cover the month of July. Optimally timed larviciding can reduce locally transmitted cases of chikungunya by 20% - 33% for a single application (dengue: 18–22%) and up to 43% - 65% if treatment is repeated four times throughout the season (dengue: 31–51%). In larger municipalities (>35,000 inhabitants), the cost of comprehensive larviciding over the whole urban area overcomes potential health benefits related to preventing cases of disease, suggesting the adoption of more localized interventions. Small/medium sized towns with high mosquito abundance will likely have a positive cost-benefit balance. Involvement of private citizens in routine larviciding activities further reduces transmission risks but with disproportionate costs of intervention. International travels and the incidence of mosquito-borne diseases are increasing worldwide, exposing a growing number of European citizens to higher risks of potential outbreaks. Results from this study may support the planning and timing of interventions aimed to reduce the probability of autochthonous transmission as well as the nuisance for local populations living in temperate areas of Europe.


Additional information on input data
.1 shows the location of the 70 trapping sites distributed over 10 municipalities in the provinces of Trento and Belluno; the location of the study area with respect to the map of Italy is shown in the inset. Figure S1.2 shows the daily average temperatures measured in the 10 municipalities for the two mosquito seasons; the mean over different traps from the same municipality is reported. Between mid-June and September, temperatures in 2015 were constantly higher than corresponding ones registered in the same dates in 2014. The figure shows that temperature variability across sites is smaller than inter-year variability.

Mosquito population model
Equations for the mosquito population model are taken from [S1] and reported below. (Eq S1) E, L, P and V represent populations in the four developmental stages of mosquitoes, i.e. eggs, larvae, pupae and adult mosquitoes respectively. Fixed model parameters are the stage-specific mortality (m) and developmental rates (d) from one stage to the next; g V , whose inverse represents the gonotrophic cycle; and the number of eggs per oviposition n E ; all fixed parameters are temperature-dependent according to functions described in [S2],except for n E , which is set to 60 [S2]. Free model parameters are the capture rate α, which is different from zero only in days where traps are active (hence the dependence on time t in the equation); and the coefficients coding density-dependence for larval mortality, a sy , which vary by site s and year y, and represent a measure of the habitat suitability. Parameters were calibrated to reproduce capture data according to an MCMC procedure based on the Poisson likelihood of captures, as described fully in [S1]. For each site, year, disease, intervention scenario and coverage value, 100 sets of parameter values were sampled from the posterior distributions of the calibrated mosquito population model; to account for model stochasticity, 100 random repetitions were run for each parameter set. Table S1 reports the mean and 95% confidence interval of the posterior distribution of free parameters, while Figure S3.1 shows a comparison between observed (black dots) and model-predicted (red dots, with 95% confidence interval) captures at all sites and in the two years.  In Figure S2.2 we show model-predicted densities of adult female mosquitoes in the absence of control interventions for the two seasons and 10 study sites.

Modeling larvicides
The population model was modified to include larvicide interventions. We denote by c the intervention coverage, i.e. the proportion of all breeding sites on which larvicide treatment is actually performed. We assume that aquatic stages (eggs, larvae and pupae) are equally distributed across treated and untreated catch basins so that, for example, E treated = c E and E untreated = (1-c) E and similarly for larvae and pupae. Larvicidal treatment is assumed to instantaneously kill existing larvae, therefore the total larval population just after treatment, L(T + ), will be given by where Tand T + are, respectively, the times immediately before and immediately after initiation of the treatment intervention. For the duration of treatment, eggs hatching in treated catch basins are assumed to die without developing into larvae; furthermore, the density-dependent mortality parameters is reduced accordingly, in order to account for the decrease in the number of viable breeding sites. Therefore, for the duration of treatment only, the equation regulating the dynamics of larvae is represented by Equation S3.
For each type of intervention (public only or supplemented by the involvement of private citizens), we assessed larvicide effectiveness under two coverage values representing a realistic range. In general, coverage can be expressed as: where q is the fraction of (public or private) catch basins which are effectively treated and b is the total number of existing (public or private) breeding sites. Equivalently, b can be expressed in terms of the breeding site density per unit area. In a large-scale survey of different types of breeding sites conducted in urban areas in northern Italy [S3], it was found that 94% of the pupal population was produced within catch basins, while other types of water-filled containers (such as plant saucers, drums and buckets) contributed marginally to the abundance of adult mosquitoes. Based on these findings, we approximated the number of breeding sites in a given area with the number of catch basins; furthermore, we did not consider the effect of control interventions directed to the removal of other water-filled containers. Vector control interventions by municipalities can be designed to cover all public catch basins; however, some catch basins may be missed or treatment may be ineffective for various reasons, e.g. flushing, dilution or rapid dissolution of the larvicide product: therefore, we assumed that q pub is between 85% and 95%. Results from a pilot study on the involvement of citizens in mosquito control from San Michele all'Adige, a small town in the province of Trento, suggest a value for q priv between 45% and 55% (F. Baldacchino, personal communication). The much lower coverage of private interventions depends on several issues, including the presence of abandoned premises, difficulties in contacting reference persons, occasional denial of collaboration, and the actual compliance of citizens nominally adhering to the program. For what concerns the density of catch basins, we base our estimates on a survey conducted in San Michele all'Adige within the above-mentioned pilot study, which found b pub = 16.8 per hectare and b priv around 30.7 per hectare (F. Baldacchino, personal communication). These results are consistent with a previous survey in northern Italy [S3], which estimated b pub between 7 and 19 per hectare and b priv at about 36.3 per hectare. With the given estimates for q and b, we obtain a range for realistic coverage values of about 30% to 50% for public interventions (where q priv = 0), and 60% to 75% for interventions including private premises.
We considered 24 intervention scenarios, which differ by the number of treatments (effort level) within a mosquito season and by starting date. We considered between 1 and 4 treatments within a season, and we assume that each re-treatment is performed 30 days after the last treatment, so that the effect of larvicide is kept constant throughout the intervention. Possible starting dates were sampled at intervals of 15 days between the 1 st of May and the 1 st of September, and we considered only scenarios whose overall effectiveness end before October 1 st ( Figure S3.1). In this way, we obtain 9 scenarios with single interventions, 7 with two treatments, 5 with three treatments and 3 with four treatments, the latter covering almost the whole mosquito season. For interventions with involvement of private citizens, we assume for simplicity that treatments in private premises are perfectly synchronized with public ones.

Transmission dynamic model
The adopted transmission dynamic model is a standard SEI-SEIR model [S2], which can be mathematically expressed as Equations S5: where V S , V E and V I are the number of vectors in different infection states (susceptible, exposed and infected respectively). We assume that vectors from all infection states have the same temperaturedependent mortality rate m V and that new female adults (given by the term ½ d p P) begin their adult life as susceptible (i.e. no vertical transmission in mosquitoes). Mosquitoes may become infected through blood meals at a rate  V , become infectious after an average time called "extrinsic incubation period" and given by 1/ . Infectious mosquitoes remain so throughout the rest of their life. The force of infection on mosquitoes  V is given by Equation S6: where k is the mosquito biting rate, is the probability that a mosquito becomes infected upon a single blood meal on an infectious human, H I is the number of infectious humans and N is the total human population (N=H S +H E +H I +H R ). H S , H E , and H R represent, respectively, the remaining infection states for humans, namely susceptible, exposed and recovered. Humans may acquire infection with a rate , become infectious after an average time called "intrinsic incubation period" given by 1/ , and become lifelong immune after recovery with a rate . The force of infection on humans  H is given by Equation S7 = (Eq S7) where is the probability that a human becomes infected upon a single bite from an infectious mosquito.
In a previous study on a Chikungunya outbreak in northern Italy, the biting rate k was estimated at 0.09 bites per adult mosquito per day [S2]. Parameter values for the two diseases were also obtained from previously published modeling studies; they are reported in Table S2. (Eq S10) where T is the temperature in Celsius degrees, T k is the temperature in Kelvin degrees and R is Avogadro's constant. Resulting seasonal averages of parameter values are reported in Table S2, while actual temporal values over the two considered mosquito seasons are shown in Figure S4.1.

Material and methods for the economic analysis
The average cost per case and DALY loss per case were derived using a decision tree approach and considering the probabilities associated to both a dengue and a chikungunya case of being symptomatic and asymptomatic, severe and non-severe and the various disease outcomes (death, hospitalization and ambulatory assistance) ( Figure S5.1).

Figure S5.1. Decision trees for the classification of cases of Dengue and Chikungunya
In particular, the average cost per case and DALY loss per case were derived as in the following and considering base case parameters' values as indicated in Tables S5.1 and S5.2: In Tables S5.1    In particular, the cost of larvicide treatments on public spaces was estimated to be around 1.17 euros per catch basin (0.80-1.70), including costs for both personnel and larvicide products. The cost for doorto-door interventions was found to be 12.66 euros per premise (4.80-30.12), including both personnel costs for home visits and costs of larvicide products. To these, we added 0.95 euros per inhabitant (CI 0.46-1.43), for the organization of an education campaign made of public meetings for residents to explain the benefits of larvicide applications and facilitate the acceptability of the intervention in private premises [S14]. The net health benefit (NHB) was derived for each intervention and uncertainty around its average value was taken into consideration through simulations from the parameters distributions presented above.

= ∆ − ∆ /
where ∆ and ∆ are respectively DALY averted and incremental costs due to intervention, and the willingness to pay (WTP) k, as mentioned in the main text, was fixed at 35.000 euros. Such value can be interpreted as the amount of money the public Italian healthcare system is willing to pay for each DALY averted.
For each intervention scenario and each site, stochastic realizations of the NHB were drawn according to the distributions of the DALY loss per case and cost per case shown in Fig S5.2 and to the distribution of intervention costs shown in Fig. S5.3.

Figure S5.2. Distributions of DALY loss and cost per case for Dengue (respectively panel A and C) and for
Chikungunya (B and D) Optimal strategies are associated with the highest NHB, therefore probabilities in Fig. 3 and Fig. 4 of the main text are computed as the fraction of simulations for which each intervention has the highest net health benefit.

Questionnaire on budgeted and sustained cost for control programs against Ae. albopictus in Trentino Alto-Adige in 2013
VOLANO X X X X 0.02  Figures 2B and 3 in the main text for chikungunya). In particular, Figure S7.1 shows the expected number of secondary cases, Figure S7.2 the probability of local transmission and Figure S7.3 the outbreak size distributions.   In the case of dengue, given that local transmission is limited to sporadic events, the reduction in the outbreak size afforded by larviciding is negligible, and the benefits of the intervention derive mainly from a reduction in the transmission probability.

Effect of larviciding in both public and private catch basins
In Figures S8.1-S8.3, we report results on the epidemiological effectiveness of larviciding in both public and private catch basins (as in Figures 2b and 3 in the main text for public-only interventions).  shows the probability of optimal strategy in terms of net health benefit, comparing no intervention vs. public-only larviciding vs. public and private larviciding, assuming optimally timed scenarios with 2, 3 and 4 larvicide applications in a mosquito season (as in Figure 5 in the main text for single treatment).

Minimum average number of imported cases to observe secondary cases
Based on our estimates on the probability of outbreak p given an imported infection in a given site and year, we can define the probability q of not observing secondary cases after n arrivals of viraemic cases as: = (1 − ) . Therefore, for each site and year we can estimate the minimum number N of importations after which the likelihood of not observing a secondary case is below a given threshold Q: where the ] . [ operator represents the ceiling function (rounding to the first integer above). In table S9.1 we report the values of N for Q=10% and different infections, sites and treatment scenarios of public larviciding. Table S9.2 reports analogous values for larviciding in both public and private catch basins. To facilitate the interpretation of the tables, we colored the cells by risk levels defined by (arbitrary) thresholds on the number of required importations: the risk was considered high for N between 0 and 15, moderate between 16 and 30, low between 31 and 50 and negligible above or equal 51.