Structural sensitivity in the functional responses of predator-prey models

To determine if the Leslie-Gower-May model is structurally sensitive, we plot a bifurcation diagram (shown in Figure 3(a)) using the carrying capacity , K, as the bifurcation parameter for each of the sub-models produced using the three functional responses (3) in model (2), and note the differences and similarities between them.

We find that the Leslie-Gower-May model is structurally sensitive: the Hopf bifurcation occurs at a very different carrying capacity for each sub-model. Figure 3(a) shows that the coexistence steady state in the model using the Holling functional response destabilizes at $K=5.55$, in the Ivlev model at $K=12.09$, and in the trigonometric model at $K=26.36$. Additionally, the Ivlev and trigonometric models have saddle node bifurcations of limit cycles at $K=10.51$ and $K=12.94$, respectively, and this bifurcation is not present in the Holling model. Following the formation of these limit cycles, bistability is observed in the Ivlev and trigonometric models. The Ivlev model is bistable over prey density values $K\in[10.51,12.09]$ and the trigonometric model is bistable over prey density values $K\in[12.09,26.36]$.

The addition of stochasticity to the functional responses makes them more like real data and increases their overlap. We add stochasticity to models (1) and (2) by rewriting the model equations in the following form:

i.e., $\xi(t)$ is an Ornstein-Uhlenbeck process and $(x,y)$ solves a system of ODEs driven by $\xi$. We use $\sigma$ to alter the amount of noise included in each time step. $F$ is the right-hand-side of equations 1a and 2a (depending on the model in question) with $f_{i}(x)+\xi(t)$ in place of $f_{i}(x)$. Similarly, $G$ is the right-hand-side of equation 1b with $f_{i}(x)+\xi(t)$ in place of $f_{i}(x)$. Note that $f_{i}(x)$ does not appear in equation 2b, so it is unchanged. Hence, model 4 is a system of one stochastic differential equation (SDE) and two ordinary differential equations (ODEs). At each time step, we sample a Gaussian random variable, $dW=\sqrt{dt}Z$ where $Z\sim N(0,1)$, and solve the SDE (equation 4a) using the Euler-Maruyama scheme to obtain $\xi(t)$. Then, we plug $\xi(t)$ into the system of ODEs and numerically solve the system using Euler’s method with a step size small enough to achieve stable solutions. Note that if the solutions result in prey or predator densities less than zero due to the stochastic element, then the prey or predator density is set to zero. For the Rosenzweig-MacArthur model, we use $\sigma=0.01$ and for the Leslie-Gower-May model, we use $\sigma=50$ unless specified otherwise. Sample stochastic functional responses are given in Figures 1(b) and S1b. The corresponding results are given in Section 4.1.

The behaviour of both the Rosenzweig-MacArthur and the Leslie-Gower-May models is stable coexistence for values of the carrying capacity, $K$, less than the destabilizing bifurcation. We hypothesize that if the model behaviour is more similar leading up to the destabilizing bifurcation, then perhaps the destabilizing bifurcation will occur for more similar values of the carrying capacity. If our hypothesis is correct, it suggests that a key portion of the functional response is in the vicinity of the coexistence steady state.

We choose an interval containing the prey density values of the coexistence steady state to perform the fit. We call this interval the fitting region and choose it such that it is smaller than the original $[0,4]$ fitting region. Note that we choose this region regardless of whether the coexistence steady state is stable or unstable. In models that do not use data, we choose a functional response to be fixed. We will refer to the sub-model where Holling is chosen as the fixed functional response as the fixed Holling case and the other sub-models follow the same naming scheme. Note that in the applied case of this method, we do not have to choose a fixed functional response since we can fit the functional responses to data. We then apply nonlinear least squares to the functional responses fit to either the fixed functional response or to data over the fitting region. The average prey density values of the coexistence steady states in the Rosenzweig-MacArthur and Leslie-Gower-May models are $x=0.05$ and $x=0.65$, respectively. So, the corresponding fitting regions are taken to be $x\in[0,0.1]$ and $x\in[0.3,1]$. The parameter values obtained from this fit to the fixed Holling functional response are given in Table 2. The fixed Ivlev and fixed trigonometric functional responses are given in Figures S1c and S4. The corresponding parameter values are given in Table S1. The Holling fitted functional response results are given in Section 4.2.

