Modelling the effect of temperature on the seasonal population dynamics of temperate mosquitoes

(cid:1) We model temperate mosquito species using delay-differential equations. (cid:1) Variable time delays incorporate effects of temperature on life stage duration. (cid:1) The effects of changing climate on mosquito seasonality are explored. (cid:1) Effects of seasonal temperature on abundance can be complex. (cid:1) An increase in mosquito abundance is expected under predicted UK warming.


Introduction
Climate change is expected to affect the distribution and seasonal dynamics of mosquito populations, with substantial implications for disease seasonality and persistence (Githeko et al., 2000). Globally, mosquito-borne diseases are a major public health concern (Gubler, 2002), and increasingly so in temperate climates, with disease caused by mosquito-borne pathogens including West Nile (Sambri et al., 2013), Chikungunya (Grandadam et al., 2011), Usutu (Weissenböck et al., 2002) and Toscana (Charrel et al., 2005) viruses being reported in Europe in recent years. A range of factors including societal, land use and habitat changes, have been linked to changes in exposure to mosquitoes, which vector these diseases (Gubler, 2002). These factors combine with widespread predictions that rising temperatures may increase mosquito population size, development rates and per host biting rate, producing increases in the incidence of mosquito-borne diseases (Mirski et al., 2012). The non-linear and opposing impacts of temperature on vector demographic rates require detailed modelling, incorporating the various biological processes at play, if we are to understand the likely impacts of climate change on vector abundance and whether predictions of increased exposure to vectors will be borne out (Rogers and Randolph, 2006).
To model disease seasonality and persistence it is essential that epidemiological models are coupled with accurate seasonal predictions of vector density (Lord, 2004). This is likely to become increasingly important in coming years, as the climate is expected to become more variable (Jones et al., 2009). Without accurate predictions of vector density, calculations of the basic reproduction number of the disease (R 0 ), which describes ability of a disease to persist in a susceptible population, will be subject to considerable error. However, many recent epidemiological studies of vector-borne disease do not account for seasonality in vector populations (Cruz-Pacheco et al., 2005;Bowman et al., 2005;Wonham et al., 2004). In this study we utilise a delay-differential equation (DDE) framework to explicitly model the effects of temperature on each of the life stages of the vector population. This gives valuable insight into how inter-annual variability in temperature may affect mosquito populations, which will have a direct effect on disease transmission.
There is considerable evidence showing that environmental drivers have a large impact on both the mosquito life and disease transmission cycles (Jian et al., 2014;Lebl et al., 2013;de Almeida Costa et al., 2006). Hence, understanding these mechanisms and resulting effects will aid our predictions for vector-borne diseases. Increasing temperature is known to decrease the length of time spent in each of the immature stages and to decrease the lifespan of adults Loetti et al., 2011). Furthermore, the death rate of the immature stages is strongly linked to temperature (Ciota et al., 2014;Madder et al., 1983;Loetti et al., 2011) and increasing temperature leads to decreases in the lengths of both the gonotrophic cycle (the time required for location of a blood meal, embryonic development and egg-laying) Vinogradova, 2000) and the extrinsic incubation period (time between an adult vector contracting a pathogen and becoming infectious, referred to as EIP) (Hartley et al., 2012). With photoperiod, temperature is also believed to control the induction and termination of diapause (Spielman and Wong, 1973;Madder et al., 1983). Overwintering behaviour is a vital aspect of disease transmission because diapausing mosquitoes may act as a pathogen reservoir between seasons when hosts are no longer infectious (Nasci et al., 2001). The winter survival of mosquitoes may therefore influence disease persistence and seasonal population dynamics. Diapause, which we consider in this study, has been ignored in previous DDE models as the focus has generally been tropical mosquito populations which do not exhibit this behaviour (White et al., 2010;Beck-Johnson et al., 2013). Clearly temperature exhibits complex and opposing effects on different parts of the mosquito life and transmission cycles. As such, mathematical models are needed which explicitly incorporate the impact of temperature on each life stage if we wish to understand its effect on seasonal abundance.
Stage-structured matrix population models allow populations to be broken down according to their life stages, with climate dependencies on each stage included, making them a popular tool for modelling mosquito populations (Lončarić and Hackenberger, 2013;Schaeffer et al., 2008). However, in many cases, these models make the assumption that development time for each stage is fixed (Schaeffer et al., 2008), when temperate species' generation times at high temperatures, can in fact exceed those observed at lower temperatures by a factor of more than four (Loetti et al., 2011). Combined with the effects of temperature on larval survival (Loetti et al., 2011), these models may substantially under-or overestimate population size. One solution is to split each stage into multiple sub-stages and allow temperature to influence transition probabilities (Lončarić and Hackenberger, 2013). This approach substantially improves upon the assumption of a fixed development time, however the transition probability remains dependent only on the current temperature and life stage, without accounting for temperatures experienced earlier in development. Similar problems are faced by models based on ordinary differential equations (ODEs) (Erickson et al., 2010;Alonso et al., 2010).
In contrast to discrete-time matrix-based approaches, and to allow the duration of each life stage to vary continuously, we can look to work carried out by Gurney et al. (1983), who advocated using a system of DDEs to model stage-structured populations. These DDEs are derived from continuous age-structured partial differential equations (PDEs) by assuming lumped age classes, which are appropriate for many insect species with distinct life stages. This formulation can then be extended as shown in Nisbet and Gurney (1983) to allow variation in stage duration dependent on biotic or abiotic factors.
The stage-structured DDE framework has recently been utilised to investigate ectotherm population viability under a climate warming scenario (Amarasekare and Coutinho, 2013) and to model the life cycle of Anopheles, a genus of tropical mosquito species that vector malaria (Beck-Johnson et al., 2013). To our knowledge, theirs is the first DDE model, applied to mosquitoes, where the total length of the immature life stages was assumed to be temperature-dependent. They used this model to make inferences about the ability of Anopheles to transmit malaria under various constant temperature scenarios. As previously alluded to, further complications arise for temperate mosquito species due to greater fluctuation in seasonal temperatures and overwintering behaviour. To simplify their model, Beck-Johnson et al. (2013) also assume that all immature life stages are affected by temperature in the same way, which allows the delay equations to be transformed onto the physiological timescale and reduced from three separate equations into one. Our review of the literature suggests that many temperate mosquito species have differential, temperature-driven rates of development and survival between life stages, so we allow each stage to have a unique relationship between temperature and vital rates.
The DDE modelling framework can be applied to any temperate mosquito species by choosing parameters and vital rate curves to fit the data available for that species. We focus our attention on modelling Culex pipiens, the most abundant potential vector of West Nile Virus (WNV) in the UK. WNV is the most significant cause of mosquito-borne disease in temperate regions including Europe and North America (Pervanidou et al., 2014;Barzon et al., 2013;Petersen et al., 2013) and substantial uncertainty remains about its seasonal dynamics (Golding, 2013).
Temperature is often considered to be the primary driver of mosquito development (Knies and Kingsolver, 2010). By developing a stage-structured model which explicitly captures its effects on each life stage we aim to understand how predicted seasonal temperature changes may affect mosquito seasonality. Whilst mosquito surveillance in the UK has increased in recent years (Townroe and Callaghan, 2014;Medlock and Vaux, 2015;Townroe and Callaghan, 2015) there remains a lack of publicly available empirical field data for UK Cx. pipiens populations, so accurate predictive modelling is needed to understand seasonal population dynamics and how vectors may be impacted by changing climate. We begin to address this issue here by providing testable theory on the impacts of temperature conditions on seasonal abundance of Cx. pipiens, to guide future laboratory or field studies and to highlight in precisely which areas this work is required. By carrying out a simulation study over a range of different temperature regimes predicted for the UK, we examine the effects of both intra-and inter-annual temperature variability on mosquito seasonal abundance to ascertain the climatic conditions that are most likely to lead to high vector abundance and hence to potential disease outbreaks. In particular, we explored the effects of varying mean temperature, amplitude of seasonal fluctuations, the timing of the summer and winter periods, and the sharpness of the summer season. Where possible, these changes were placed in the context of predicted changes for the UK climate. In addition, we use a 30-year UK temperature dataset from the North Kent marshes, a known habitat of Cx. pipiens and Cx. modestus, to investigate sensitivity of vector populations to subtle changes in seasonal temperature regimes as well as changes in mean temperature levels. The analysis highlights the importance of the explicit incorporation of temperature in mosquito abundance models, as we have shown it has a profound effect on mosquito seasonality. By including over winter dynamics in the model, we test the value of the DDE approach in stage structured modelling for understanding dynamics of temperate arthropod vectors.

