Optimal control of irrupting pest populations in a climate-driven ecosystem

Consumer–resource model

The underlying consumer-resource model is the best-fit model developed by Holland et al. (2015). Relative mouse abundance M(t) is quantified by an index: captures per 100 trap nights (C/100TN) (Ruscoe, Goldsmith & Choquenot, 2001). The rate of change of M over time t (years) is given by (1)${\displaystyle \frac{dM}{dt}=\left(\alpha g\left(F\right)-{\mu }_1-{\mu }_{2}M-B\left(t\right)\right)M,}$

Food availability F (seeds m⁻²) is predicted by the functional response g(F) and α is the demographic efficiency of mice (i.e., efficiency at converting food into recruitment for the mouse population). The total mortality rate is μ₁ − μ₂M − B(t), where the parameters μ₁ and μ₂ are density-independent and density-dependent mortality rates respectively. They may both be positive or negative depending on non-food related processes, e.g., predation, social interactions, Allee effects. In this paper, we extend the previously published model by adding B(t), which is a time-dependent, density-independent rate of mortality due to control by bait application (see below).

The original model tested four candidate models for the food availability functional response. The best fit was a Holling II (Ivlev) function where c (seeds m⁻² mouse⁻¹ year⁻¹) is the maximum per capita feeding rate and e ((seeds m⁻²)⁻¹) is a measure of foraging efficiency: (2)${\displaystyle g\left(F\right)=c\left(1-exp\left(-eF\right)\right).}$

The rate of change of available food over time is modelled by (3)${\displaystyle \frac{dF}{dt}=S\left(t\right)-hF-g\left(F\right)M,}$