Here, rather than fitting the functional responses to each other, we modify the functional responses so that they are exactly the same over portions of the prey density domain $[0,4]$. Each functional response has three intersection points with each of the other two functional responses. All three functional responses intersect at (0,0), since there is zero predation at zero prey density. The values for the non-trivial intersections of the Rosenzweig-MacArthur and the Leslie-Gower-May functional response are given in Table 3. As a result, we can consider piecewise functional responses made up of contributions from two functional responses. These piecewise functional responses give us another way to ask if there is a key part of the functional response that determines the model behaviour. That is, we can examine how each portion of the functional response, (i.e., low, medium, and high prey density) affects the model’s bifurcation structure.

We consider three cases for each functional response for each model (18 cases total). To form each case we choose two functional responses and interchange the functional response we use for each of the three regions. The intersections noted earlier become the transition points between the regions. Sample piecewise functional responses are given in Figures 1(d) and S1d. In these figures, we introduce notation that is used from here forward to describe the piecewise identical sub-models. This notation is a three letter code that corresponds to each of the three prey density regions we have chosen to investigate, with Holling, Ivlev, and trigonometric functions each represented by their first letter. For example, HII represents the functional response that uses the Holling functional response for low prey densities, and uses the Ivlev functional response for medium and high prey densities. The corresponding results are given in Section 4.3.

In the case of the Rosenzweig-MacArthur model, we added Gaussian stochasticity with $\sigma=0.01$. We choose this amplitude because stochasticity of smaller amplitude had virtually no effect on the bifurcation structure. $\sigma=0.01$ stochasticity adds some variance to the time series but does not significantly alter the qualitative behaviour for carrying capacities less than the Hopf bifurcation in the deterministic model. That is, given initial conditions that lead to the stable coexistence steady state in the deterministic case, the stochastic model also results in trajectories that converge to the coexistence steady state. For values of the carrying capacity greater than that of the bifurcation in the deterministic model, the lower side of the limit cycle reaches very low prey densities of approximately $10^{-60}$ at times ($10^{-58}$ times smaller than the stochasticity). As a result, the stochasticity causes extinction within 100 time steps in most of the time series of the stochastic model for carrying capacities greater than that of the deterministic model’s destabilizing bifurcation. These results mean that we are unable to clearly observe the stochastic behaviour and we cannot conclude whether or not stochasticity reduces structural sensitivity in the Rosenzweig-MacArthur model.

In the case of the Leslie-Gower-May model, we added Gaussian stochasticity with $\sigma=50$. As in the case of the Rosenzweig-MacArthur model, we chose this amplitude of stochasticity because it is large enough to affect the bifurcation structure, but not so large that all of trajectories rapidly terminate in extinction of the predator or both populations. Results obtained from both smaller and larger amplitudes of stochasticity are provided in the Supplementary Material (Figure S3).

The bifurcation diagrams (quasi-steady state behaviour) for the stochastic Leslie-Gower-May model are given in Figure 3(b). This figure shows better similarity in the bifurcation structures with the addition of stochasticity. The bifurcation structures are obtained by plotting the maximum and minimum of the stochastic time series (after removal of transients). Thus, these windows show the smallest value of the carrying capacity for which the coexistence steady states may destabilize. Limit cycles in our stochastic model appear to have greater amplitude than those of the deterministic model. We find that the stochastic bifurcation structures show better similarity than the deterministic model’s bifurcation structures. This similarity suggests that there is a decrease in structural sensitivity of the model.

Plots of time series (Figure S2) indicate that bistability remains a property of the stochastic Ivlev and trigonometric models. In the bistable regions, the solution trajectories of the Ivlev and trigonometric models alternate randomly between the coexistence steady state and the stable limit cycle, and also spend some time at the unstable limit cycle. This ability to alternate between stable and unstable steady states is unique to the stochastic model, and has been observed in other work (Sharma et al., 2015; Abbott and Nolting, 2017). These observations occur for small to intermediate stochastic amplitudes, larger amplitudes lead to rapid extinction, and more generally, to the stochastic behaviour washing out the deterministic structure.