Mathematical framework
The mosquito life cycle is composed of four main life stages: egg, larval, pupal and adult. Adult females lay rafts of eggs on the surface of pools of water. Eggs hatch into larvae which continue their development to become pupae. Pupae metamorphose into adults which mate and, if female, locate a host to obtain a blood meal and eventually lay more egg rafts. Many temperate mosquito species, including those of the genus Culex, overwinter through inseminated adult females entering diapause. The four stages of the mosquito life cycle can be incorporated into a sequence of DDEs, which are solved to give the abundance of individuals in each life stage through time . Building on the Nisbet and Gurney (1983) framework, one can make the straightforward extension that variation in stage duration may occur as a result of changes in a range abiotic environmental drivers. An outline of the model is presented in this section, with more details on the derivation and implementation in Supplementary Material S.1-S.3. The four state equations which correspond to eggs, E(t), larvae, L(t), pupae, P(t) and adults A(t) at time t, are  where δ () ( = ) ti ELPA ,,, i represents the stage-specific, densityindependent, temperature-driven, mortality rate, δ (() ) π Lt represents the larval mortality rate due to external predation and R i (t) and M i (t) number of individuals recruited into and maturing from stage i, respectively. It was assumed that predation is the main source of density-dependent population regulation. The maturation equations are defined as  with b(t) as the number of eggs laid per female per day at time t and τ ( ) t i , g i (t) and S i (t) as the stage duration, development rate and probability of survival of individuals in stage i (= ) iE L P ,, at time t respectively. Recruitment into the egg stage occurs solely as a result of egg laying by adult females and recruitment into all other stages is a result of maturation of individuals from the stage before. We assume a 1:1 sex ratio between male and female Cx. pipiens (Agnew et al., 2000), which is reflected in the egg recruitment equation. The temperature-dependent growth rate, g i (t), defines the speed at which individuals progress through the stage, allowing the duration of the immature stages to vary continuously. The term . Though the vital rates are temperature-dependent, the equations are written showing a time dependence. The temperature, T, varies as a function of time, such that =( ) TT t . The proportion of individuals which survive from recruitment into one class, to maturation to the next, is defined by the following sequence of DDEs (see Supplementary Material S.1 for derivation), , .
Lastly, the rate of change of the duration of each life stage is given by Here the development rate, g i (t), is dependent on temperature, although the model formulation is sufficiently flexible to allow a range of environmental drivers to be incorporated simultaneously. We solve the system of Eqs. 1(a)-(d) in Fortran 90 using the DDE solver (DDE_SOLVER) written by Thompson and Shampine (2004). The code for the model described here can be found at Ewing et al. (2016). The initial conditions for the equations are described in Supplementary S.2. A full description and derivation of the delay equations required to code the model using a DDE solver is presented in Supplementary Material S.3.

