The importance of temperature fluctuations in understanding mosquito population dynamics and malaria risk

Temperature is a key environmental driver of Anopheles mosquito population dynamics; understanding its central role is important for these malaria vectors. Mosquito population responses to temperature fluctuations, though important across the life history, are poorly understood at a population level. We used stage-structured, temperature-dependent delay-differential equations to conduct a detailed exploration of the impacts of diurnal and annual temperature fluctuations on mosquito population dynamics. The model allows exploration of temperature-driven temporal changes in adult age structure, giving insights into the population’s capacity to vector malaria parasites. Because of temperature-dependent shifts in age structure, the abundance of potentially infectious mosquitoes varies temporally, and does not necessarily mirror the dynamics of the total adult population. In addition to conducting the first comprehensive theoretical exploration of fluctuating temperatures on mosquito population dynamics, we analysed observed temperatures at four locations in Africa covering a range of environmental conditions. We found both temperature and precipitation are needed to explain the observed malaria season in these locations, enhancing our understanding of the drivers of malaria seasonality and how temporal disease risk may shift in response to temperature changes. This approach, tracking both mosquito abundance and age structure, may be a powerful tool for understanding current and future malaria risk.


Supplemental Methods
The equations for the temperature-dependent stage-structured model are shown below for convenience. A full description of the model, assumptions and parameters can be found in Beck-Johnson et al. 2013 [22].
The state equations are (egg (E(t)), larva (L(t)), pupa (P(t)), adult (A(t))): where R i (t) (i = E, L, P, or A) is the recruitment into or out of a state and d i (t) represents the per capita, stage-specific mortality rate. The recruitment into a stage is dependent on the recruitment into the previous stage according to, R A (t) = R P (t t P (t))S P (t) h P (t) h P (t t P (t)) where, t i (t) is the duration of stage i at time t. The egg-laying rate (the number of eggs per female per day) is given by b(t) and the temperature-dependent, stage specific development rate is h i (t).
The length of the delay, t i (t), is determined by the temperature-dependent development rate h i (t) of stage i at time t. S i (t) represents the survival through stage i and expands as follows, The stage-specific per capita mortality rate equations (d i (t)) are given by, where µ j ( j = 0, 1, 2, 3, 4, or 5) is a scalar and g i (i = E, L, or P) is the proportion of the juvenile life cycle that the eggs, larvae and pupae take up, respectively. The extra mortality experienced by the larvae because of density-dependence is denoted by s.