We only consider the fixed Holling case here as the fixed Ivlev and fixed trigonometric cases are similar and given in the Supplementary Material (Figures S4 and S5).

The Rosenzweig-MacArthur fitted functional response sub-models yield the bifurcation plots shown in Figure 2(b). The bifurcation values for these sub-models are given in Table 4. These fits result in better similarity of the bifurcation structures and less structural sensitivity as compared to the original model (Figure 2(a)). For example, when the Ivlev functional response is fit to the fixed Holling functional response, the Ivlev fit to Holling sub-model shows behaviour similar to the Holling only sub-model as there is a leftward shift in the destabilizing Hopf bifurcation. This result holds across each of the fits. Hence, considering the prey densities in the chosen fitting region appears to be a better metric to describe functional responses in the Rosenzweig-MacArthur model; in the sense that model behaviour is less sensitive to small changes in the functional response using this method.

The Leslie-Gower-May fitted functional response sub-models yield the bifurcation diagrams given in Figure 3(c) and the value of the bifurcation parameter for each bifurcation is given in Table 4. The resulting sub-models show better similarity in bifurcation structure when compared to the fixed functional response. For example, Ivlev fit to Holling shows a leftward shift of the Hopf bifurcation as compared to the Ivlev only bifurcation structure and the disappearance of the saddle bifurcation of limit cycles. So, the Ivlev fit to Holling sub-model shows behaviour similar to the Holling only sub-model. The same result holds for each of the other cases, and there is better similarity in model behaviour when compared to the fixed functional response. Therefore, as in the Rosenzweig-MacArthur model, fitting the functional responses over the chosen fitting region appears to decrease structural sensitivity in the Leslie-Gower-May model.

We would like to note that these theoretical results are subjective to which functional response is chosen to be fixed. That is, the model behaviour will always become more similar to the fixed model. In the case of the results shown in this section, the model behaviour always becomes more similar to the Holling sub-model behaviour. Our method when applied to a model that uses data to fit the functional responses, however, is not subjective. We will discuss this application of our method to data further in Section 5.2.

Suppose that medium prey densities determine the bifurcation structure of a model. Then, according to our hypothesis, we expect YXY and ZXZ to show model behaviours more similar to the X model behaviour (be that H, I, or T). Based on the results in Section 4.2, we might actually expect low prey densities to be the most important, in which case we expect XYY and XZZ to have better similarity to the X model behaviour. However, our results show that none of the three regions we consider separately determine the model behaviour. In this section, we describe only the piecewise identical functional responses made up of the Holling and trigonometric functional responses. The bifurcation diagrams for these cases of the Rosenzweig-MacArthur and Leslie-Gower-May sub-models are given in Figures 2(c) and 3(d). The rest of the cases are similar and are given in the Supplementary Material (Figures S6-S8).

In the case of the Rosenzweig-MacArthur model, there is not better similarity of the Holling to trigonometric sub-models (i.e., sub-models of the form HTT, THT, or HHT) behaviour towards the trigonometric behaviour. This result suggests that both the first and second regions determine the resulting bifurcation behaviour of the Rosenzweig-MacArthur model and that these regions are the most structurally sensitive.

However, the Rosenzweig-MacArthur trigonometric to Holling piecewise identical sub-models (i.e., sub-models of the form THH, HTH, or TTH) are slightly different. In the THH plot, there is a small leftward shift of the saddle bifurcation of limit cycles from the original value (Table 4) at $K=2.644$ to $K=2.431$. There is, however, no change in the location of the Hopf bifurcation. In the HTH plot, the same leftward shift occurs. Lastly, the TTH plot has a resulting bifurcation structure that matches that of the trigonometric only functional response model behaviour. Overall, the trigonometric functional response has a strong effect on the bifurcation structure. Further, when the Holling functional response takes up two of the three regions (in the THH and HTH cases), the Holling functional response has an effect, though limited, on the bifurcation structure. This effect causes the aforementioned shift in the saddle bifurcation of limit cycles. The piecewise identical functional responses seem to affect the resulting model behaviour without following any clear rules. Hence, we are unable to make any claims on which prey density regions determine the Rosenzweig-MacArthur model behaviour.