Functional forms of demographic parameters with temperature
We carried out a literature review to identify and parametrise appropriate functional forms for the relationship of development and death rates of Cx. pipiens with temperature. There has been some debate on the exact structure and taxonomy of the Cx. pipiens complex, which is considered by many to be unresolved. A thorough discussion of this debate and of recent decisions made about the taxonomy of Cx. pipiens is presented in Harbach (2012).W h e r e possible we have used data from Cx. pipiens, which is widespread in the UK, but in one instance a lack of available data meant we supplemented this with data from the closely related Cx. modestus (also present in the UK) and Cx. quinquefasciatus which is known to mix with Cx. pipiens in the Southern USA and to be a vector of WNV.
The fitted functional forms representing development, death, and predation follow similar forms as previous models and so their fitting is presented in detail in Supplementary Material S.4. The formulation of the diapause process and its effects on the egglaying rate were newly devised for our model and so are presented in the main text. A table of the fitted parameter values is presented in Supplementary Material S.7. The parameter values were fitted using the nonlinear least squares fitting procedure in the statistical software package R (Core Team, 2013).

Diapause
The proportion of active adults at any point in time, ζ ( ) t , regulates the number of eggs laid per female (() ) bt and is dependent on both photoperiod, ψ ( ) t , and temperature ζ ξψ and ξ S is the spring photoperiod threshold for which greater than 50% of the population have left diapause (used when ψ ( ) t is increasing) and ξ A is the autumn threshold for the population entering diapause (used when ψ ( ) t is decreasing). This means that individuals will enter and remain in diapause provided the photoperiod is below a given threshold and the temperature is below°10 C. At these low temperatures the death rate of adults of very low, which allows for survival through the winter period (details in Supplementary Material S.4). The°10 C temperature threshold was chosen because predicted thermal development thresholds are approximately −°103 4 C (Almirón and Brewer, 1996;Loetti et al., 2011) and photoperiod threshold values were chosen to coincide with induction and termination of diapause (Sulaiman and Service, 1983). The effect of photoperiod on diapause behaviour is strongly dependent on rearing conditions (Vinogradova, 2000). We chose to base our estimates on Sulaiman and Service (1983) because they carried out a UK-based study and we are interested in predicting for the UK. Photoperiod is calculated in hours using the CBM model described in Forsythe et al. (1995), which is a function of latitude. For our simulations we chose the latitude of the North Kent marshes (°) 51 N because it is a habitat for both Cx. pipiens and Cx. molestus (Golding et al., 2012). This area may also run a higher risk of WNV introduction than other parts of the UK, due to possible introduction of mosquitoborne viruses through nearby ports (Golding et al., 2012).

Egg-laying rate
Various factors including type of blood meal (Richards et al., 2012), number of gonotrophic cycles (Richards et al., 2012) and size of the female (Colless and Chellapah, 1960;Cochrane, 1972) have been shown to affect the number of eggs in a raft. For this model we have chosen to use a constant egg raft size of R (200), which we selected to coincide with the average egg raft sizes found in the literature (Jobling, 1938). To quantify the frequency with which egg rafts are laid we divide this egg raft size by the length of the gonotrophic cycle, which is temperature-dependent. The relationship between gonotrophic cycle and temperature was chosen by fitting a function of the following form to the data Vinogradova, 2000) , is the rate of progression of the gonotrophic cycle, at time t. The function means that the gonotrophic cycle halts below°1 0 C, which is consistent with the literature (Vinogradova, 2000). The egg-laying rate per adult female per day, b(t), was then calculated according to

