A dynamically structured matrix population model for insect life histories observed under variable environmental conditions
Renewal processes represent development under variable conditionsThe consequence of a drastic environmental change can be demonstrated by introducing a shift in development time during the process. For demonstration, we consider a scenario where a group of individuals enter into a favourable environment reducing development time from (40pm 5) time units to (20pm 5).We show, in Fig. 1, that our dynamic pseudo-stage-structured MPM yields a gradual stage completion with an average development time of approximately (30pm 5) steps (solid dark lines) when conditions shift at ({tau }=20) (each step corresponds to 1 time unit). The target Erlang-distributed development trajectories without the shift are shown as dashed gray lines. The snapshots of the population structure, represented by the development indicator q, taken at each time step, show that half of the development is complete at the time of the switch and the switch accelerates the accumulation of q (Fig. S1).Figure 1Response to change in development time. The number of developing individuals is simulated by using the cumulative development process and compared to (a) the age-dependent development process, (b) an ODE representation, (c) an LCT representation, and (d) a DDE representation. Solid dark lines show the cumulative development and thick blue lines show the alternative models. Dashed gray lines mark the two target trajectories before and after the shift in development time (marked with red crosses).Full size imageIn age-dependent development, a sharp transition, instead of a gradual one, is observed at the (20^{th}) step (Fig. 1a). The switch results in the majority of individuals reaching target development age immediately at the time of switch. Previous work, reported in Erguler et al.59 and Erguler et al.55, aimed at modelling population dynamics under variable conditions, based on this dynamic age-dependent framework. Our results suggest that cumulative development might improve the fit to the data, prediction accuracy, and applicable geospatial range of these models.We see in Fig. 1b that the canonical ODE framework represents an exponentially distributed development time and a shift in rate at (t=20). The LCT extension to the framework helps to incorporate time dependence and represent the long and short development time distributions (Fig. 1c). The resulting model accommodates change in the rate parameter (gamma ) (Eq. 8), e.g doubling of (gamma ) changes development time from (40pm 5) to (20pm 2.5). However, to accommodate the required shift, the model needs to be transformed from a 66-dimensional system to an 18-dimensional one, which is beyond the scope of this work. We argue that in cases where development time distribution is fixed a priori (excluded from model calibration), the LCT framework provides a significant advantage over canonical ODEs. Although the framework has been used in the field of infectious disease epidemiology64,65, it has recently been applied to the modelling of vector population dynamics30.The DDE framework also yields a gradual development trajectory with an intermediate duration (Fig. 1d). However, the distribution tends towards the longer development trajectory compared to the one achieved with cumulative development. The canonical DDE framework assumes a homogenous cohort, where all individuals react in the same way to variations in development rate. The assumption gives rise to sharp stage transitions within a single generation if all individuals are introduced at the same time. As a potential workaround, it has been proposed to generate a plausible population history, through variable entry times, until the required (or observed) developmental variation builds up31,32. Variation in development rates then acts upon the population and results in the modification of the existing age-structure. It is worthwhile to mention that a recent extension to the DDE framework to accommodate trait variation in population dynamics34 might also accommodate changing development rates within a single stage; however, it has not yet been employed at this scale.Cumulative development is in agreement with the widely known degree-day (DD) framework, where development time is predicted by the heat accumulating in organisms46. Although the rate of accumulation in response to environmental conditions varies considerably, the DD framework implies that the combination of two different rates yields an average development time (also seen with cumulative development in Fig. 1). Experimental evaluation of this will be the topic of future research.It is worth mentioning that our dynamically structured renewal process-based MPM follows the assumption of random population heterogeneity9,11; namely, at the individual level, the future behaviour of an organism is not affected by its historical behaviour. However, trait variation within a population is prevalent in many species, and is known to impact population dynamics and species interactions34,66,67. Future development of our framework will consider improving upon this limitation.Environmental variation transformed into development timesSeveral non-linear relationships have been proposed to represent the temperature dependence of insect development68. A common feature is the presence of low and high temperature thresholds beyond which development is prohibitively slow. Often, there exists an optimum between the thresholds where the process is most efficient. A typical relationship between temperature and development rate, reported in Briere et al.50, is seen in Fig. 2a. Mean development time, given by the reciprocal of rate in Fig. 2b, exhibits the two thresholds and the optimum.Figure 2Development under environmental variation. In (a), development rate (Eq. 9) is shown with (alpha =1.5times 10^{-5}), (T_L=0^oC), and (T_H=50^oC). In (b), mean development time is shown together with the probability densities of three temperature regimes ((rho _L), (rho _M), and (rho _H)). In (c), the number of individuals completing development at each step are shown with respect to the three temperature regimes. Solid lines indicate the median, shaded areas indicate the (90%) range of 1000 simulations, and thick lines indicate simulations with the expected values of each regime.Full size imageTo investigate how temperature variation is transformed into cumulative development time, we assumed three variation regimes at relatively low, medium, and high temperatures ((rho _L), (rho _M), and (rho _H), respectively). Densities of the corresponding Gaussian probability distributions are plotted in Fig. 2b. Accordingly, each variation is transformed by a slightly different region of the rate function (Eq. 9). Eventually, the three development time distributions emerge as shown in Fig. 2c.We found that the output of (rho _H) is skewed towards longer durations compared to what we would otherwise obtain if we simulated the process under constant conditions with the mean of (rho _H). The impact of variation in the middle range, (rho _M), is similar to that of (rho _H), but less pronounced. Conversely, the output of (rho _L) is skewed towards shorter durations. Our results suggest that, when development is already highly efficient, variation in temperature causes frequent encounters of longer (but not shorter) development durations, eventually extending the overall duration of the process. In the low efficiency range, development takes long to complete, but frequent encounters of relatively short durations—especially as the process approaches its optimum duration—triggers completion earlier than in the case of no variation.Overall, our model predictions are in agreement with the rate summation effect, which states that the different outcomes obtained under constant and varying temperatures is due to the non-linear relationship between temperature and development rate (the Kaufmann effect)16. Furthermore, acceleration of development in insects subjected to varying high temperatures, its retardation at varying low temperatures, and low variability of development time in the linear range of the rate curve have been widely discussed69. Several groups have reported evidence in support of this effect, which is also in agreement with our results. For instance, Vangansbeke et al. (2015) reported for three insect species, Phytoseiulus persimilis, Neoseiulus californicus, and Tetranychus urticae, that varying temperatures with a lower mean yields faster development compared to the yield at mean constant temperatures70. However, observations of this phenomenon might result in different responses for different species at similar temperatures due to the difference in rate curves. Identification of the optimum temperature range may facilitate comparison. For instance, Carrington et al. (2013) assumed (26^oC) as optimum based on the high dengue incidence in Thailand, and showed that large variations around (26^oC) increases development time for the dengue vector, Aedes aegypti71. Wu et al. (2015) demonstrated that development is faster at around (26^oC) compared to (23^oC) for the fly, Megaselia scalaris, and found that varying temperatures at around (23^oC) accelerates the process47. Finally, in a modelling study employing DDs, Chen et al. (2013) reported that larger diurnal temperature ranges relate to additional DD accumulation and faster development in grape berry moth, Paralobesia viteana72. Under the realistic non-optimum field conditions, where these simulations had been performed, a decrease in development time is expected in response to varying temperatures according to our results.We note that the variation in development times is due to temperature since we ignore intrinsic stochasticity to demonstrate the impact of (rho ) in isolation. The deterministic setup removes the upper limit in the number of distinct pseudo-stage indicators: a different q emerges from each k, and a different k emerges from each (rho ). Since the number of pseudo-stages quickly exhausts the computational resources, we set the precision of q to the nearest 100(^{th}) decimal point, effectively capping the number of pseudo-stages at 100 (see Accuracy of the pseudo-stage approximation). As shown in Fig. S2, the approximation has a negligible impact on accuracy.Environmental dependency extracted from life tables under constant conditionsHaving discussed the importance of environmental variability in development, in this section, we employ a well-established experimental method to unravel the relationship between temperature and development time in a common mosquito species. In contrast to invasive vectors, which effectively render new territories suitable for disease transmission, Culex species pose an imminent threat with their wide distribution and ornitophilic (Cx. pipiens biotype pipiens), mamophilic (Cx. pipiens biotype molestus), and intermixed (their hybrids) blood feeding behaviour. Here, we investigate the temperature dependencies of mortality and development of Cx. quinquefasciatus, the southern house mosquito, which is an important disease vector, widely distributed across the tropics and sub-tropics73,74.To infer the dependencies, we used a generic temperature-driven insect development model, described in Methods, and the life history observations performed at five constant temperatures (15, 20, 23, 27, and (30,^{circ })C) under laboratory conditions60,61. As a result of the inverse modelling procedure, detailed in Methods, we found that the generic model yields an overall match between the simulations and observations. In Fig. 3a, we present a comparison of observed and simulated maximum production and the stage-emergence times for pupae and adults. Here, we define the stage-emergence time as the time taken from the beginning of an experiment to the time when half of the maximum production of a stage (pupa or adult) is observed. In addition, in Fig. S3, we present the comparison of time trajectories separately for each temperature.Figure 3Inverse modelling of Cx. quinquefasciatus environmental dependency. The comparison of observed and simulated maximum pupa (P) and adult (A) production and the corresponding stage-emergence times is given in (a). Observations are represented with dots and simulations with box plots. The environmental dependency of larva and pupa development time (b) and mortality (c), derived by the posterior mode sample (Theta _q), is shown in (b,c). Solid lines represent the median and shaded areas represent the (90%) range.Full size imageWe found that the generic model faithfully replicates the observed development times of larvae and pupae. On the other hand, stage mortalities are predicted well at three temperatures, but are overestimated at 20 or (27,^{circ })C. The impact of temperature on mortality might be more complex than it is captured by the quartic equation (Eq. 11). Optimum survival seen at (27,^{circ })C suggests that the relationship might be non-symmetrical or multimodal. In addition, the observed variability in mortality suggests that the mismatch could also be due to experimental error or the intrinsic stochasticity of the biological processes.We extracted the functional forms of temperature dependence from the posterior samples, shown in Fig. 3b, c, and found that the data inform the model as expected within the temperature range of the experiments ((15{-}30,^{circ })C). Stage durations are well informed, and reflect the low variability seen in the data (the standard deviation is less than 1.5 days at all temperatures for both stages). Accordingly, pupae develop in less than 4 days, which is much shorter than the larva development time (between 10 and 20 days above (20,^{circ })C). The model predicts that the minimum temperature at which development occurs (from the larva stage) is (10.5,^{circ })C, which is close to (10.9,^{circ })C, reported in Grech et al.75.The observed variability in pupa and adult production suggests that survival is a highly stochastic process regardless of the controlled laboratory conditions. A deterministic model, such as the one used in this context, represents the mean of such processes but does not capture their variability. The simulated variability is a result of the uncertainty in parameter estimates. Model parameters contribute unequally to the output as a result of the model structure and the functional forms of temperature dependence, and the data inform certain parameters better than others76,77. For instance, daily mortality, shown in Fig. 3c, is more constrained for larva than pupa, which is likely due to the short duration of the pupa stage—changes in daily mortality have larger consequences as development time increases.We note that a well-informed model yields predictions in the form of verifiable hypotheses; however, these are not necessarily accurate predictions. Model accuracy is assessed when such hypotheses are experimentally tested as part of the cyclic process of model development78. Here, we demonstrated that our modelling framework can be used to derive biologically meaningful inferences and to help improve the understanding of the temperature dependence of Cx. quinquefasciatus.Greater information content of semi-field experimentsThe number of experiments required to test a range of conditions, including different combinations of multiple drivers, may quickly exhaust available resources. Moreover, variable conditions may have a previously unaccounted impact on development and mortality. In this section, we demonstrate that observations performed under variable conditions are valuable sources of information for our modelling framework, which is capable of representing the dynamics under such conditions.Cx. pipiens, the northern house mosquito, is a competent disease vector, widely distributed across the temperate countries in North America, Europe, Asia, and North and East Africa74,79. Unlike Cx. quinquefasciatus, Cx. pipiens biotype pipiens is known to enter a reproductive diapause phase, where adult females arrest oogenesis during harsh winter conditions80,81. When larvae are exposed to short photoperiods and low temperatures during development, they emerge as adults destined to diapause. Although Cx. pipiens biotype molestus has lost the ability to diapause, its immature stages have been reported to retain metabolic sensitivity to photoperiod82,83.To reveal the environmental dependence of the molestus biotype, we exposed its eggs to variable temperatures in semi-field conditions until adult emergence (or loss of cohort). The numbers of viable larvae, pupae, and adults observed in different experimental batches are given in Fig. S4. We employed the extended model with both temperature and photoperiod dependence (see Methods), and calibrated the model against seven of the semi-field experiments, performed in March, May, June, July, August, and September (Fig. S4(a), (b), (d), (f), (g), (i) and (j)).As a result, we found that the model replicates the patterns of abundance emerging in the observations, e.g. stage timing and maximum adult production, reasonably well in most of the experiments, regardless of the times during which they were performed (Figs. S5 and S6). Quantitative evaluation of the agreement reveals that the observed and simulated adult emergence times are less than a week apart (Table 1).Table 1 Comparison of observed and simulated adult emergence time and the total number of adults produced. Simulation output is given in terms of the median and (90%) range.Full size tableOn the other hand, we found that egg and larva mortalities, and also, pupa and adult production are highly variable in the observations (see Fig. S4(c), (f), and (g)). Spikes of larva mortality are seen in Spring and Autumn (especially in May, September, and October). Despite this variability, the difference between the predicted and observed adult production was around 11 or less, except in the case of the experiment E7, which unexpectedly yielded only one pupa and no adults.We obtain relatively large mismatches when predicting larva abundances, specifically where egg mortality is not predicted well (E5, E7, E8, E10, E11, E12). We hypothesise that the stress associated with rearing lab-grown specimens under variable conditions might elevate egg mortality, induce premature hatching, or affect the survival of the larvae produced. Since egg development starts inside gravid females, i.e. under the optimum conditions of the laboratory, the observable part of development subjected to variable conditions remains mainly the hatching behaviour. Consequently, we observed rapid and synchronous completion of the egg stage in all experiments (see Figs. S5 and S6). Being exposed to a narrow range of temperatures, relatively less information can be obtained on the environmental dependency of the egg stage. As a potential improvement, we recommend that future adaptations of the semi-field experiments consider using field-captured adult female mosquitoes as the source of eggs.In addition to egg mortality, we observed spikes of larva mortality in May (E3), July (E8), and in Autumn (E14, E15, and E16). A likely cause of such transient high mortality is brief temperature shifts towards the extremes. However, the rarity of such events prevents the inverse modelling procedure from adequately capturing their impacts on life processes. As a potential improvement, we recommend that the experiments are performed in overlapping time frames, increasing the likelihood of observing the impact of an extreme event at different times during development. We note that the early decline in larva abundance seen in Autumn could be a result of insufficient food supply due to the increased nutritional requirements. According to the proposed metabolic response to short photoperiods, larvae would require additional food to accumulate fat reserves in preparation for diapause, the state where adult females endure several months without feeding. This implies that development takes longer than it would at long photoperiods when subjected to similar temperature regimes.Using the extended model and the semi-field data, we identified the environmental dependencies shown in Fig. 4. The data informed about the temperature dependency of each life stage as well as the photoperiod dependency of larvae. As expected, the overall variability in the inferred dependencies is higher for Cx. pipiens compared to Cx. quinquefasciatus (Fig. 3). We found that the larva and pupa development times closely match the observations reported by Spanoudis et al.62 at long photoperiods (see Fig. S7). However, the development times reported in Kiarie-Makara et al.84 at short photoperiods and moderate temperatures do not suggest a significant impact of daylight, which could be due to the particular strain of Cx. pipiens used in these experiments. As expected, the temperature dependency of egg development was not well informed by the data in the current configuration of the model and the functional forms of environmental dependence.Figure 4Environmental dependency of Cx. pipiens development and mortality inferred from semi-field life table experiments. Solid lines represent the median and shaded areas represent the (90%) range.Full size imageWe found that the photoperiod dependency is significantly non-linear with an average threshold of 13.7 hours of daylight (Fig. 4c). Photoperiod-driven extension in development time (about 1.7 times more at 13:11 h L:D than at 15:9 h L:D) contributes to improving the accuracy of predictions at the end of the high season (Fig. S8). The critical photoperiod (CPP) agrees well with the ones identified for Cx. pipiens biotype pipiens85,86. For instance, Sanburg and Larsen reported that there is an exponential relationship between follicle sizes in adult females (signifying commitment to diapause) and the photoperiods they were exposed to during immature stages85. We inferred a similar (but reverse) gradient between photoperiod and the extension of larva development time from 15 to 12 hours of daylight (Fig. 4c).Risk assessment with annual development profilesWe extrapolated the development dynamics of Cx. pipiens over the calendar year by setting up a hypothetical experiment at the beginning of each week. We simulated the subsequent development dynamics and obtained the annual development profile as shown in Fig. 5. Accordingly, the immature stages begin development from late February and the first adults emerge in May (adults emerging late in May start developing in the experiments set up late in March). The profile is consistent with the regular Cx. pipiens high season in the region.Figure 5Annual development profile of Cx. pipiens in Petrovaradin, Serbia, in 2017. The outcome of each hypothetical semi-field experiment is plotted vertically along the y-axis at the date when the experiment is initiated. The maximum number of adults produced is given in blue, and the time it takes (from the date indicated on the x-axis) to produce half of the maximum is given in green. Solid lines represent the median and shaded areas represent the 90% range of model predictions. Outcomes of the semi-field experiments (dots) are plotted together with the model predictions. The time points marked with circles indicate the experiments used to calibrate the model. Estimated time of first adult emergence is given in the inset.Full size imageAs seen in Fig. 5, predicted adult emergence times agree well with the observations throughout the high season. However, there is a greater variation in the maximum number of adults than the times of emergence (extending to almost (40%) of the possible outcomes in early August). A greater variability (almost (80%) in August) is seen in the corresponding observations, which we transformed into the percentage of eggs emerging as adults (where available) to facilitate comparison. According to the model, variation in adult production is associated with the variation in both development times and mortality during immature stages. We recall that the uncertainty in the informed environmental dependencies is high around relatively less frequently encountered values—especially the lower and higher temperature extremes (Fig. 4). Specifically, egg development times cannot be identified precisely, but immediate hatching of the larvae is predicted between 20 and 25 °C. Consequently, we found that frequent exposure to temperatures outside the well-informed range have a significant impact on the variation in adult production (Fig. S9).We adopt the time of first adult emergence as a proxy of the first generation of adults in the season. According to our model, early adult emergence is a result of shorter development times and higher success rates, which indicates that the temperature conditions allow for an early first generation of adults. An early first generation greatly contributes to an early peak of adult abundance, which may increase the risk of vector-borne disease transmission in humans. For instance, an early peak of abundance may cause an early start of West Nile virus circulation and amplification in Culex pipiens and their avian hosts, which increases the likelihood of virus spillover to humans51,87. Anecdotal evidence shows that the anomalously hot April and May that occurred in 2018 in Serbia shifted the peak of Cx. pipiens abundance forward by more than one month (Petrić et al., unpublished). Similarly, 2018 was the year with the largest number of autochthonous West Nile virus infections throughout Europe (more than the total of the previous seven years together)88,89.In summary, our results showed that the semi-field experiments, when used in combination with our dynamic pseudo-stage-structured MPM, help to develop predictive models and inform over a wide range of environmental conditions. We developed a predictive model of Cx. pipiens biotype molestus development and gained insights into the specifics of temperature and photoperiod dependencies by reducing the need of extensive laboratory data. We used life history observations from 7 experiments performed under semi-field conditions and employed a generic model structure, largely uninformed on the specific environmental dependencies of the species. The cumulative development framework we introduced applies broadly to poikilotherms subjected to highly variable environmental conditions. Although the generic model structure helps to develop exploratory models and identify potential environmental dependencies, accuracy can be improved by customising the models for the known dependencies of particular species. With a straightforward extension of the development model to cover the complete life cycle (with egg laying and density dependence), it is possible to incorporate field observations of eggs or adult mosquitoes, and develop an environment-driven population dynamics model. More