In the case of the Leslie-Gower-May model, the HTT and THT sub-models show better similarity when compared to the trigonometric bifurcation structure while the resulting bifurcation structure for the HHT sub-model matches that of the Holling only functional response. We found that in the Leslie-Gower-May model, the functional response that is used across more prey density regions is the one that determines where the Hopf bifurcation occurs and whether there exists a saddle bifurcation of limit cycles. However, the functional response used across only one prey density region causes a shift in the saddle bifurcation of limit cycles (if it exists) for every case, except surprisingly the IIH sub-model where there is no shift from the Ivlev behaviour (Figure S8). This result means that the sub-models with the Holling functional response describing two prey density regions result in only Holling model behaviour, but the sub-models with either the Ivlev or trigonometric functional response describing two prey density regions result in mixed dynamics.

Structural sensitivity is an intrinsic feature of a wide range of ecological models (Adamson and Morozov, 2013, 2014). In our work, we focus on predator-prey models and investigate the structural sensitivity that occurs with respect to the predation function. Our work has attempted to answer how to fit data in such a way to gain confidence in the results, regardless of which curve is used. We found that only the sub-models that use the functional responses fit around the coexistence steady state consistently decrease structural sensitivity across both the Rosenzweig-MacArthur and Leslie-Gower-May models. Using this result, we also give suggestions as to which intervals of prey density are most relevant when fitting a functional response to data and which intervals of prey density are most relevant for ecologists to collect data.

The main application of our work is in models that use data to fit the functional response. As we noted in Section 4.2, the choice of the fitted functional response is subjective if our method is used merely theoretically. However, in the case of applied models, our method is no longer subjective since we fit the functional responses to data instead of a subjectively chosen fixed functional response. Once all three functional responses are fit using our method to make their bifurcations structures more similar, the resulting bifurcation structure should show the behaviour of the data rather than a single chosen functional response.

Ecological data, however, is typically sparse as it is often difficult to collect (Ovaskainen and Soininen, 2011). Furthermore, the data is also typically very noisy. As a result, it is important that ecologists know what data is most important to collect to build accurate models that can be used with confidence to predict future population dynamics. In Section 4.2, we found that model predictions are most consistent if the functional responses match over the coexistence steady state prey density region. The implication is that better predictions can be obtained if there is additional data near the coexistence steady state. We believe that this region of prey densities is most important in fitting the functional responses because it corresponds to the steady state behaviour for values of the carrying capacity less than destabilizing bifurcation. The idea being that if the model behaviour is more similar before the destabilizing bifurcation, then the destabilizing bifurcation may occur for more similar values of the carrying capacity. This result opposes previous work that states saturation is the most important part to match of the functional responses (Aldebert and Stouffer, 2018).

To apply our result to ecological data and to make predictions in predator-prey models, we propose the following procedure:

1. 1.

   Collect data as per normal methods,

2. 2.

   Fit the three functional responses to this data and determine the parameter values for each $a_{i}$ and $b_{i}$,

3. 3.

   Use the functional responses to choose a fitting region around the stable or unstable coexistence steady state prey density values across all three models,

4. 4.

   Do one or both of the following:

   1. (a)

      Collect more data for prey densities in the coexistence steady state prey density region,

   2. (b)

      Give more weight to the data in the coexistence steady state prey density region,

5. 5.

   Fit the three functional responses again and use these functional responses to make predictions regarding the system.

This algorithm takes into account our result that collecting data for prey densities around the coexistence steady state appears to be most important while also considering the difficulty of collecting data (step 4b). Note that we have only considered structural sensitivity with respect to the functional response in this work. Models that are structurally sensitive with respect to other terms, such as prey growth, may require a different algorithm.