Effects of temperature variation
To determine the effect of changing seasonal temperature on mosquito populations a modified cosine wave was used to approximate the annual temperature curve for the UK (Fig. 1(a)), where μ represents the annual mean, λ the amplitude of annual fluctuations, ϕ the timing of the temperature peak and γ the sharpness of the summer season, which acts as a measure of the  (d)) appropriate for the UK climate. We saw that the distribution of peak temperature dates was centred around the 1st of August, so all ϕ values are reported as þ/ À days from this central point. The summer sharpness parameter, γ, takes values in the range from 0.5 to 2.1 with small values leading to a longer warm period and large values leading to a longer cold period ( Fig. 2(a)). Upon examining the correlation coefficients between these parameters we saw a moderate negative correlation between mean and amplitude (= − μλ r 0.55) i.e. years with lower average temperatures tend to have a larger amplitude of temperature fluctuations. This is observable from the opposing skews of Fig. S.4(a) and (b) and suggests that some of the most extreme high temperatures encountered by the model may not occur in reality, as high mean temperatures and large seasonal fluctuations are unlikely to occur simultaneously. We also saw a strong positive correlation between mean temperature and sharpness of peak (= ) μγ r 0.80 telling us that more peaked summer periods are more likely to coincide with higher temperatures across the year. These temperature parameters were co-varied and effects on mosquito seasonal abundance patterns were analysed and discussed.
Air temperature was used in all parts of the model because information on water temperatures experienced by the aquatic stages was not available. Investigations of the relationship, and discrepancy, between air and water temperature have been carried out, primarily for tropical habitats. However, we are unaware of any appropriate functional forms to link air to water temperature. Furthermore, due to the large range of meteorological factors at play, results cannot be extrapolated from one climatic region to another. Studies examining temperate habitats are rare, though a recent study on British container-breeding mosquitoes by Townroe and Callaghan (2014) found that water temperatures tended to slightly exceed air temperatures, particularly in urban settings, though the difference was generally negligible and never exceeded°2 .5 C at any point in the year. Due to the lack of a defined functional relationship between the two temperatures and the infrequent periods in which there was a substantial difference approximating all temperatures by the air temperature is reasonable for this study.
The effect of covariation of the temperature parameters on seasonal abundance of mosquito populations was described using three mosquito summary statistics ( Fig. 1(b)). Our discussion is focussed around effects on adult mosquitoes, as it is adults which transmit disease. However, example outputs for all life stages are shown in Fig. 2. The peak adult abundance, κ, describes the success of the mosquito population in a given year with a high κ value signalling favourable conditions. Minimum abundance, θ, gives a measure of overwinter survival with small values clearly indicating a small population surviving through winter. Length of the biting season, β,defines the maximum length of time over which the mosquitoes are potentially active in transmission, in the event of a pathogen introduction. Initiation of the biting season occurs when more than 50% of the adult population have left diapause. The end of the biting season is defined as occurring one full life cycle (time required for completion of the gonotrophic cycle and egg, larval and pupal development) before adult emergence of the autumn progeny, which are assumed not to bite before going into diapause i.e. it ends with the biting of the late-summer adults that give rise to this autumn generation. This approximation is made because most Cx. pipiens which overwinter, do so having developed in preparation for diapause through the immature stages and do not take a blood-meal on emergence but instead feed on nectar to increase lipid reserves (Mitchell and Briegel, 1989). It is therefore assumed that the main biting season stops the length of one cycle before these diapausing individuals emerge. Within the biting season, the intensity of transmission will depend upon the ratio of adult vectors to hosts which fluctuates through the season.

