Modelling can reduce contamination from mosquito population control

Environmental contamination due to pest control in general, and mosquito control in particular, is an important issue expected to increase with climate change. We use a validated model for population dynamics of mosquitoes and historical environmental data to explore performance of larvicidal, adulticidal, and combined treatments. Results show that depending on treatment timing, larvicidal treatments can induce very good results, or have negative outcomes that increase overall mosquito population. Combined larvicidal and adulticidal treatments, however, exhibit much lesser dependence on timing, and therefore give the greatest chance of positive outcomes if environmental conditions are not known. Based on the results, we argue for adaptive mosquito management, in which weather data and forecasts are used to drive a model that identifies best intervals for insecticide use. Such an approach can have considerably better results than static, calendar-driven management and, therefore, considerably reduce environmental contamination. Adaptive management could consider larvicidal treatment because it gives good results if the timing is correct. Static management should, however, combine larvicidal and adulticidal treatments for the greatest chance of success.


Introduction
Environmental contamination from pesticides is pervasive in mosquito population control, and is only expected to intensify with climate change (World Health Organization 2021; Rocklöv and Dubrow 2020). Mosquitoes are highly effective disease vectors annually responsible for more than one million deaths (Caraballo and King 2014). The range of severe disease-bearing mosquitoes is expected to increase as climate change enables invasions by non-native species, and higher temperatures promote earlier and longer activity of native species (Couper et al. 2021;Colón-González et al. 2021) in temperate climate zones. Therefore, the importance of mosquito population control is only expected to increase. Mosquito control typically requires either destruction of mosquito-bearing habitats, or use of insecticides. Insecticides preserve a more natural ecosystem but may accumulate in the environment, thus posing a long-term risk.
Many adverse effects of insecticides have been documented. Non-target effects of pyrethrins and pyrethroids commonly used for adult mosquito control have been found in small-bodied arthropods (Boyce et al. 2007), lobsters (De Guise et al. 2005;Zulkosky et al. 2005), and some benthic fauna species (Kingsbury and Kreutzweiser 1979); adverse effects of malathion, also commonly used, were found in non-target flying insects (e.g., bees in Johansen (1972)) and soil invertebrates (Kuperman et al. 1999;Stepić et al. 2013). Although no deleterious effects to non-target organisms have been observed at expected field concentrations of commonly used larvicides (Lagadic et al. 2014;Lawler 2017), increased frequency of adulticide use in response to adult mosquito population invasions could present a risk to some ecosystems (Bonds 2012). Additionally, Chemical mosquito population control, as it is currently practiced in most of the world, is of limited effectiveness (Fernandes et al. 2018), thus incentivizing the use of alternative population control measures that are more harmful to the environment. Reducing insecticides would reduce pressures on the environment by minimizing the introduction of toxic substances, and reducing nontarget activities by increasing feasibility of more expensive, targeted control measures.
Here we propose an adaptive strategic framework designed to reduce pesticide use and cost of mosquito population control without sacrificing effectiveness. The framework is based on a matrix population model of flood mosquitoes developed by Lončarić and Hackenberger (2013) for the region of Eastern Croatia. We simulate and compare 33 mosquito population control strategies focusing on a) treatment start day, b) treatment duration, c) treatment type, and d) consecutive repetition of treatment patterns (see Fig. 2 for an example).