Supplemental Results
The amplitude of the predicted seasonal mosquito population cycle increased as the magnitude of the annual temperature range increased with each of the four annual fluctuations (4, 8, 12, 16 C), around mean temperatures ranging from 16-40 C ( Figure S7 & S8). This pattern was most distinct at the edges of the temperature range, with large amplitude cycles populating the middle range of temperatures only as the seasonality increased. Mean temperatures in the middle of the temperature range experienced smaller amplitude cycles because the fluctuations shifted the model over parameter space where mortality is low (e.g. Figures 1 & 2); as the mean temperature shifted away from the middle of the range, fluctuations moved the parameter space between higher and lower mortality rates, creating large cycles (e.g. Figure 2A). The cycle length corresponded to the annual temperature driver with which we forced the model. Population abundances at some of the temperatures, particularly those in the middle of the tested range, displayed overcompensation cycles caused by density-dependence. Low amplitude overcompensation cycles in the larval population were predicted, which were damped or absent in the adult population when driving the model with either the 8 or 14 C diurnal fluctuation, around mean temperatures ranging from 16 to 40 C. These overcompensation cycles did not occur when mean temperatures were between 22 and 24 C for the 8 C diurnal fluctuation ( Figure S9). This type of low amplitude overcompensation cycle at certain temperatures was discussed in our previous study [22]. The amplitude of the overcompensation cycles was much lower than the amplitude of cycles resulting from temperature fluctuations and will likely be overshadowed in the field by other population perturbations. We will therefore focus on the cycles resulting from temperature variability. The population abundance in general decreased when moving from smaller to larger diurnal temperature ranges. The range of mean temperatures at which mosquito populations are predicted to persist was reduced for both the diurnal temperature range explored when compared to the annual temperature range outputs.
A combination of the results from the separate annual and diurnal runs were produced when the model was driven with each combination of the annual and diurnal temperature ranges (Figures S10-S13). Annual fluctuations occurred, increasing in amplitude as the annual temperature range increased. The temperatures at which populations are predicted to persist corresponded to those seen in runs with only the diurnal driver when the annual temperature ranges were 4 or 8 C. When the annual temperature ranges were 12 or 16 C, the range of temperatures at which populations are predicted to persist was contracted when compared to the diurnal driver. This was seen most dramatically when the largest diurnal range (14 C) was combined with the largest annual temperature range (16 C) ( Figure S13).   Figure S3: Annual population dynamics with seasonal or diurnal fluctuation. Temperature function (dark gray line or band) and the equilibrium abundance of adults (blue line) and potentially infectious mosquitoes (red line) across the two diurnal drivers and the four annual drivers. The abundance is number of mosquitoes per liter of larval habitat. The columns going from left to right are 8, 14 C diurnal temperature range and 4, 8, 12, and 16 C annual temperature range. The rows correspond to the mean temperature with A) 18 C B) 22 C C) 26 C and D) 30 C. The radians correspond with months of the year, going clockwise from 1 to 12. Note the axis scale is different for the temperature function and the mosquito population lines are different. The top left plot shoes the upper limit of the axes for the three lines in the corresponding color; they all start at 0, the temperature function plot goes to 42 C and the two mosquito lines go to 120. The axes are the same for all the radial plots. Figure S4: Annual population dynamics with seasonal and diurnal fluctuation. Temperature function (dark gray band) and the equilibrium abundance of adults (blue line) and potentially infectious mosquitoes (red line) across the combinations of the four seasonal drivers and the 8 C diurnal driver. The abundance is number of mosquitoes per liter of larval habitat. The columns, going from left to right, have an 8 C diurnal temperature range combined with a 4, 8, 12 or 16 C annual temperature range. The rows correspond to the mean temperature with A) 18 C B) 22 C C) 26 C and D) 30 C. The radians correspond with months of the year, going clockwise from 1 to 12. Note the axis scale is different for the temperature function and the mosquito population lines are different. The top left plot shoes the upper limit of the axes for the three lines in the corresponding color; they all start at 0, the temperature function plot goes to 42 C and the two mosquito lines go to 120. The axes are the same for all the radial plots. Figure S5: Annual population dynamics with seasonal and diurnal fluctuation. Temperature function (dark gray band) and the equilibrium abundance of adults (blue line) and potentially infectious mosquitoes (red line) across the combinations of the four seasonal drivers and 14 C diurnal driver. The abundance is number of mosquitoes per liter of larval habitat. The columns, going from left to right, have a 14 C diurnal temperature range combined with a 4, 8, 12 or 16 C annual temperature range. The rows correspond to the mean temperature with A) 18 C B) 22 C C) 26 C and D) 30 C. The radians correspond with months of the year, going clockwise from 1 to 12. Note the axis scale is different for the temperature function and the mosquito population lines are different. The top left plot shoes the upper limit of the axes for the three lines in the corresponding color; they all start at 0, the temperature function plot goes to 42 C and the two mosquito lines go to 120. The axes are the same for all the radial plots.   Figure S7: Equilibrium abundance of model runs using a annual temperature driver. The x-axis is the mean temperature around which the temperature is fluctuated. The vertical bars are the abundances (number of mosquitoes per liter of larval habitat) that the system cycles through at equilibrium. The line going through the bars is the average abundance. Panel A, B and C are the larval, total adult and potentially infectious adult abundance with a 4 C annual temperature range, panel D, E and F are the larval, total adult and potentially infectious adult abundance with an 8 C annual temperature range.

Mean Temperature (ºC)
16º ATR 12º ATR Equilibrium Abundance Figure S8: Equilibrium abundance of model runs using a annual temperature driver. The x-axis is the mean temperature around which the temperature is fluctuated. The vertical bars are the abundances (number of mosquitoes per liter of larval habitat) that the system cycles through at equilibrium. The line going through the bars is the average abundance. Panel A, B, and C are the larval, total adult and potentially infectious adult abundance with a 12 C annual temperature range, panel D, E, and F are the larval, total adult and potentially infectious adult abundance with a 16 C annual temperature range. At some combinations of high temperature and large temperature fluctuations the solver was too stiff to handle the rapid changes in temperature and could not accurately capture the dynamics of the potentially infectious population.