where the second term, hF, describes the change in available food that happens throughout the year at a constant rate h (year⁻¹) unrelated to mouse density and the third term, $g\left(F\right)M$, describes the rate of seed consumption by mice. The first term, S(t), describes the rate at which food is delivered as a function of time, with (4)${\displaystyle S\left(t\right)=\left\{\begin{array}{c}\frac{F_y}{0.25}\mathrm{if}0\le t-\text{floor}\left(t\right)<0.25\hfill \\ 0\text{otherwise}.\hfill \end{array}\right.}$

The floor function rounds t down to the largest integer smaller than t. Thus, during the yth annual cycle, a total amount of seed F_y is produced, which is delivered at a constant rate over the first quarter of the year (nominally February–April). At the start of each year we set $F\left(t\right)=0$, i.e., seed is not carried over between years. The annual seedfall amount F_y, was determined by a climate induced seedfall model (see below). All parameter values are given in Table 1 and were those determined as best-fit parameter values by Holland et al. (2015) (Table 2 in Holland et al. (2015)). These were chosen by modelling mouse density over 25 years using observed annual seedfall data from the OV (starting February 1972) as the annual values of seed F_y. Mouse density at the start of each quarter was extracted from the continuous-time model prediction. These predicted values were compared to observed quarterly mouse density data from the OV collected over the same time period (quarterly, February 1972 –November 1996) and best-fit parameter estimates chosen by minimising the root mean square error. Holland et al. (2015) showed that with these best-fit parameter values the model predicted all major outbreaks in mouse density occurring in the 25 year observed mouse density (C/100TN) time series, although it tended to slightly under-predict the magnitude of outbreaks.

During a control year, bait is applied as an impulse function at a specific time point $t_i^{\ast }$, such that $t-t_i^{\ast }$ describes the time since the ith bait application. After application, bait degrades according to a simple decay function $B_{0}exp\left(-d\left(t-t_i^{\ast }\right)\right)$, where d is the decay rate. The constant B₀ governs peak bait availability, i.e., at the time of application ($t-t_i^{\ast }=0$). Therefore, bait availability $b_i\left(t\right)$ from the ith application is described by a piecewise function (5)${\displaystyle b_i\left(t\right)=\left\{\begin{array}{c}{B}_{0}exp\left(-d\left(t-t_i^{\ast }\right)\right),\mathrm{if}{t}_i^{\ast }\le t<t_{i+1}^{\ast }\hfill \\ 0,\text{otherwise},\hfill \end{array}\right.}$

and overall bait availability B(t) at time t is given by (6)${\displaystyle B\left(t\right)=\sum _{i=1}^{n-1}{b}_i\left(t\right),}$

for n bait applications. In the absence of mice, bait is considered to be inactive after one month, so we choose d = 50 (meaning that b_i(t) has decayed to <2% of its original size one month after the ith application, and <0.02% of its original size after two months. It is presumed that, compared to this decay rate, the effect of mouse predation on bait levels is negligible. Note that the actual value of B₀ is defined later in terms of the kill success rate.

Climate-induced seedfall model

A 1,000 year normally distributed temperature time series T₁, T₂, …, T_1,000, was generated where (7)${\displaystyle T_y\sim N\left(14,1\right).}$

This represents historical mean summer temperatures (daily average for the three month period January to March) in the Orongorongo Valley (hereafter OV), 1972–2014. Randomly generating time-series for the OV in this way was shown to be a valid approach by Holland & James (2015). From this, a time series was calculated using the ΔT model of Kelly et al. (2013), where (8)${\displaystyle \Delta T_y=T_y-T_{y-1}.}$

Annual seedfall predictions were made using the following linear relationship fitted to observed OV beech seedfall data (1972–1996) by Holland et al. (2015): (9)${\displaystyle {log}_{10}{F}_y=0.33+0.97\Delta T_y+{ϵ}_y.}$

The noise term was chosen to have distribution ϵ_y ∼ N(0, 1.3) to give a correlation between change in temperature and seedfall of r² ≈ 0.54 corresponding to the findings of Kelly et al. (2013). These seedfall time series were used as annual inputs of F_y to the mouse model above. Mouse density (C/100TN) was simulated for each of the control scenarios above, with $M\left(0\right)=1.0$.

Plague definitions

We define the mouse density plague level, M_p = 2.02, to be the maximum mouse density in a median seedfall year if the initial mouse density is one, i.e., $M\left(0\right)=1$, (see Fig. 1). Using this definition, if the mouse density was 1 at the start of the year ($M\left(0\right)=1$) 50% of years would be defined as plague years i.e., the mouse density ‘just’ reached plague level at some point during the year. In a longer time series where no control measures are applied (Fig. 2) and the density is continuous across years, i.e., mouse density at the start of the year is the density at the end of the previous year, the plague level definition is not changed and the proportion of years that are plague years is much higher, 85%. Specific thresholds for what constitutes a mouse ‘plague’ or eruption in New Zealand have not been formally defined in the literature; it is difficult to measure exact population densities, and it is not known exactly at what threshold level mice may have an impact on native biodiversity. We therefore use the term ‘plague’, to mean greater than average population abundance for an extended period of time, i.e., demonstrably not an undetectable population and therefore likely to have some impact, and a convenient reference point with which to compare scenarios. In the 25 year time-series of observed mouse abundance from the OV, when mouse abundance was above 2.02 the population tended to be undergoing one of seven larger eruptions. In addition, mouse abundance was above the 2.02 level in 80% of years, suggesting that our definition for plague level and estimate of proportion of plague years are reasonable in this context. We also define plague time, the proportion of time during which the mouse density is above the plague level. For example, in Fig. 1 the plague time for the mast year trajectory is 0.89. In the long term time series of Fig. 2 the plague time is 0.71. Finally we define the plague size, the highest mouse density during the plague period. In the Fig. 1 example, the plague size for the mast trajectory is 6.93. In the long term series (Fig. 2), the expected plague size, given that a plague occurs, is 5.7. Higher values of M_p could be used with the same plague definitions given here, with qualitatively similar results.

Control definitions

The strength of the control impulse is governed by the parameter B₀, which is the value of B(t) at the time of bait application $t=t_i^{\ast }$ (Eqs. (5)–(6)). The absolute value of the control impulse is of little practical use in modelling terms, though operationally it relates to control effort and a parameter value could be calibrated for a given operation. A more appropriate measure of control size is the control success, typically assessed in terms of percentage kill, defined here as the relative decrease in the mouse density one month after the control impulse. The control success will change depending on the current mast level and mouse density but percentage kill is still a useful and widely used measure (Elliott & Kemp, 2016; Innes et al., 1995). In the example of Fig. 1 (black dotted line), where B₀ = 150 and control is applied at the start of the mast season, control success is 88%. Even with this relatively high level of control (consistent with historic broad-scale aerial poison operations for rodents (Innes et al., 1995)), the mouse density at the end of the year remains high as the population recovers, and the mouse density rises above the plague level during September. However, both the plague time (0.31) and the plague size (2.18) are much lower than in the uncontrolled mast year. Crucially, for the success of native species that rely on beech mast for food, the seedfall level is much higher in the controlled scenario compared to the uncontrolled.

Scenarios

To obtain precise estimates for expected plague time, plague size and proportion of plague years under different control scenarios, we ran the extended model (Eqs. (1)–(6)) or a 1,000 year long simulation with M(0) = 1, using a seedfall time series generated with (Eqs. (8)–(9)) from a temperature time series generated using Eq. (7). In this longer simulation, the initial mouse density is not reset to 1 each year but continues with the value at the end of the previous year. Seedfall is reset to zero at the beginning of each year to reflect the inedibility of the previous year’s seed. Note that all the control scenarios used the same underlying seedfall series.

Optimal control timing

Initially, we assume that control can only be applied once each year and we examine three levels of control: low, medium and high, B₀ = 50, 100, 150, respectively. For each level of control, we used the 1,000 year weather/seedfall time series, applying control at a range of times throughout the year. We calculated the proportion of plague years, proportion of time above plague level (i.e., plague time), and average plague size for each control application time.

Tailored control

We now consider how other control options with similar costs could offer greater benefit if tailored. An alternative to annual control could be a biannual-biennial regime, i.e., control is applied twice a year every two years. In some cases this may be more cost effective as resources only need to be acquired every other year; it may also be socially advantageous if the application of poison is somewhat controversial and its use needs to be limited. As operational costs may be considerable (particularly labour and transport/flight-time), and the amount of bait applied contributes relatively little to the overall cost of each control dose, we do not consider control via very frequent smaller doses here. We used the 1,000 year seedfall simulation to compare three biannual-biennial regimes with no control and the annual control regime described above:

1. Regular biennial control: Control occurs every second year and is executed in early September and a month later in early October, i.e., straddling the optimal control period (see ‘Results’).

2. Seedfall determined control: Control occurs in years when seedfall is above the median (i.e., with the same long term average frequency as regular biennial control), in early September and early October.

3. Climate determined control: Control occurs in years when the seedfall temperature driver ΔT is above the median (i.e., with the same long term average frequency as regular biennial control but with more opportunity for error in true/false seedfall prediction), in early September and early October.

Regime two could be used if the seedfall could be measured early enough to plan a control application that year. In cases where this was too late to muster a control application then regime three may be of use.

Optimal control timing

As the control time is changed across the year the plague time varies (Fig. 3A). The different control levels (low, medium and high) have a much stronger effect on the plague time than the control timing. For example, under the low control regime, where the control success rate is between 50 and 60%, the optimal control timing to reduce plague time is spring (mid-September) before the mast season starts when plague time is 0.36. Conversely the lowest reduction is seen in autumn (late April) at the end of the mast season which reduces the plague time to 0.42. The optimal control time at medium and high control levels is also around mid-September. Medium control, which has an 80–85% success rate, reduces the plague time to 0.20, while high control has a 92–94% success rate and reduces the plague time to 0.045.

In contrast, the expected plague size (Fig. 3B) appears to be strongly affected by control timing, in particular for low and medium control levels. The optimal timing to reduce plague time gives the least reduction in expected plague size. Initially this seems counter-intuitive but less so after recalling the definition of expected plague size (the maximum mouse density given that there is a plague). On examination of the time series (Fig. 4) for control during the mast season (March—Figs. 4A, 4C, 4E) and after (September—Figs. 4B, 4D, 4F) we see that the small reduction in plague time between the control timings comes from the years where the plague threshold is only just exceeded. For low control effort both timings have relatively little impact on the maximum mouse density reached for large plagues, but in post mast control (September) the timing is appropriate to reduce mouse density to below the plague threshold for smaller plagues. This means that, for low level control post mast season small plagues are often avoided but if control is during the mast the small plagues still occur. Larger plagues still occur and reach similar maximum densities as in the absence of control for both timings. When calculating an expected plague size over the entire time-series, smaller plagues will reduce the average. Therefore, while the larger expected plague size for low control may seem counter-intuitive, the reduction in small plagues driving this is actually a desirable outcome. Using an alternative metric of plague severity, for example the expected density over the entire time-series, loses this subtlety and gives results similar to the plague time metric.

Tailored control

Seedfall- and climate-determined biennial-biannual control regimes are more effective in reducing proportion of plague years and time above plague levels than a regular biennial-biannual control regime (Table 2). For both the low and medium control levels a seedfall driven biennial-biannual control regime is more effective than annual control (Table 2). If the control trigger is based on climatic variables rather than actual seedfall, the biennial-biannual regime is slightly less effective than annual control, though the differences in plague time between these methods is relatively small.