UKCIP climate change projections
Temperature profiles used for the simulations were chosen to coincide with the nature of climate warming predicted by the UK Climate Impacts Programme (UKCIP) (UKCIP, 2010). The report details three warming scenarios corresponding to low, medium and high CO 2 emissions, with temperature changes by the 2020s, 2050s and 2080s reported relative to a baseline from 1961-1990(MetOffice, 2010. We present the expected value under the medium emissions scenario for each date range for South East England, the area with the highest estimated risk of WNV introduction from migratory birds (Bessell et al., 2014). Predictions under other climate warming scenarios can easily be obtained by examining the UKCIP projections but are not presented here as we only wish to give an overview, rather than a full region-by-region analysis accounting for a range of emissions scenarios.

Exploring sinusoidal temperature approximation
Many studies incorporate seasonal forcing of either parasite transmission (Altizer et al., 2006) or vector abundance (Smith et al., 2004;McLennan-Smith and Mercer, 2014) through the use of a sinusoidal wave. Indeed, our analysis of how changing temperature affects seasonal abundance of temperate mosquitoes relies on the approximation of the annual temperature profile by a modified cosine wave. We explore how this approximation may affect our predictions of vector abundance and to understand the effects on population parameters of incorporating inter-annual variation in seasonal forcing. We compared results obtained when interpolating between mean daily temperature data points to those found using the modified cosine approximation, both incorporating and excluding variability between years. The model was run for the full 1961-1990 baseline period (with the results from 1961 discarded to allow the equations to stabilise) interpolating between the mean daily temperature values in the North Kent marshes (MetOffice, 2014). The model was then run 29 further times with each year replaced by the best fitting modified cosine wave for that year. Finally, 29 more runs were carried out with each year's daily temperature series replaced by a modified cosine wave which had been fitted to the full baseline time series simultaneously.

Results
An analysis of the effects of changing the temperature parameters on the mosquito summary statisticsmaximum abundance, minimum abundance and length of biting seasonwas carried out to understand how different temperature regimes affected seasonal abundance. These results were also placed into the context of potential changes to the UK climate using the UKCIP predictions. A sensitivity analysis was then carried out to clarify which of the temperature parameters were expected to be the most influential drivers of mosquito abundance. Finally, we examined the impacts of including intra-and inter-annual variability in temperature on mosquito seasonal abundance to help understand its potential influence on epidemiological models.

Maximum abundance (κ)
The effects of changing the temperature profile on maximum abundance, κ, are shown in Fig. 3. In all cases we see that increasing mean daily temperature, μ, leads to an increase in maximum abundance, κ ( Fig. 3(a,b,f)). Similarly, increased amplitude of seasonal temperature fluctuations, λ, consistently leads to increased maximum abundance (Fig. 3(c,d,f)). All cases also show that maximum abundance increases as the peak temperature shifts later in the year (ϕ increases) within the studied range ( Fig. 3(a,c,e)). In all cases the model shows high population growth late in the season, leading to a peak in adult abundance in mid-October. This is because decreasing larval numbers at the end of the year causes a release from density-dependence, resulting in higher larval survival, increased pupal numbers (Fig. 2(c)) and consequently high adult abundance. This effect would not be captured by host-seeking trapping methods as the emergent adults are programmed for diapause (Mitchell and Briegel, 1989) and so would not be caught by traps. Sulaiman and Service (1983) observed mosquitoes entering diapause throughout October in some years, lending some support to this prediction.
Finally, when varied with mean and amplitude, decreased sharpness of the summer season (smaller values of γ), leading to longer periods of high temperature, also give larger maximum abundance ( Fig. 3(b,d)). However, when sharpness of the summer season, γ, and time of peak temperature, ϕ, are co-varied we see that (when considering a late August peak of +20 days) more sharply peaked temperature functions lead to the largest abundances (Figs. 2 and 3(e)). When considering an earlier peak temperature, we see the opposite, with longer periods of high temperature leading to larger abundance.

Minimum abundance (θ)
The effect of environmental conditions on spring starting population size, θ, are shown in Fig. 4. As with peak abundance we see that increasing annual mean temperature or amplitude of fluctuations always led to increased spring population, due to a combination of increased population growth in summer and decreased mortality over shorter winters (Fig. 4(a-d,f)). In co-varying mean and peak temperature timing, the effects of changing peak time were relatively small and difficult to discern in comparison to those of increasing mean and generally appeared to be negligible ( Fig. 4(a)). When varying timing of peak temperature with amplitude we saw increased spring populations as amplitude increased and the temperature peak moved later in the year (Fig. 4  (c)). Similarly, varying peak temperature timing with summer sharpness suggested that later temperature peaks tended to lead to larger spring populations (Fig. 4(e)). We also observe a trade-off between slightly earlier peak temperatures minimising the winter duration and slightly later peaks leading to larger summer   populations. Interestingly, conditions which led to higher summer abundance did not necessarily lead to a higher number of throughwinter survivals. This was particularly evident in examining the relationship between timing of peak temperature and duration of high temperatures as lengthening summer, and correspondingly shortening winter, led to larger spring populations in some cases, whilst also giving smaller summer populations (Fig. 3(c) and 4(c)). Fig. 5 illustrates the effect of temperature on the length of the biting season, β. We see that all four temperature parameters can affect the biting season, provided either the termination or onset of diapause is governed by temperature, via cooler spring or autumn period, rather than via a photoperiodic queue. When photoperiod determines both the start and end of diapause, temperature has a negligible role in determining the length of biting season. Amplitude of seasonal fluctuations was seen to have a smaller effect on the biting season than the other parameters because it had the least potential to affect the duration over which the temperature exceeded°10 C (Fig. 5(c,d,f)). Increasing the mean temperature led to notably longer seasons when the°10 C temperature threshold occurred within the range of the photoperiod thresholds ( Fig. 4(a,b,f)). We see this clearly in the divides between the yellow regions (governed by photoperiod) and the darker regions (governed at least partly by temperature). Shifting the peak temperature later in the year led to a decreased length of biting season because the onset of the season is delayed by low temperatures and the end of the season remains governed by photoperiod for the entire range we explore ( Fig. 4(a,c,e)). Decreasing the sharpness of the temperature peak consistently led to longer biting seasons because the°10 C thresholds move apart and development times generally decrease ( Fig. 4(b,d,e)).

Sensitivity analysis
A sensitivity analysis was carried out to understand the effects of changes within the expected range of variability for each of the four parameters of the modified cosine wave (Fig. 6). We calculated the coefficient of variation for each of the four temperature parameters according to the data from the North Kent marshes and varied the parameters accordingly. The plots clearly show that peak mosquito abundance is most sensitive to changes in mean temperature, μ. This is as expected because both development and survival of larvae increase rapidly with rising temperatures, relative to those currently experienced in the UK. We also note that µ− µ+ λ− λ+ φ− φ+ γ− γ+  6. (a-c) Show the effect of a change in each of mean temperature, μ, amplitude of fluctuations, λ, phase shift, ϕ and sharpness, γ, on peak abundance, κ, minimum abundance, θ and length of biting season, β, respectively. The size of the changes shown are chosen according to the magnitude of the coefficient of variation (+ − / 10.6%, 11.8%, 17.1% and 28.1% for mean temperature, amplitude of fluctuations, timing of peak temperature and sharpness of summer period, respectively).
an increase in mean temperature has a larger effect than an equivalent decrease because development rates are concave upwards ( Fig. S.1). The effect of amplitude, λ, mirrors that of mean, μ, but is less marked because this change does not result in as large an increase in summer temperatures and diapausing mosquitoes experience less sensitivity to small temperature changes. The effect of peak temperature timing, ϕ, is sensitive to which date the shift is relative to. However, we see that in this instance shifting the timing of the peak temperature (a shift of approximately þ/ À 6 days) affects peak mosquito abundance on a similar scale to varying amplitude within its range of variation. Sharpness of the summer season γ, has quite a substantial effect on peak abundance and, akin to our observations with mean temperature, we see that lengthening summer period has a greater impact than an equivalent shortening.
The effects of mean, amplitude and sharpness on spring population size, θ, mimic their effects on peak abundance. Timing of peak temperature, ϕ, is slightly more influential in determination of post-winter population size than on peak abundance because it has a notable effect on the length of the biting season. This is because earlier peak temperatures tend to lead to an earlier biting season start date but the termination of the season occurs as a result of photoperiod. Biting season is affected almost equally by an increase in mean temperature or by a lengthening of the summer season, though decreasing mean has a larger impact than an equivalent shortening of the season. Both amplitude and timing of peak temperature have quite small effects on the length of the biting season, producing changes of about 5% or lower.  show percentage changes, in peak abundance, κ, and minimum abundance, θ, respectively, in predictions from interpolation between mean daily temperatures and from the two modified cosine waves (pinkfixed, blue -variable). (c) shows the actual change in the length of the biting season in moving from the interpolation between mean temperature values and the two cosine waves.

Exploring modified cosine temperature approximation
The mosquito summary statistics (maximum abundance, κ, minimum abundance, θ, and length of biting season, β)w e r e calculated for each year, interpolating between exact temperature values and using cosine waves which either included or excluded inter-annual variability (Fig. 7(a-c)). It is clear from Fig. 7(a) that the cosine wave incorporating inter-annual variability generally predicts a peak abundance, κ, closer to that estimated using the interpolation method than the cosine wave which is fixed across the years. The observation is supported by the fact that the mean absolute percentage difference between the interpolation prediction and the modified cosine prediction is 8.6% (95% CI (5.5%, 11.6%)) when inter-annual variability is included and 14.7% (95% CI ( 11.0%, 18.5%)) when excluded. A Wilcoxon signed ranks test shows that this difference is statistically significant with a p-value of 0.006. The differences in minimum abundance, θ, estimation are less pronounced with an average mean absolute difference of 6.7% (95% CI (5.3%, 8.1%)) when incorporating interannual variability and 8.8% (95% CI (6.6%, 11.0%)) in the fixed case, which was not found to be statistically significant (pvalue ¼ 0.066). When comparing the length of the biting season for the three methods it is clear that the interpolation method leads to greater predictions for the length of the biting season. This is because when we smooth over the true temperature values by fitting a cosine wave the date at which the season begins is pushed later. There was no statistically significant difference found between the estimated lengths of the biting season for the two cosine waves (pvalue ¼ 0.915). The residual sum of squares for each cosine fitw a s not found to be linked with the magnitude of the changes between the interpolation predictions and those from the modified cosine waves. Similarly there was no statistical evidence of a link between any of the individual modified cosine parameters and the percentage changes in abundance.

Discussion
Whilst it has long been understood that a range of environmental drivers, including temperature, have a large impact on mosquito physiology and thus on development, recruitment and death rates (Couret et al., 2014;Ciota et al., 2014;Lyimo et al., 1992), the effect of changing climate on seasonal dynamics has been unclear (Semenza and Menne, 2009). By developing a mechanistic model which explicitly incorporates the effects of temperature on each of the mosquito life stages, we have taken important steps towards understanding how predicted warming can be expected to influence seasonal population abundance of temperate mosquitoes. Our model supports expectations that increased temperatures could lead to substantial increases in mosquito numbers and estimates that, given predicted increases in both mean annual temperature and amplitude of seasonal temperature fluctuations for the UK, Cx. pipiens abundance will increase in coming years. Our results also indicate that even small changes in seasonal temperature, in particular mean annual temperature, could lead to sizeable changes in mosquito seasonal abundance and phenology (Mirski et al., 2012). This sensitivity highlights the need for detailed incorporation of the effects of environment on vector dynamics in epidemiological modelling. In interpreting these findings, it is important to be aware that hydrological processes, such as the effects of rainfall and humidity on adult and larval survival, have yet to be incorporated into the model. Inclusion of these processes will be crucial to fully understand vector seasonality, as high temperatures in temperate climates often lead to lower humidity and rainfall, which could counteract the positive effects of temperature on abundance.
Our results show that the size, both of peak mosquito populations and of populations surviving over winter, will increase with increasing amplitude of seasonal temperature fluctuation. One may have expected that, with increasing amplitude causing both higher summer temperatures and cooler winters (when annual mean temperature remains fixed), effects on peak vector abundance would trade-off against one another. This trade-off is not observed because the model is much more sensitive to temperature changes in the summer months than through the winter. This is partly because there is little empirical data on responses of mosquitoes to cold temperatures, though diapausing mosquitoes are believed to be relatively insensitive to temperatures changes within the region of 0 to°10 C, which is typical for the UK winter. Extreme cold and very low relative humidities have been shown to drastically increase mortality in diapausing Cx. pipiens (Rinehart et al., 2006) but these effects are not included in our model (Eq. S.7) because adults typically overwinter in protective shelters and are therefore likely to be buffered against extreme temperatures and humidities. This buffering was seen in populations of Cx. pipiens, which thrive through much of Russia, where external temperatures ranged from −22 to −°4Cwhilst temperatures in diapause shelters ranged from −11 to −°1.1 C (Vinogradova, 2000). It is rare for temperatures cold enough to induce exceptionally high mortality rates to occur in the UK (3 consecutive days below −°5C in Rinehart et al., 2006). However, even subtle effects of low temperature on overwinter survival have altered disease persistence in the US (Wimberly et al., 2014). Very low temperatures may be limiting for key transmission processes which need to be considered in addition to subtle survival effects when the model is extended to predict disease transmission.
The model also shows that a "late" summer, where temperatures peak in mid-to-late August, leads to a larger peak mosquito abundance than an "early" summer, with peak temperatures in July. We expect this is due to the peak temperature aligning with the development of a larger number of mosquitoes, since abundance in the model steadily increases through the summer months. When peak temperature occurs in late August a shorter, more sharply peaked, period of high temperature leads to higher abundance than a longer warm period, centred at the same time, with the same peak temperature (Fig. 2). This is because more sharply peaked temperature curves delay the exit from diapause due to the cooler temperatures early in the year. As a result, the°10 C threshold required for egglaying (Eq. 3) is reached later in the year, once the photoperiod no longer places restrictions on adult emergence (Eq. 2). Consequently, the exit from diapause is more synchronous, with all individuals leaving immediately once the required temperature is reached. This synchronous exit from diapause causes a large population growth rate early in the biting season ( Fig. 2(b) black line). Conversely, when the temperature peak is broader, the°10 C temperature threshold is reached earlier in the year, when photoperiod still limits adult emergence. This leads to more of a gradual emergence of adults from diapause resulting in a smaller population growth rate in the early part of the year and a smaller population overall ( Fig. 2(b) red line). This highlights that differences in the timing and rate of emergence from diapause may be very influential in determining abundance well beyond the early part of the season. Many standard approaches model seasonality using sinusoidal seasonal forcing to describe parasite transmission and thus fail to capture daily and inter-annual variability (Altizer et al., 2006). We compared seasonal abundance predictions using an interpolation between mean daily temperatures, a series of annually fitted modified cosine waves, and a series of cosine waves ignoring interannual variability. This highlighted that mosquito populations can be very sensitive to small changes in the seasonal temperature profile. We can more closely mimic the results obtained using the observed daily temperatures when incorporating inter-annual variability in cosine curve parameters than we can when excluding such variability (Fig. 7(a)). This finding is crucial because current climate projections for the UK suggest that inter-annual variability in weather patterns will increase in coming years (Jones et al., 2009). We did not find any evidence that failure to either capture temperatures at a daily scale or to include inter-annual variability would lead to systematic over-or under-prediction of seasonal abundance. This contrasts with what we may have expected from Jensen's inequality (Jensen, 1906), which would suggest that by smoothing through the points with a cosine wave we would overestimate development times. Given these results we deduce that the approximation of a modified cosine can still be reasonable for a simulation-based experiment, as presented here, because it does not appear to show bias in results. Indeed, such approaches are necessary if we wish to understand the impacts of specific seasonal temperature characteristics, such as late summers or particularly short summers, as these effects can be difficult to tease out when using observed temperature values. This technique is also useful for forecasting future changes, when precise daily variability in temperatures will be unknown but predictions of how annual means and timings of peak temperatures may shift are available. However, when using the modified cosine wave, estimates of peak mosquito abundance may regularly differ from those we would obtain interpolating between mean daily temperature values by as much as 20% and these margins of error should be considered in the formulation of any transmission model that uses these values as input.
Rainfall patterns are likely to have a direct effect on the number and size of larval habitats, which will affect larval density. Many studies have shown that intra-specific competition, as well as predation, impacts larval development and death rates (Alto et al., 2012;Loetti et al., 2011;Peters and Barbosa, 1977;Beketov and Liess, 2007), with recent work showing that the interaction between temperature and competition may strongly influence seasonal dynamics (Amarasekare and Coutinho, 2014). These impacts will be strongly linked to hydrology and can be included in models using a variety of functional forms discussed by White et al. (2011). The importance of such processes on seasonal abundance was clearly displayed by Morin and Comrie (2013), who showed that both abundance and timing of biting season of Cx. quinquefasciatus across the Southern United States was strongly latitude and longitude dependent, due to the interaction of temperature and rainfall throughout the season. Our model results indicate that population growth is strongly affected by our assumption of a fixed larval habitat, as increased larval density in midsummer causes population growth to level off, due to predation rates approaching their maximum, before increasing again in the latter part of the season, particularly when there is a long summer (Fig. 2 (b) red line). As larval numbers dwindle in the late summer months, density decreases and larval survival increases again, leading to a late season abundance increase. In the field, due to contraction of breeding sites as they dry out over the summer, the crash in larval survival due to predation, and in future extensions of the model competition, may continue into the late summer period. However, larval habitat unrelated to natural precipitation, such as people storing water butts in their gardens or using hosepipes, may also be key when considering abundance of Cx. pipiens in the UK, due to the decreased dependence on rainfall for larval habitat through the drier months, particularly in urban settings (Townroe and Callaghan, 2014). To account for these effects a similar approach to that taken by Tran et al. (2013) could be taken where larval habitat was composed of both permanent and temporary reservoirs. For these reasons, in the future hydrological processes must be explicitly incorporated into the model if we are to more fully understand how changing climate may influence seasonal abundance of temperate mosquitoes.
We have shown that DDEs can be an effective and flexible tool in modelling mosquito seasonal dynamics. To develop and refine this work, advancements must be made in three key areas of research. Firstly, understanding of how predation on larvae and competition between larvae interact to regulate population size must be improved, as these areas have, to date, been studied in isolation of one another (Quiroz-martinez and Rodriguez-Castro, 2007;Madder et al., 1983;Agnew et al., 2000;Alto et al., 2012). Secondly, humidity has been shown to have an effect on mosquito populations (Lebl et al., 2013;Carrieri et al., 2014), however the impacts on individual mosquitoes are poorly understood, which restricts model parametrisation. Finally, the interaction between temperature and photoperiod and its effect on diapause behaviour are strongly dependent on region and rearing conditions (Vinogradova, 2000;Spielman, 2001) and whilst numerous investigations have been carried out to understand these processes in Cx. pipiens populations in Russia and surrounding countries (Vinogradova, 2000) and in North America (Spielman, 2001;Bailey et al., 1982;Madder et al., 1983;Sanburg and Larsen, 1973), this work remains to be done for European populations.
Vector-borne disease models can only be effective if the estimates of vector density, required as inputs, are accurate (Lord, 2004). The extrinsic incubation period (EIP) of a pathogen, mosquito biting rate and vector competence may all vary substantially with temperature (Ciota and Kramer, 2013), meaning the identification of periods where peaks in vector abundance may overlap with high biting rates and short EIPs will be crucial in determining transmission risk. Our model predicts these timings given environmental conditions and shows that vector abundance and phenology are highly variable, which could interact positively or negatively with disease dynamics, causing sporadic disease outbreaks. By explicitly modelling the vector population in all its life stages and incorporating the effects of environmental variability on the vector life cycle we have developed a tool which can be used in tandem with existing epidemiological models to improve seasonal predictions of vector-borne disease transmission, using R 0 models like that presented by Hartemink et al. (2011). Recent studies investigating the risk of WNV introduction or spread, such as Bessell et al. (2014), have been forced to make the assumption that mosquito populations are fixed, without accounting for the effects of environment on mosquito seasonality and behaviour. By coupling our predictions of seasonal vector abundance with these epidemiological models we could substantially improve model predictions and our understanding of disease transmission risks. In particular, from our assumption of a 1:1 sex ratio we can estimate the size of the females population and through the inclusion of photoperiod we segregate the adult population into both active and diapausing mosquitoes. This segregation can be applied to the predicted adult abundances using the temperature and photoperiod predictions to estimate the size of the active population, the key subset of the adults required for epidemiological models.