Mean Temperature (ºC)
14º DTR 8º DTR Equilibrium Abundance Figure S9: Equilibrium abundance of model runs using a diurnal temperature driver. The x-axis is the mean temperature around which the temperature is fluctuated. The vertical bars are the abundances (number of mosquitoes per liter of larval habitat) that the system cycles through at equilibrium. The line going through the bars is the average abundance. Panel A, B, and C are the larval, total adult and potentially infectious adult abundance with an 8 C diurnal temperature range, panel D, E, and F are the larval, total adult and potentially infectious adult abundance with a 14 C diurnal temperature range.  Figure S10: Equilibrium abundance of model runs using a both an annual and diurnal temperature driver. The x-axis is the mean temperature around which the temperature is fluctuated. The vertical bars are the abundances (number of mosquitoes per liter of larval habitat) that the system cycles through at equilibrium. The line going through the bars is the average abundance. All panels in this figure have a 4 C annual temperature range. Panel A, B, and C are the larval, total adult and potentially infectious adult abundance with an 8 C diurnal temperature range, panel D, E, and F are the larval, total adult and potentially infectious adult abundance with a 14 C diurnal temperature range. At some combinations of high temperature and large temperature fluctuations the solver was too stiff to handle the rapid changes in temperature and could not accurately capture the dynamics of the potentially infectious population.  Figure S11: Equilibrium abundance of model runs using a both an annual and diurnal temperature driver. The x-axis is the mean temperature around which the temperature is fluctuated. The vertical bars are the abundances (number of mosquitoes per liter of larval habitat) that the system cycles through at equilibrium. The line going through the bars is the average abundance. All panels in this figure have a 8 C annual temperature range. Panel A, B, and C are the larval, total adult and potentially infectious adult abundance with an 8 C diurnal temperature range, panel D, E, and F are the larval, total adult and potentially infectious adult abundance with a 14 C diurnal temperature range. At some combinations of high temperature and large temperature fluctuations the solver was too stiff to handle the rapid changes in temperature and could not accurately capture the dynamics of the potentially infectious population. Figure S12: Equilibrium abundance of model runs using a both an annual and diurnal temperature driver. The x-axis is the mean temperature around which the temperature is fluctuated. The vertical bars are the abundances (number of mosquitoes per liter of larval habitat) that the system cycles through at equilibrium. The line going through the bars is the average abundance. All panels in this figure have a 12 C annual temperature range. Panel A, B, and C are the larval, total adult and potentially infectious adult abundance with an 8 C diurnal temperature range, panel D, E, and F are the larval, total adult and potentially infectious adult abundance with a 14 C diurnal temperature range. At some combinations of high temperature and large temperature fluctuations the solver was too stiff to handle the rapid changes in temperature and could not accurately capture the dynamics of the potentially infectious population.  Figure S13: Equilibrium abundance of model runs using a both an annual and diurnal temperature driver. The x-axis is the mean temperature around which the temperature is fluctuated. The vertical bars are the abundances (number of mosquitoes per liter of larval habitat) that the system cycles through at equilibrium. The line going through the bars is the average abundance. All panels in this figure have a 16 C annual temperature range. Panel A, B, and C are the larval, total adult and potentially infectious adult abundance with an 8 C diurnal temperature range, panel D, E, and F are the larval, total adult and potentially infectious adult abundance with a 14 C diurnal temperature range. At some combinations of high temperature and large temperature fluctuations the solver was too stiff to handle the rapid changes in temperature and could not accurately capture the dynamics of the potentially infectious population.

Victoria Falls, Zimbabwe
Abundance 1969 1975 1979 1985 Time in Years  Figure S15: Mean annual mosquito population dynamics, filtered by water availability. This plot shows the model predictions for the potentially infectious mosquito population for Mozambique from 1980 -1990. The left axis corresponds to the mosquito population abundance (number of mosquitoes per liter of larval habitat) and the right y-axis corresponds with precipitation. The blue dashed line is the mean monthly precipitation and the horizontal green dotted line shows the 80mm precipitation threshold. The mosquito population line has been color coded to correspond with the availability of water. Recall that the model is not directly run with precipitation information, but is rather filtered once the temperature-dependent results are established. The red portions of the line correspond to populations in which the eggs were laid when precipitation was equal to or above the 80mm per month threshold. The orange portions of the mosquito population line shows the populations where precipitation has fallen below the threshold but mosquitoes that were already in the potentially infectious class are surviving. The grayed out portion of the line shows what the population would do if water was available year round. The black horizontal time index bars at the top of the panels mark the beginning and end of the malaria season. Note that the temperature data is much less complete after 1990 so the results from 1990 -2000 are not as robust as the first 10 years.