The Role of Climatic and Density Dependent Factors in Shaping Mosquito Population Dynamics: The Case of Culex pipiens in Northwestern Italy

Culex pipiens mosquito is a species widely spread across Europe and represents a competent vector for many arboviruses such as West Nile virus (WNV), which has been recently circulating in many European countries, causing hundreds of human cases. In order to identify the main determinants of the high heterogeneity in Cx. pipiens abundance observed in Piedmont region (Northwestern Italy) among different seasons, we developed a density-dependent stochastic model that takes explicitly into account the role played by temperature, which affects both developmental and mortality rates of different life stages. The model was calibrated with a Markov chain Monte Carlo approach exploring the likelihood of recorded capture data gathered in the study area from 2000 to 2011; in this way, we disentangled the role played by different seasonal eco-climatic factors in shaping the vector abundance. Illustrative simulations have been performed to forecast likely changes if temperature or density–dependent inputs would change. Our analysis suggests that inter-seasonal differences in the mosquito dynamics are largely driven by different temporal patterns of temperature and seasonal-specific larval carrying capacities. Specifically, high temperatures during early spring hasten the onset of the breeding season and increase population abundance in that period, while, high temperatures during the summer can decrease population size by increasing adult mortality. Higher densities of adult mosquitoes are associated with higher larval carrying capacities, which are positively correlated with spring precipitations. Finally, an increase in larval carrying capacity is expected to proportionally increase adult mosquito abundance.

The lengths of the developmental periods associated to different mosquito life stages (i.e. for any ∈ = { 1 , 2 , 3 , 4 , }) at different temperatures were calibrated according the following procedure. Given the length of developmental period for temperatures ∈ = {7°C, 10°C, 15°C, 20°C, 25°C, 30°C, 33°C} as observed in [1], we assume that ( ) = ( ; ) + where ( ; ) is a parametric function of the temperature ( indicates the set of free model parameters defining ) in a suitable set of functions, comprising exponential, parabolic and logistic functions, and is a random sample of a 0 mean normal distribution with unknown variance 2 . For each considered life stage, we calibrate the function ( ; ) by minimizing the square error between predicted and observed length of the period which is defined as = ∑ ( ( ) − ( ; )) 2 ∈ . Uncertainty of estimated parameters (i.e., ) was computed following a bootstrap procedure similar to the one adopted in [2,3]. In particular, we simulated 100 different { ( )} ∈ by adding an error sampled from a normal distributed (0, 2 ) to the best interpolation ( ;̃) where the variance 2 was taken as the average of the estimated residuals associated to the best interpolation of the model i.e., the average of the quadratic differences ( ( ) − ( ;̃)) 2 . Finally, for each simulated { ( )} ∈ we repeated the optimization procedure described above. Obtained estimates of ( ; ) were used to compute the rate of development as ( ; ) = 1/ ( ; ). The same technique was applied to estimate the probability p for a fully developed pupa to become a diapausing adult as function of the daylight duration using the data presented in [4] and to estimate the mortality rates of all immature stages as functions of temperature.
However, since mortality data for different mosquito life-stages were available as survival probabilities, an additional step was required to estimate the associated mortality rates. The survival probability for different mosquito life-stages observed in [1,5] was obtained in lab conditions by following a cohort of n individuals until all of them would either die or develop in the subsequent life stage.
In our analysis, we estimate the mortality rates associated to eggs, different larval instars and pupae by maximizing the likelihood of observing the number of surviving individuals k obtained in the lab experiments at fixed temperatures, given the initial number of individuals n and a known developmental rate for each temperature and each mosquito life-stage.
For instance, specializing model M to the case of a cohort of pupae kept at a fixed temperature ̅ and starting at time 0 with the initial value ( ̅ ), one sees that, following the approach described in [6], the probability that a pupa will eventually develop into an adult (instead of dying) is ( ̅ ) = ( ̅) ( ̅) + ( ̅) , so that the number of pupae developing into adults will be a binomial of parameters ( ̅ ) and ( ̅ ).
A simple computation shows then that the maximum likelihood estimate of ( ̅ ) is Estimates of developmental times and mortality rates at different temperatures and the diapause rate associated with different daylight durations are presented in Fig. A. Our analysis suggests that at higher temperatures, survival probabilities of all stages decrease and developmental times shorten, while longer daylight durations reduce the probability that a fully developed pupa becomes a diapausing adult. In the case of eggs and pupae, lower temperatures (i.e., below 15°C) can increase the death rate as well.

Transition probabilities
Different life-stages of mosquito population are updated from time t to time t+1 according to the following procedure. Given a specific population class c we denote the set of possible transition from the class c to other population classes as Τ where each specific transition Τ occurs at a specific rate . In our simulations at each time step t the number of individuals k leaving the class c is drawn from a binomial distribution ℬ( , ∑ Τ Δ ), where n is the number of individuals in class c at time t. The number of individuals following the different transitions ̃Τwill then be computed from a multinomial with parameters k and ̃/ ∑ Τ .
For instance, in case of eggs, at each time t eggs can either develop into larval instar at a rate or die at a rate . If E is the number of eggs at time t, first one will obtain the total number of exits, k, by drawing a number from a binomial distribution ℬ( , Δ + Δ ); then the number of new larvae will be drawn from a binomial of parameters k and /( + ) (correspondingly the number of dead eggs will be a binomial of parameters k and /( + )).

Additional results
The proposed model is able to well reproduce the annual variations and the high heterogeneity observed in mosquito population dynamics among different breeding seasons (see Fig. D).