Mosquito matrix population model
This work implements the discrete-time, stage, and agestructured climate-dependent matrix population model of Aedes vexans described in Lončarić and Hackenberger (2013). The mosquito population is divided into three aquatic stages: eggs, larvae, pupae, and six adult stages: one nulliparous stage (non-reproducing females) and five gonotrophic cycles (reproducing females). Moreover, each stage is further divided according to the maximal duration of each stage at low mean daily temperatures (Kamura 1959). Hence, the egg, larval, pupa, and adult stages are divided into durations of 20, 26, 5, and 8 days, respectively ( Fig. 1), resulting in a projection matrix of dimensions 99 9 99. Since mosquito reproduction is not limited by the number of males (Takken et al. 2006), and only females are involved in blood-feeding and transmission of vector-borne diseases, the model focuses on females only. The matrix population model was developed for the city of Osijek, located in the continental part of Croatia 10 km southwest of a marshland area and Kopački Rit Nature Park. Flooding dynamics, i.e., the quantity of water originating primarily from the Danube river and to a smaller extent from the Drava river, determine ecological features of this marshland area. The area is an ideal floodwater mosquito breeding site, and the floodwater mosquito species A. vexans dominates in prevalence and abundance (Merdić et al. 2010). Given the biology of the species and the area's ecological features, the main drivers of mosquito population dynamics are photoperiod, temperature, rainfall, and flooding (Danube and Drava water levels). Thus, the A. vexans matrix population model has the general form where t is time in days, n is the vector containing the number of individuals in each stage and ''age''-group, TDPM denotes a three-dimensional projection matrix that is a nonlinear function of mean daily temperature (T), rainfall (R), photoperiod (Php), and the Danube (WL Da ) and Drava (WL Dr ) water level. The TDPM has dimensions 99 9 99 x number of days for which the mosquito population dynamic is simulated. Additional information on the structure of the TDPM can be found in SI1. Factor q -1 is the reciprocal of Leslie's density-dependent factor based on the number of larvae present at time t. The densitydependent factor (q) is set only for larvae and is modified to account for the variable carrying capacity arising from the flooding dynamics. The carrying capacity coefficient (K f ) is modelled using an empirically obtained exponential function dependent on the Danube's water level (WL Da ) in [cm], which serves as a proxy for the surface of the flooded areas: Overview of the modeled mosquito life-cycle. The life-cycle is divided into nine stages, each of which is further divided into days of duration according to the maximal duration of each stage at low mean daily temperatures. The egg stage (E) is divided into a 20-day duration, the larval (L) stage into a 26-day duration, the pupal (P) stage into a 5-day duration, and the nulliparous (N) stage and all five gonotrophic cycles (GCs) into 8-day durations. Depending on the environmental conditions, possible transitions between developmental stages are also indicated on the life-cycle graph. The figure is adapted from Lončarić and Hackenberger (2013), and used here with the authors' permission The carrying capacity for larvae at each time interval (t) is calculated using where K 0 is the carrying capacity when there is no flood (i.e., the Danube's water level is below the flooding threshold of 200 cm) and K f is the carrying capacity coefficient at time t.
The density-dependent factor (q) is calculated as follows: where K is the carrying capacity for larvae at time t, k is the dominant eigenvalue of the projection matrix at time t, and n LAR is the number of larvae present at time t. For more details on modelling the variable carrying capacity and its implementation in the A. vexans model, please consult Lončarić and Hackenberger (2013). Since A. vexans mosquitoes overwinter as eggs, the model includes/implies removal of all life stages except eggs with the onset of unfavorable conditions (low temperature and photoperiod) at the end of the year. Therefore, at the beginning of the following year, only eggs are present, which hatch when being flooded under favorable environmental conditions. Please see Lončarić and Hackenberger (2013) for details and discussion on model construction, parameters, performance, calibration, and validation.

Mosquito population control simulations
The effects of chemical mosquito population control interventions were simulated for an 11-year period, 2005-2015. Mean daily temperature, rainfall (measurement station: Osijek), Danube (measurement station: Batina), and Drava (measurement station: Osijek) water levels for the period were obtained from the Croatian Meteorological and Hydrological Service in Zagreb, Croatia. Mosquito population dynamics were calculated independently for each year, with the same initial number of eggs.
Effects on adult population abundance of three insecticide treatment types were investigated: larvicide (L), adulticide (A), and combined larvicide-adulticide (C). A total of 11 basic treatment strategies have been considered. Each basic treatment strategy consisted of N R repeats of a treatment pattern, which consisted of N T consecutive days of treatment days, and a pause of N P days, i.e. N P consecutive days without treatment (see Fig. 2 for an example). The treatment pattern parameters varied as follows: (i) N R was set to one, two or three; (ii) N T was set to one, three, seven, 14 or 28; and (iii) N P was set to seven or 14.
Combined with treatment types (L/A/C), this resulted in 33 treatment strategies in total (Table 1).
Each of the 33 treatment scenarios in Table 1 was simulated for all possible days of treatment start. Thus, for each treatment strategy, a total of 365-N R * (N T ? N P ) simulations were run for each year (Fig. 2), ensuring that a particular control strategy and starting date were selected from uniform distributions (i.e., all control strategy-start day combinations were equally likely). As the mosquito population control interventions were simulated for an 11-year period, the outcome of a particular treatment strategy-start day combination was represented by 11 ''samples''. During simulations, segments of the population vector corresponding to the targeted stage(s) were reduced by a fixed percentage (efficacy; eff) every day of the treatment. The efficacy of all treatment types (L, A and both treatments in C) was fixed to eff = 55%, a mean efficacy level recorded in previous mosquito population monitoring studies (e.g., Taylor and Schoof 1971).
However, using the same efficacy for combined (C) as for individual-stage (L and A) treatment types may have introduced a bias into the analysis. Retaining efficacy of 55% for both adulticidal and larvicidal treatments implies that twice as much insecticide has been used in C than in L and A treatment types; this might affect comparisons between treatment type efficiencies: C may be more efficient simply because it uses more insecticide than either L or A treatment types. To distinguish between the effects of combined treatments from the effects of using more pesticides, we also simulated C treatments with half the efficacy for each of the two components, eff = 27.5%.
In the example above ( Fig. 2), in L treatments, all adult larval classes (21-46) were reduced by 55% each day of the treatment (grey days in Fig. 2); in A treatments, all adult classes (52-99) were reduced; in C treatments, both larval (21-46) and adult classes (52-99) were reduced by 55% initially, and 27.5% for analysis of potential bias. Each day of active mosquito population control (grey squares in Fig. 2), the population vector (n) was modified according to the treatment type as follows: The outcomes of treatment strategies were compared by the relative change in the abundance of adult mosquitoes. For any given treatment, the relative change in adult abundance, R TS_d , was calculated as: where N B is the total number of mosquitoes with no mosquito control interventions, and N TS_d is the total number of adult mosquitoes for treatment scenario (TS) starting at day d.

Computations and data analysis
The discrete-time stage and age-structured climate-dependent matrix population model of A. vexans were implemented in Python programming language (Python Software Foundation), with packages NumPy (Harris et al. 2020), numba (Lam et al. 2015), and pandas (McKinney 2010) used for model construction. All data were analyzed in R (R Core Team 2021).
Probabilities of favorable mosquito population control outcomes for all treatment start days and all treatment strategies (reduction of the abundance of adult mosquitoes; P(R TP_d \ 0%)) were determined based on the respective empirical cumulative distribution functions computed using the ecdf function in R. In line with the preparation of mosquito population control simulations described in Sect. 2.2., the probabilities of favorable mosquito population control outcomes assume that the control strategies Fig. 2 Schematic representation of the implementation of mosquito population control strategies shown on the example of the treatment strategy L/A/C-3R-7 T-7P. The treatment pattern consisting of seven treatment days (N T = 7) and seven pause days (N P = 7) is repeated three times (N R = 3) to generate the general treatment strategy. The general treatment scenario is then projected onto the vector representation of the year, with a one-day shift in each iteration. Individual treatment patterns are highlighted using bolded borders. Grey squares represent treatment days, while white squares correspond to days without mosquito population control (pause) Table 1 List of mosquito population control strategies simulated in this study. L, A, and C denote respectively larvicidal, adulticidal, and combined types of treatment.
N R represents the number of consecutive repeats of treatment patterns, N T is the number of days treatment is (consecutively) implemented within a treatment pattern, and N P is the number of days without treatment within a treatment pattern (pause). N P is not defined for strategies with a single repeat of treatment patterns and has been set to zero and their respective starting dates were selected from uniform distributions (i.e., all control strategy-start day combinations were equally likely). Differences in overall effectiveness of treatment strategies implemented on various days of the year (R TP_d ) across treatment types (L/A/C) were calculated using robust one-way ANOVA (t1way function) and determined based on linear contrasts (lincon function; a = 0.05) using the WRS2 package (Mair and Wilcox 2020). Statistically significant differences from normality were detected using the Shapiro-Wilk test. Simulation results were visualized using the matplotlib (Hunter 2007) package in Python and ggplot2 (Wickham 2016) in R.
The computer code used for simulating mosquito population treatment scenarios is available on GitHub (SI2).

General overview of mosquito population control outcomes
Treatments could have both negative and positive effects on adult abundances. The relative change in the abundance of adult individuals resulting from all treatment strategies ranged from -59.39% (for treatment strategy C-3R-7 T-14P starting at day 126 in 2006.) to ? 58.79% (for treatment strategy L-1R-14 T-0P starting at day 237 in 2008.) (Fig. 3).
Although treatments had a slight generally positive outcome, chances of any particular treatment to be beneficial in all years is extremely low. The median relative change in adult mosquito abundance for all treatments simulated during the mosquito emergence and breeding season (days 60-310) was -2.57%, while the mean was -5.40%. The probabilities of favorable outcomes (P(R TP_d \ 0%)) during the mosquito emergence and breeding season ranged from 48.48% to 54.56%, with a mean of 51.07%. The probabilities of other treatment outcomes were as follows: 28.45% (range 24.39%-32.45%) for a relative change in abundance of adults of -5% or more (P(R TP_d \ = -5%)); 17.95% (range 15.65%-20.89%) for a relative change in abundance of adults -10% or more (P(R TP_d \ = -10%)); and 7.57% (range 5.68%-9.76%) for a relative change in abundance of adults of -20% or more (P(R TP_d \ = -20%)).
Ranges and median values for all treatments as functions of start day suggest that early treatments are more likely to give positive outcomes (Fig. 4). Further analyses distinguishing treatment types and strategies show additional detail. Sections 3.2. and 3.3. compare outcomes of different treatment types, and different treatment strategies, respectively. Results in Sect. 3.4. demonstrate that combined treatments (C) are superior even when efficacy is halved, i.e., that simulation results hold even after accounting for higher insecticide use in combined treatments.

Outcomes by treatment types
Combined treatments are considerably more likely to have positive outcomes. The relative change in abundance of adults due to larvicide treatments ranged from -55.82% to ? 58.79%, with a median of -1.93% and a mean of -4.16% (Fig. 5a). For the adulticide treatments, the median of relative change in abundance of adult mosquitoes was -0.48%, and the mean was -0.78%, with a range between -36.74% and ? 50.78% (Fig. 5b). The relative change in the abundance of adult individuals resulting from combined treatments ranged from -59.39% to ? 20.28%, with a median of -7.99% and a mean of -11.37% (Fig. 5c). Favorable outcomes had mean probabilities of 62.31% (SD: 21.70%), 56.64% (SD: 19.67%), and 90.52% (SD: 15.85%), respectively for all simulated larvicide, adulticide, and combined treatments.
Interestingly, a clear contrast emerges between treatments of identical total duration (and thus the equal amount of insecticide applied) but a different schedule. Namely, a consistently lower probability of a favorable outcome was observed for all possible treatment start days for 1R-3 T-0P (treatment performed at three consecutive days, once in a year) than for the treatment strategy 3R-1 T-7P (three repetitions of treatment pattern consisting of one treatment day and seven pause days). A similar result emerges when Fig. 4 Relative change in the abundance of adult individuals for all treatments as a function of the treatment start day. The black line corresponds to the median relative change in adult mosquito abundance resulting from all treatments simulated across all treatment types, strategies, and years for a given start day, while the grey ribbon denotes the ranges (minimum and maximum) of the relative change in abundance of adults resulting from all treatments simulated across all treatment types, strategies, and years for a given start day. The red dashed line represents a 0% change in adult abundance comparing probabilities of population reduction corresponding to treatment strategies 1R-14 T-0P and 2R-7 T-7P (Fig. 6).
Combined treatments (C) have a consistently higher mosquito population reduction than single larvicidal (L) or adulticidal (A) treatments for all treatment strategies combined treatment (C). All treatments use 55% efficacy. The black line corresponds to the median relative change in adult mosquito abundance resulting from all treatments simulated across all treatment strategies and years for a given treatment type and start day, while the grey ribbon denotes the ranges (minimum and maximum) of the relative change in abundance of adults resulting from all treatments simulated across all treatment strategies, and years for a given treatment type and start day. The red dashed line represents a 0% change in adult abundance Fig. 6 Relative change in abundance of adult individuals for two different mosquito population control strategies with the same total duration but a different treatment strategy: a relative change in abundance of adult individuals resulting from the application of treatment strategy consisting of 14 consecutive days of treatment with one repetition (strategy 1R-14 T-0P) across all treatment types and years for a given start day; b relative change in abundance of adult individuals resulting from the application of a treatment strategy consisting of seven treatment days and seven pause days, repeated two times (strategy 2R-7 T-7P) across all treatment types and years for a given start day. The black line corresponds to the median relative change in adult mosquito abundance resulting from the respective treatment strategy simulated across all treatment types, and years for a given treatment start day, while the grey ribbon denotes the ranges (minimum and maximum) of the relative change in abundance of adults resulting from the respective treatment strategy simulated across all treatment types, and years for a given treatment start day. The red dashed line represents a 0% change in adult abundance (Fig. 7). The probability of a favorable outcome is also generally higher (75-100%) for combined than for larvicidal and adulticidal treatments of the same treatment strategy. The probability of a favorable outcome for combined treatments is consistently high throughout the year.
Regularities in the likelihoods for larvicidal and adulticidal treatments also emerged. Larvicide treatments exhibited high probabilities of population reduction (75-100%) if started at days 80-120 and 230-270. The highest chances of a favorable outcome (60-90%) for adulticide treatments appeared for treatment start days 110-150 and 260-290.

Effects of insecticide amounts in combined treatments
Halving treatment efficacy to compensate for increased insecticide use in combined treatments does not affect results (Fig. 8). In combined treatments, both larvicides and adulticides are used, effectively doubling the amount of insecticide used if the standard efficacy of 55% applies.
To counter the doubling, we halved the efficacy to 27.5%. Qualitatively, results remain the same: combined treatments at 27.5% efficacy were better than either larvicidal or adulticidal treatments at 55% efficacy. Quantitatively, however, some scenarios did have less positive outcomes. This effect was most evident for treatment strategy 1R-14 T-0P, with treatments of efficacy level 27.5% started at days 130-150 and 200-220 having a greater probability of a favorable outcome (85-100%) compared to the same treatment strategy implemented at 55% efficacy (65-80%). Generally, an increase in the efficacy (insecticide quantity) results in greater variability of the relative change in adult abundance. Therefore, the outcome of the combined treatment becomes less predictable in some scenarios (Fig. 8). However, mean probabilities of adult population reduction were over 90% for both efficacy levels simulated. Fig. 8 Relative change in the abundance of adult individuals resulting from the application of combined treatments implemented at two different efficacy levels for treatment strategy 1R-14 T-0P: a 27.5%, b 55%. The black line corresponds to the median relative change in adult mosquito abundance resulting from the respective treatment efficacy level simulated across all years for a given treatment start day, while the grey ribbon denotes the ranges (minimum and maximum) of the relative change in abundance of adults resulting from the respective treatment efficacy level simulated across all years for a given treatment start day. The red dashed line represents a 0% change in adult abundance

Discussion
We used a validated matrix population model for mosquito population dynamics to investigate the efficiency of population control scenarios. Since the discrete mosquito population dynamics model implemented in this research considers each mosquito developmental stage and the transitions through the life cycle, it allows accurate prediction of variations in mosquito abundance due to control measures targeting different developmental stages.
The need for optimization of mosquito population control interventions in terms of timing, duration, type and treatment strategy has been recognized by the World Health Organization (World Health Organization 1982). Furthermore, there is an increasing emphasis on ''integrated vector management'' approaches and decisionmaking considering the biology and ecology of the target species, as well as environmental factors driving mosquito population dynamics across both space and time (AMCA 2021;Cai et al. 2021;Ferguson 2018;Fouet and Kamdem 2019).
Consistent with previous research, simulations clearly indicate greater effectiveness of combined larvicide-adulticide treatments than treatments targeting only one of these stages (Fig. 7). For example, White et al. (2011) demonstrate that combining interventions such as the application of larvicidal or pupacidal agents that target the aquatic stages of the mosquito life cycle with vector control interventions directed against adult mosquito stages (longlasting insecticide-treated nets or indoor residual spraying) can lead to substantial reductions in adult mosquito density. Kiware et al. (2017) developed a model of mosquito population dynamics to optimize the impact of vector control intervention combinations to suppress malaria transmission. Based on model simulations, the strategy of ''attacking mosquitoes on multiple fronts,'' i.e. simultaneous implementation of different vector control interventions, is more effective than performing a single vector control intervention in isolation. A similar conclusion was drawn by Rafikov et al. (2015) who found that combination of three vector control strategies (namely insecticide control, sterile insect technique and reduction of the environmental carrying capacity) not only increased the efficacy of vector control compared to scenarios involving individual vector control strategies, but also significantly reduced the costs of mosquito population control.
Reducing the amount of insecticide used in combined treatments did not change results. In principle, combined treatments use more insecticides than treatments targeting a single developmental stage simply because multiple chemicals are used. Therefore, a greater probability of positive outcomes than targeting single developmental Fig. 7 Relative change in abundance of adult individuals resulting from the application of larvicide a, d, adulticide b, e, and combined c, f treatments implemented across all years for treatment strategies 1R-14 T-0P (first row) and 3R-1 T-7P (second row). The black line corresponds to the median relative change in adult mosquito abundance resulting from the respective treatment strategy simulated across all years for a given treatment start day, while the grey ribbon denotes the ranges (minimum and maximum) of the relative change in abundance of adults resulting from the respective treatment strategy simulated across all years for a given treatment start day. The red dashed line represents a 0% change in adult abundance stages could simply result from the higher dosage. Our results, however, suggest this is not the case: reducing insecticide, modeled as halving the efficacy of treatments, does not change the positive outcome patterns (Fig. 8).
Counterintuitively, treatments can increase the mosquito population throughout the year. In mosquito populations, density-dependent growth occurs in the larval stage, with the development of individuals limited by available space and, in some environments, food (Jian et al. 2014;Legros et al. 2009). If this natural population regulation mechanism is disrupted at the wrong time, treatment can result in a significant increase in adult mosquitoes compared to no mosquito control interventions (Fig. 3b). As shown in Eqs. 1-4, the model used in this study implies limiting the number of larvae not only based on the number of individuals in the larval stage, but also under the influence of environmental conditions (flooding dynamics). Targeting the larval stage during a phase of accelerated development supported by favorable environmental conditions disrupts the density-dependence mechanism at its peak (when it is the strongest), thus triggering a strong compensation of destroyed larvae and a rapid increase in adult population abundance.
A similar mechanism operates when targeting the adult stage. Since the number of adult individuals directly determines the number of larvae, treatments targeting the adult stage during a strong density-dependence of the larval stage also trigger compensation in the larval stage. The compensation causes a rapid increase in adult abundance, with a delay corresponding to the time required for development from the larval to the adult stage. Combined treatments significantly reduce the likelihood of adverse outcomes driven by reduced activity of density-dependence mechanisms because combined treatments simultaneously act on the adult stage, thus reducing the population rebound.
Generally, optimal timing depends on the treatment type. Larvicidal treatments perform better early in the year, starting around day 100 (Fig. 7). This can be explained by considering the life cycle of the mosquito. The development of the first larval population of Aedes species is (given an adequate photoperiod) triggered by flooding (Porphyre et al. 2005;Shaman et al. 2002). First widespread flooding, therefore, synchronizes larval development, with most of the mosquitoes undergoing the sensitive stage simultaneously. Treatment at that particular time yields great results. Since treatments early in the season are more likely to correspond to the initial flooding correctly, such treatments also give better outcomes. Later reproductive cycles lose synchronicity, thus reducing the impact of larvicidal treatments.
On the other hand, adulticidal treatments are only effective when there is a significant number of adults. This, however, happens only later in the season and highly depends on the history of environmental drivers during the year. Therefore, only late (days 260-290) adulticide treatments consistently produce positive outcomes (Fig. 7). Additionally, there is a shift of about one month between periods favorable for larvicidal and adulticidal treatments. The shift corresponds to the time from the development of a critical number of larvae to the emergence of a critical number of adults.
Combined treatments offer the best chance of getting a positive outcome, and yield the best treatment result (Fig. 3a). At the same time, a larvicidal treatment gave the worst result (Fig. 3b). As discussed above, this is due to the interplay between treatment timing and density dependence: if timed correctly and early in the season (during the first flooding, when emergence of the first larval population is expected), an appropriate treatment can almost completely suppress the population. If timed badly, the treatment negates larval density dependence, and the population can actually grow to larger than without any treatment. Larvicidal treatments result in a wide range of outcomes. Combined treatments, on the contrary, tend to have positive outcomes because they effectively bet-hedge by treating both larvae and adults, thus avoiding negation of density dependence.
Shorter, pulsed treatments offer a distinct advantage over longer treatments of the same total duration (Fig. 6). Even though the total amount of insecticide introduced is the same, pulsed treatments can affect the (sub)population during its recovery from the previous pulse or environmental disturbance, thus providing a better overall outcome even though the reduction in mosquito population during the first pulse in pulsed treatment may be much smaller than the reduction caused by a prolonged treatment. This finding is supported by Cailly et al. (2012), who used a continuous, mechanistic, climate-driven model of seasonal mosquito population dynamics to compare the efficiency of mosquito control strategies targeting larvae. They found larvicide spraying at regular time intervals is more efficient than spraying only when the abundance of host-seeking females reaches a given threshold.
Our results also indicate that the outcome of mosquito population control is not necessarily proportional to the amount of pesticide applied. The implementation of mosquito population control is based on applying a constant amount of larvicide and adulticide agents, most often expressed in units of volume or mass per area (e.g., mL ha -1 or g ha -1 ) (Boubidi et al. 2016). Although such an approach is not practiced (as far as the authors know), there is a need to optimize the amount of pesticide applied according to the type and timing of treatment, thus reducing environmental contamination and mosquito control costs.
We have shown that timing is crucial, particularly timing relative to important environmental cues such as flooding. However, these general results are becoming less relevant due to climate change. In the past, weather patterns, and therefore these cues, were somewhat predictable. Climate change, however, is rapidly affecting weather patterns, reducing the relevance of past experiences. Unless easy-to-use tools for optimizing insecticide treatments are developed, the efficacy of insecticide use will decline significantly (Rocklöv and Dubrow 2020;Colón-González et al. 2021). The decline will either prompt an increase in insecticide use and the related environmental toxicity and risk, or cause an increase in mosquito populations, thus negatively affecting human wellbeing.
There have been significant advances in the use of available datasets and targeted data collection to predict mosquito populations, particularly in preventing the spread of mosquito-borne diseases. Joshi and Miller (2021) give a comprehensive overview of state-of-the-art. However, there is a lack of software solutions that would serve as auxiliary tools and a basis for decision-making for local authorities to optimize mosquito population control and change its nature from reactive to preventive. We, therefore, call for a significant additional effort in spatially explicit modelling of mosquito populations and insecticidal treatments, which should rely on weather forecasts and flood modelling to inform adaptive mosquito population control measures.

Conclusions
The effectiveness of mosquito population control depends on the time of onset of treatment, its duration, type, and treatment strategy. When designing static mosquito management, i.e., one driven by the calendar rather than actual environmental conditions, combined treatments are most likely to give best results. If timing of treatment is right, however, larval treatment could yield superior results. Based on real-time environmental data and short-term weather forecasts, models can correctly predict mosquito development and population dynamics , thus identifying the best period for larvicidal treatment a few days ahead. Switching to an adaptive management framework based on such model predictions can therefore considerably improve outcomes. Static mosquito control strategies should therefore use combined treatments; model-driven adaptive strategies can, however, also use larval treatments for considerable improvement in outcomes and, therefore, reduced environmental contamination. We hypothesize that other pest management strategies can benefit from the modelling approach as well.

Conflict of interests
The authors have no relevant financial or nonfinancial interests to disclose.