# Community dynamics and sensitivity to model structure: towards a probabilistic view of process-based model predictions

## Abstract

Statistical inference and mechanistic, process-based modelling represent two philosophically different streams of research whose primary goal is to make predictions. Here, we merge elements from both approaches to keep the theoretical power of process-based models while also considering their predictive uncertainty using Bayesian statistics. In environmental and biological sciences, the predictive uncertainty of process-based models is usually reduced to parametric uncertainty. Here, we propose a practical approach to tackle the added issue of structural sensitivity, the sensitivity of predictions to the choice between quantitatively close and biologically plausible models. In contrast to earlier studies that presented alternative predictions based on alternative models, we propose a probabilistic view of these predictions that include the uncertainty in model construction and the parametric uncertainty of each model. As a proof of concept, we apply this approach to a predator–prey system described by the classical Rosenzweig–MacArthur model, and we observe that parametric sensitivity is regularly overcome by structural sensitivity. In addition to tackling theoretical questions about model sensitivity, the proposed approach can also be extended to make probabilistic predictions based on more complex models in an operational context. Both perspectives represent important steps towards providing better model predictions in biology, and beyond.

### 1. Introduction

With the need for more accurate predictions in biology and environmental sciences [1–3], two philosophically different streams of research have been growing, statistical inference and mechanistic modelling. While the former aims to make predictions based on uncovering statistical relationships in large datasets, mechanistic modelling aims to make predictions based on causal mechanisms that explain observed patterns. In practice, the pros of one approach are the cons of the other, so a promising way forward would be to combine them in a ‘symbiotic relationship’ [4]. Here, we provide an example of such cross-fertilization. Specifically, we use Bayesian statistics to present probabilistic predictions of a deterministic mechanistic model—built around empirical data—in a way that takes into account uncertainty both in model construction and model parameterization. We, therefore, improve the model’s predictive capability by including prediction uncertainty while maintaining the explanatory power of a mechanistic model. As our focal model, we use the Rosenzweig & MacArthur [5] predator–prey model, which is known to be ‘structurally sensitive’ (as defined by Cordoleani *et al*. [6]); that is, apparently minor changes in model formulation can lead to a dramatic change in both quantitative (predicted biomasses over time) and qualitative predictions, such as prey–predator oscillations or the coexistence of alternative stable states [7–9].

Structural sensitivity is a common phenomenon that emerges in mechanistic biological models, which usually aim to summarize multi-level processes into equations after adopting simple assumptions regarding the complexity of the biological system of interest. As the entire complexity can rarely, if ever, be taken into account, available empirical data may be insufficient to statistically discriminate between alternative models [6,10]. Alternative models are known to make different predictions when the uncertain process is the infection in a host–pathogen system [11], the colimited uptake of nutrient [12], or predation in predator–prey and food-web models [6–8,13–20]. Some studies also indicate that major ocean-scale predictions, such as the dominance of phytoplankton groups, primary production and export to the deep ocean, can be deeply affected by this form of sensitivity [21–23]. In a way, structural sensitivity can be considered an extension of classical parameter sensitivity; moreover, Cordoleani *et al.* [6] and Adamson & Morozov [14] have provided examples where a change in model formulation has a higher effect on model predictions than an equivalent change in parameter values. However, as far as we know, such an analysis has never been statistically performed simultaneously across alternative models and all plausible parameter values.

Here, we ask whether the predictions made by a biological model are more sensitive to the uncertainty in its parameterization or to its mathematical formulation. To answer this question, we present a probabilistic view of predictions made by a deterministic predator–prey model, the Rosenzweig & MacArthur [5] model (§2). In particular, we explore the consequences of alternative predator functional responses in changing the behaviours predicted. We achieve the shift from determinism (§3) to probabilism (§4) by merging bifurcation theory for model analysis [24–26] and Bayesian statistics [27–29], thus harnessing the latter’s ability to include uncertainty in a model and to propagate it forward into predictions. Intriguingly, our probabilistic approach—applied on three example datasets—indicates that parametric uncertainty is regularly overcome by model uncertainty (§5), an observation that has broad implications across various other challenges faced in biology.

### 2. A simple predator–prey model

We study the Rosenzweig–MacArthur model [5]:

*N*follows logistic growth with

*per capita*growth rate

*r*and environmental carrying capacity

*K*in the absence of predators. The predator population

*Y*has a linear

*per capita*mortality rate

*μ*, conversion efficiency

*e*, and functional response

*f*(

*N*). All parameters and variables are strictly positive to ensure their biological meaning.

In line with the biology, a well-behaved functional response is expected to be strictly increasing, concave, saturating and vanishes only at 0. Here, we consider three possible alternatives for *f* ∈ *F* = {*f*^{(H)}, *f*^{(I)}, *f*^{(t)}}, all of which only depend upon two parameters:

*a*

^{(H)}and

*h*

^{(H)}are predator attack rate and handling time. For Ivlev, 1/

*h*

^{(I)}and

*a*

^{(I)}

*h*

^{(I)}are predator maximum digestion rate and satiation coefficient. The hyperbolic tangent is actually a phenomenological model without underlying biological assumptions, which nevertheless can at times provide a more accurate description of data than others [10].

To drop one parameter and ease the analysis without loss of information, we re-scale model (2.1): *y* = *Y*/*r*, *t* = *r*τ, *m* = *μ*/*e*, *ɛ* = *e*/*r*, and write *x* = *N*. The re-scaled model reads:

*K*, (ii) the scaled predator mortality

*m*, (iii) the timescale factor

*ɛ*and (iv) the two parameters and from the functional response that must be estimated from data. These first two are of biological interest as they might be affected by external factors (e.g. environmental degradation, additional predator mortality due to harvesting).

### 3. Probabilistic predictions

#### 3.1. Fitting functional responses on data

Empirical data were not considered in previous structural sensitivity analyses in Rozensweig–MacArthur model [8,9]. As our baseline, we therefore use three examples of functional-response data to provide an overview of our approach and its practical use. To introduce the approach, we first focus on one of these data sets (until §5). Specifically, we use data from an experiment where two individuals of *Gazella thomsoni* were eating hand-assembled grass swards, and their grazing rate as a function of grass biomass (i.e. a functional response) was estimated [33, fig. 1]. This consumer–resource system can be modelled with generic assumptions on population dynamics: (i) grass biomass follows a logistic growth limited for instance by space and nutrients availability; (ii) the growth rate of *Gazella thomsoni* population is proportional to its grass intake which is a strictly increasing, concave, saturating function of grass biomass that vanishes in absence of grass; (iii) *Gazella thomsoni* has a linear *per capita* mortality rate (e.g. ageing, harvesting). All these assumptions define the Rozensweig–MacArthur model, which is often refered as a predator–prey model but is generic enough to also describe a variety of consumer–resource systems.

Previous studies on structural sensitivity compared a deterministic analysis of model predictions based on best-fitted functions ([8,9,11,18,20], among others). Fitting a function to data implies uncertainty in the parameters being inferred, which means that there is additional insight to be gained by considering both parameter and model uncertainty. To perform a probabilistic analysis of model predictions, we first estimate the parametric uncertainty while fitting each functional response to data. This is achieved by a Hamiltonian Markov chain Monte Carlo (HMCMC) algorithm computing the probability *P*(** θ_{f}**) that the set of functional-response parameters

**= (**

*θ*_{f}*a*

^{(f)},

*h*

^{(f)}) of function

*f*predict observed data (figure 1

*a*–

*c*), assuming a Gaussian noise around

*f*([29], details in electronic supplementary material, appendix). The probability

*P*(

**) is proportional to the likelihood of**

*θ*_{f}**, and the set of probabilities corresponding to the inferred parameter values is called the posterior distribution in Bayesian statistics. This posterior distribution of**

*θ*_{f}**is estimated by the HMCMC algorithm by learning from available data (figure 1**

*θ*_{f}*a*–

*c*). The maximum of this distribution occurs at parameters having the maximum-likelihood and thus giving the best fit to data based on this criterion. Therefore, the HMCMC approach gives the same best-fitted functions as a maximum-likelihood estimation, but it adds information about parametric uncertainty for each function by considering the whole posterior probability distribution (figure 1

*d*–

*f*).

Based on the posterior distribution for each function, we use the widely applicable information criterion (WAIC) to derive Akaike weights which correspond to the relative predictive accuracy of each function—i.e. the probability *P*(*f*) that function *f* ∈ *F* gives the best fit to new data, conditional to the alternative functions that we consider [29]. This relative accuracy can be estimated using any other suitable method than WAIC, without altering our generic framework. For the first dataset, the hyperbolic tangent has the highest probability of 0.572 compared to 0.350 for Ivlev and 0.079 for Holling (exact values may change due to algorithm stochasticity, table 1). Based on that, one could decide to focus only on the hyperbolic tangent as it seems to be the best function. However, doing this neglects the fact that the alternative functions can be a better description of new data with a significant probability (0.428). Our approach hence aims to keep all the alternative functions, and to merge their predictions based on their respective likelihood through the use of model averaging.

functional response | WAIC | difference with best WAIC | weights P(f) |
ranking |
---|---|---|---|---|

Holling | −31.5 ± 9.2 | 4.0 ± 2.6 | 0.079 | 3 |

Ivlev | −30.5 ± 9.1 | 1.0 ± 0.9 | 0.350 | 2 |

hyperbolic tangent | −27.5 ± 9.1 | 0.0 | 0.572 | 1 |

#### 3.2. Overview of model predictions

Let us first present an overview of model predictions for a given functional response *f* and its plausible parameter values ** θ_{f}**. We performed a bifurcation analysis based on analytical and numerical results (mathematical details in electronic supplementary material, appendix) that reveals all the qualitative asymptotic dynamics predicted by the model for any environmental condition (

*K*,

*m*). We limited the range of these conditions to

*K*∈ ]0,

*x*

_{max}], where

*x*

_{max}≈ 197 is the maximal prey abundance in functional response data, and

*m*∈ ]0,

*m*

_{max}], where is the average mortality rate (among all functional responses) that allows predator survival for the considered range of prey abundance. The choice of

*m*

_{max}restrains our analysis to the range of parameter values that will give interesting results. Any

*m*>

*m*

_{max}will lead to predator extinction in most cases, and this is a trivial result that presents no interest for our study as all models necessarily give identical predictions.

Figure 2*a*–*c* presents model predictions based on each functional response and including parametric uncertainty. The effect of parametric uncertainty will be introduced in the next subsection and is shown by the blurred colours in the figure. For the moment, the reader can have an overview of model predictions for given parameter values ** θ_{f}** by looking at the areas with sharp colours. In this model, there is a trivial extinction equilibrium

*E*

_{(0)}= (0, 0) without any species, which can be reached only if the prey is initially absent. If both species are initially present, different asymptotic dynamics are predicted depending on the parameters and form of the functional response. An extinction equilibrium

*E*

_{(1)}= (

*K*, 0) without the predator always exists, but it is only stable in the white area in the figure. In other areas, prey carrying capacity is high enough to sustain the predator population and a coexistence equilibrium

*E*

_{(2)}= (

*x*

^{(2)},

*y*

^{(2)}) exists. It starts to exist for higher carrying capacity values when predator mortality (i.e. the loss that must be overcome to ensure predator survival) is higher. This equilibrium is stable only in the blue and red areas. In the green area found at high carrying capacity, the paradox of enrichment [34] has destabilized the coexistence equilibrium and there are stable prey–predator oscillations. Stable oscillations and the stable coexistence equilibrium are both possible in the red area, which is only predicted with the hyperbolic tangent at low predator mortality and high carrying capacity. In this bistability area, the model predicts that the system will converge to one of the two alternatives depending on its initial state. The presence of alternative stable states is of particular interest to study ecosystem resilience when facing external disturbances [18,35]. This other form of the paradox of enrichment has been found in predator–prey models incorporating a more detailed description of organisms biology [19]. All of the above qualitative predictions are independent of the timescale

*ɛ*, except for the bistability area that we computed numerically for

*ɛ*= 1; this area starts at lower (higher)

*K*for higher (lower)

*ɛ*, which does not affect our conclusions.

description | probability | Holling | Ivlev | hyp. tangent | average |
---|---|---|---|---|---|

function weight | P(f) |
0.258 | 0.326 | 0.416 | |

predicting qualitative dynamics | P(`extinction′|f) |
0.000 | 0.000 | 0.000 | 0.000 |

P(`equilibrium′|f) |
0.004 | 0.502 | 1.000 | 0.748 | |

P(`oscillations′|f) |
0.996 | 0.498 | 0.000 | 0.252 | |

P(`bistability′|f) |
0.000 | 0.000 | 0.000 | 0.000 | |

source of predictive uncertainty: | |||||

— parametric uncertainty | U_{param}(f) |
0.001 | 0.175 | 0.000 | 0.061 |

— model uncertainty | U_{model}(f) |
0.373 | 0.162 | 0.126 | 0.158 |

— total | U_{tot}(f) |
0.373 | 0.337 | 0.126 | 0.219 |

sensitivity index | S(f) |
0.997 | −0.038 | 1.000 | 0.637 |

As a general conclusion, using Holling and Ivlev functional responses leads to similar patterns of predictions and always one stable state whereas using the hyperbolic tangent allows us to predict bistability. Notably, the transitions between different qualitative predictions occur in different regions of parameter space for different models and correspond to the following biological phenomena: predator invasion (transcritical bifurcation: extinction – coexistence), onset of oscillations through enrichment paradox (supercritical Hopf bifurcation: coexistence – oscillations), and potential catastrophic shifts (i.e. tipping points) with a change in the number of alternative stable states (subcritical Hopf bifurcation: bistability – oscillations; and limit point of cycles bifurcation: bistability – coexistence). Note that if the number of alternative stable states is affected by structural sensitivity, it becomes hard to estimate ecosystem resilience under disturbances [18].

#### 3.3. Introducing parametric uncertainty into predictions

To introduce parametric uncertainty in our analysis, we now look at the probability that model (2.3) with function *f* predicts the different qualitative dynamics *X*, which can be any of ‘extinction’, ‘equilibrium’, ‘oscillations’, or ‘bistability’. For fixed values of the model parameters ** α** = (

*K*,

*m*, ε) that do not relate to

*f*, the probability that functional response

*f*leads to the prediction

*X*is:

*M*(

*X*,

**,**

*α**f*,

**) = 1 if the model based on function**

*θ*_{f}*f*and with parameter values (

**,**

*θ*_{f}**) predicts dynamics**

*α**X*and 0 otherwise, and the integral is over the range of the posterior distributions of the two parameters in

**. In practice, this probability (4) is estimated by performing model analysis for a finite sample of**

*θ*_{f}**values (here 1000 parameter sets) drawn randomly from the posterior distribution**

*θ*_{f}*P*(

**) [29]. Computing the probability (3.1) for different model parameters**

*θ*_{f}**therefore provides a probabilistic bifurcation analysis of model predictions based on one functional response. Translating the probabilities of each predictions into colour gradients gives the ‘fuzzy’ transitions between qualitative dynamics (figure 2**

*α**a*–

*c*), showing how parametric uncertainty propagates into predictions.

#### 3.4. Introducing model uncertainty into predictions

To go further, model uncertainty can be introduced through a weighted mean over the alternative functions. This allows us to compute the probability *P*(*X*) that the model generally predicts the different qualitative dynamics *X*:

*d*. Not surprisingly, the average prediction looks most similar to that made with the hyperbolic tangent since that model has a probability of 0.572. For low mortality values, the figure looks greener, with red-green and blue-green areas. In these areas, the hyperbolic tangent leads to bistability or a stable equilibrium, respectively, whereas the two other functions lead to oscillations and have a cumulative probability of 0.428. A representative example to illustrate the utility of these calculations occurs at and (point ‘+’ in figure 2), where all functional responses predict completely different qualitative dynamics: (i) the Holling model predicts oscillations, (ii) the hyperbolic tangent predicts equilibrium coexistence, and (iii) the Ivlev predicts both of these dynamics with roughly equal probability (table 2, upper half).

### 4. Identifying the sources of uncertainty

The probabilistic predictions presented in figure 2 present a form of predictive uncertainty, in the sense that predictions can change due to the uncertainty in the choice of a function and parametric uncertainty. As is commonly done in other studies, let us assume that we have chosen one function *f* among the possible ones. We will now explore the consequences of this choice, both in terms of parametric uncertainty and model uncertainty.

For the parametric uncertainty, consider two parameter sets ** θ_{f}** for function

*f*drawn independently from their posterior distribution. We know from (3.1) the probability

*P*(

*X*|

*f*) that one parameter set predicts the dynamics

*X*. Moreover, predicting dynamics

*X*or not is a binary outcome. Therefore, the joint probability that one set of parameters predicts the dynamics

*X*and the other does not follow a binomial distribution and is

*P*(

*X*|

*f*)[1 −

*P*(

*X*|

*f*)]. This probability is maximal (0.25) if there is an equal probability that

*f*predicts dynamics

*X*or not (

*P*(

*X*|

*f*) = 1 −

*P*(

*X*|

*f*) = 0.5), which corresponds to the highest uncertainty on predicting dynamics

*X*. Conversely, the probability is null if there is no uncertainty in predicting dynamics

*X*, i.e.

*P*(

*X*|

*f*) = 0 (

*f*never predicts

*X*) or

*P*(

*X*|

*f*) = 1 (

*f*always predicts

*X*). Thus, the total probability that these two parameter sets ever predict different dynamics is:

**of function**

*θ*_{f}*f*and one parameter set of any of the other alternative functions:

*g*is weighted in order to take into account the respective weight of each alternative to function

*f*. The probabilities (4.1–4.2) quantify the probabilities of making different predictions knowing that we are looking at the parametric or the model uncertainty. This implies that we can define:

*f*is the best function to describe new data, and:

*f*is the best function to describe new data. Note that this second probability is averaged over the |

*F*| − 1 functions other than

*f*, to give an equal weight to parametric and model uncertainty if all the alternative functions are equally plausible, independent of their number. In other words, we are comparing

*f*against one of its alternatives, not

*f*against all its alternatives, because in the latter case the approach would always indicate a high model uncertainty if the number of alternative functions is sufficiently high. If one chooses the Holling function in our example, the parameter uncertainty as defined in (4.3) is small in comparison to model uncertainty defined by (4.4), as the Holling function is very unlikely to be the best function (probability of 0.079) in comparison to the two alternative functions (figure 3

*a*,

*b*). According to this analysis, the total uncertainty associated to the choice of function

*f*is:

*c*.

To have a global overview, one can look at the average uncertainty across all the alternative functions that are considered (figure 3*d*):

*f*with:

*U*

_{tot}(

*f*) > 0). It is positive if model uncertainty is greater than parametric uncertainty and negative otherwise. In addition,

*S*(

*f*) = 1 if the function is not a suitable candidate (

*P*(

*f*) = 0) or its parametric uncertainty does not affect the predicted dynamics, which are affected (in any amount) by model uncertainty. Conversely,

*S*(

*f*) = −1 if predictions are affected by parametric uncertainty, and there is no model uncertainty i.e.

*P*(

*f*) = 1. Finally,

*S*(

*f*) = 0 can occur if (i) all alternative functions have an equal likelihood and they all lead to the same uncertain predictions (i.e.

*P*(

*X*|

*f*) =

*P*(

*X*|

*g*) for all dynamics

*X*and alternative functions

*g*); or more generally (ii) the function

*f*is more likely than others but others lead to predictions that are different enough to get

*U*

_{param}(

*f*) =

*U*

_{model}(

*f*).

Based on the index (4.7) computed for each function, one can look at its average value across the alternative functions to get a global overview of the source of uncertainty in the model predictions:

*a*–

*c*shows the index for the three alternative functional responses and our first dataset. If Holling is chosen, model uncertainty is always higher than the parametric uncertainty, as alternative equations are more plausible and lead to different predictions. Conversely, if the hyperbolic tangent is chosen, parametric uncertainty is higher on average than model uncertainty (

*S*(

*f*

^{(t)}) < 0 over 66.6% of the parameter space). This equation has the highest likelihood, which explains why its own parametric uncertainty appears to be more important. The Ivlev represents an intermediate case where model uncertainty is higher on average than parametric uncertainty (

*S*(

*f*

^{(I)}) > 0 over 69.0% of the parameter space). On average, choosing one equation creates a higher predictive uncertainty due to model uncertainty than parametric uncertainty in 66.3% of the parameter space. It is worth noting that parametric uncertainty is higher for high mortality rates, which are close to the predator maximum growth rate that itself is proportional to . Thus, this parameter value has a strong impact on the predicted dynamics, explaining why parametric uncertainty is higher. As a representative example, we show detailed calculations of these probabilities for in table 2.

### 5. Overview of the method’s possible outcomes and use

The underlying idea behind our approach is that the cost of choosing one model to make a prediction increases if (i) there is a high probability that another model also fits available data well and (ii) the two models make different predictions. Thus, this approach extends earlier studies on structural sensitivity ([8,9], among others) by considering the fact that, even if alternative models predict different dynamics, one model may outperform others in fitting available data. Again with the Rosenzweig–MacArthur model, we illustrate this possibility in figure 5 by using data on the ingestion rate (i.e. a functional response) of copepods (*Calanus pacificus*) feeding on diatoms (centric sp.) to parameterize the alternative functional responses ([36, fig. 4]). One functional response—the hyperbolic tangent—fits the data better than others, with a high probability (0.853) to be the best description of new data (table 3, upper rows). Thus, the predictions made with the hyperbolic tangent functional response strongly drive the average predictions across alternative functional responses (figure 5*b*). Also, the areas of highest predictive uncertainty in figure 5*c* are due to the parametric uncertainty of the dominant function, as indicated in blue in figure 5*d*. Conversely, predictive uncertainty is low in areas where it is mostly due to model uncertainty (in red in figure 5*d*). This example illustrates the idea that, if a change in equations leads to different predictions, it mostly matters whether or not alternative equations are actually likely to be chosen.

prey | functional response | WAIC | difference with best WAIC | weights | ranking |
---|---|---|---|---|---|

centric sp. | Holling | −39.4 ± 7.7 | 12.4 ± 2.9 | 0.002 | 3 |

Ivlev | −48.2 ± 7.7 | 3.5 ± 1.2 | 0.145 | 2 | |

Hyp. Tangent | −51.7 ± 8.2 | 0.0 | 0.853 | 1 | |

Thalassiosira fluviatilis |
Holling | −22.6 ± 4.8 | 0.4 ± 2.2 | 0.299 | 3 |

Ivlev | −22.8 ± 4.5 | 2.0 ± 1.6 | 0.332 | 2 | |

Hyp. Tangent | −23.0 ± 3.5 | 0.0 | 0.370 | 1 |

As a third and last example, we fit the functional responses to the data on the ingestion rate of starved copepods *Calanus pacificus* feeding on diatoms (*Thalassiosira fluviatilis*) [36, fig. 2]. Here, in contrast to the earlier examples, all three alternative functional responses have roughly the same likelihood (table 3, lower rows). As a result, the predictive uncertainty of the population model is mostly due to the model uncertainty (figure 6). This last example shares similarities to the one we used to introduce our approach. In figures 1–4, two alternative equations had a high likelihood in comparison to the third. Thus, one may imagine removing the most unlikely function from the analysis, and to keep only the equations giving an almost equally good fit. Doing so will end up either in one of two cases. First, something like figure 6 where some functions remain equally likely; or second, something like figure 5 where one of the alternative functions has a high likelihood in comparison to all others, and we can imagine only considering this best-fit function. In this last case, one may remove all functions except the one giving the best fit. By doing so, our approach simplifies in a parameter sensitivity analysis, as the scientist decides that the uncertainty in model construction—defined for the set of alternative equations that we choose to consider—can be neglected. Of course, deciding to neglect the uncertainty in model construction and choosing the best function (according to available data) can only be done after considering alternative models to fit to data, and having determined that one seems better than others.

### 6. Discussion

Here, we present probabilistic predictions of system dynamics in the classical Rosenzweig–MacArthur model that take into account both model and parametric uncertainty. Overall, we show that uncertainty in model formulation regularly leads to a larger predictive uncertainty than does usual parametric uncertainty. Moving forward, it is worth noting that uncertainty likely also involves other parameter values and processes. For predator mortality and prey carrying capacity, the resulting prediction uncertainty can be estimated directly from figure 3. Conversely, a known uncertainty based on data on the timescale parameter *ɛ* or the intrinsic dynamics of the prey (here specified as logistic growth) can be taken into account by also sampling their posterior distributions simultaneously with those of the functional response and its parameters. Knowing the respective contribution of different biological processes to the resulting predictive uncertainty would be a major step forward. To go this way, sampling a higher-dimensional parameter space is feasible with existing algorithms [29], but performing each model analysis might be computationally prohibitive. Nevertheless, the cost of conducting such an analysis may ultimately be lower than the cost of making the wrong biological/environmental decision due to an unreasonable faith in model predictions.

Decreasing predictive uncertainty can be achieved by decreasing the uncertainty in either model construction, model parameterization, or both. To do this, the naive idea of improving the experimental data collection might help to greatly decrease the model uncertainty. In our three examples, we are studying a process that is assumed to be a saturating function of prey abundance. However, if this saturation is not present in the collected data (figure 6*a*), it necessarily becomes harder to identify the best model. Indeed, the model parameter defining the plateau is less constrained by data (in comparison to figures 1*d*–*f*, 5*a*), increasing the parametric uncertainty of each alternative equations. Thus, alternative equations tend to have highly overlapping confidence intervals, subsequently increasing the uncertainty in model selection. Therefore, our study highlights the importance of designing experiments in a way that maximizes the constraints on parameter values of alternative models to decrease the uncertainty in selecting the ‘best’ model.

Finding the ‘best’ model in a given situation may imply arguments beyond the simple fit of models to data, such as the fact that one of the alternative equations is a well-established model in the literature that is based on valuable theoretical arguments (e.g. underlying mechanisms). Though we did not do so here, such non-quantitative arguments can be taken into account in our quantitative analysis thanks to the concept of prior probabilities. Prior probabilities are used in Bayesian statistics to give an *a priori* weight to some parameter values, or here alternative models, before estimating their likelihood from data [29]. With prior probabilities, one could thus give more weight to the well-established Holling functional response instead of the phenomenological hyperbolic tangent. As a consequence, this function might become the best model in the example in the table 3 (lower rows), where all functions are equally likely based solely on the data. Conversely, the Holling functional response might remain the worst model in the example in table 3 (upper rows), where the hyperbolic tangent is about 400-fold more likely according to data.

Our framework can also be easily extend to consider quantitative as well as qualitative differences in predictions (e.g. equilibrium versus oscillations). Quantitative differences are important for topics such as resource management, disease and pest control, or species invasions. Here, we combined Bayesian statistics with a qualitative analysis of model predictions (bifurcation analysis). Depending on model complexity and the question at stake, one can combine Bayesian statistics with quantitative predictions (e.g. numerical simulations) or include quantitative aspects in the bifurcation analysis. Quantitative aspects of structural sensitivity have already been considered in earlier studies [6,37]. However, including quantitative aspects in the bifurcation analysis might become computationally prohibitive. Thus, getting the quantitative predictions of the model of interest is the challenge in extending our framework, which is well-suited to consider uncertainty in both quantitative and qualitative predictions.

The predictions made with the predator–prey model we used as an example could not be compared to data on the temporal dynamics of the studied system. Indeed, we extracted functional-response data from studies on grazing/ingestion rates of a few predator individuals—functional response *sensu stricto*—but these studies did not follow the temporal dynamics of a prey–predator system over many generations (i.e. the scale of the population model). Conversely, some experiments on population dynamics are not combined with functional response experiment, and functional-response parameters are optimized so that the predicted system dynamics fit the temporal data (e.g. [38,39]). However, if one has access to both types of data, the additional information on temporal dynamics can be used to constrain the probabilistic analysis. One way to do so is to remove candidate models that never predict the observed qualitative dynamics. Another complementary way is to perform the parameter estimation with constraints coming from fitting both the functional-response data and data on the temporal dynamics. This approach might be the best way to solve the issue of structural sensitivity. However, it is worth noting that this is a truly idealistic case. Indeed, for organisms with a long lifetime (months, years), collecting temporal data on their population dynamics would require a long-term monitoring (at least many years). Therefore, making probabilistic predictions prior to such an experiment would still be of interest, given the time needed to acquire data.

Though a similar idea has been used to improve parameter inference [40], as far as we know this is the first time the full approach proposed here has been used. An intriguing application of our probabilistic approach would be to conduct structural-sensitivity analysis for additional data sets corresponding to other prey–predator species. This might allow the classification of organisms with population dynamics that are inherently more or less predictable, either because of parameter values or the uncertainty around their functional-response data. This approach can also be extended to different models, for example, to test the recent hypothesis that mass-balanced prey–predator models with maintenance are less structurally sensitive [19]. More generally, our proposed approach may benefit from further cross-fertilization with the approach of partially specified models [14–16]. That approach provides generic results but cannot currently make probabilistic predictions of the sort we make here. Another way to include uncertainty would be to use stochastic differential equations where the uncertain process is randomly drawn to include all data variability into model predictions.

To conclude, we have conducted a proof-of-concept study outlining a novel approach that considers both parametric uncertainty and uncertainty in biological model construction while presenting model prediction. We achieved this by bringing elements of Bayesian statistics into the analysis of deterministic dynamical systems. Here, the prey–predator system considered is small enough to get analytical results on bifurcations, as this provides a comprehensive overview of qualitative model predictions. More complex models can also be studied by the proposed approach, by adapting the automatic model analysis to model complexity. For instance, a full overview of qualitative model predictions can be obtained if a numerical bifurcation analysis can be performed, or a sample of numerical simulations of different scenarios can be used for models as large as global ecosystem models. Some of these models are known to make different large-scale predictions such as the global dominance of phytoplankton groups or the ocean primary production [21–23]. Those results were based on comparisons between alternative models, and the incorporation of our approach would allow probabilistic predictions to be made based on different plausible biological models. Indeed, these small changes in model construction affect predictions at the ocean-scale, but also the coexistence of alternative stable states and the predicted resilience at the scale of a few populations in interactions. Therefore, extending our approach to complex operational models, together with theoretical analyses to rank the processes and species according to the level of predictive uncertainty they produced at the ecosystem level, would be critical advances towards a better knowledge of the uncertainty and forecast horizon (as defined in [2]) of model predictions in environmental sciences.

### Data accessibility

Experimental data can be found in the references cited in the manuscript. The codes used to perform the analysis are available upon request to C.A.

### Authors' contributions

Original idea by C.A. and D.B.S. C.A. performed the research and wrote a first version of the article. Both authors contributed to discussion, wrote the final manuscript and contributed to revisions.

### Competing interests

We have no competing interests.

### Funding

C.A. received funding from European FEDER Fund under project 1166-39417. D.B.S. received support from a Rutherford Discovery Fellowship and the Marsden Fund Council from New Zealand Government funding, both of which are managed by the Royal Society Te Apārangi (RDF-13-UOC-003 and 16-UOC-008).

## Acknowledgements

We acknowledge David Nerini, Jean-Christophe Poggiale, Mathias Gauduchon, Claude Manté, Gerald Gregori and Andrew Letten for stimulating discussions. We also thank the anonymous referees for constructive comments that helped to improve the manuscript.

### References

- 1.
Anderson TR . 2010 Progress in marine ecosystem modelling and the ‘unreasonable effectiveness of mathematics’.**J. Mar. Sys.**, 4-11. (doi:10.1016/j.jmarsys.2009.12.015) Crossref, Web of Science, Google Scholar**81** - 2.
Petchey OL *et al.*2015 The ecological forecast horizon, and examples of its uses and determinants.**Ecol. Lett.**, 597-611. (doi:10.1111/ele.12443) Crossref, PubMed, Web of Science, Google Scholar**18** - 3.
Pennekamp F, Adamson MW, Petchey OL, Poggiale JC, Aguiar M, Kooi BW, Botkin DB, DeAngelis DL . 2017 The practice of prediction: what can ecologists learn from applied, ecology-related fields?**Ecol. Complexity**, 156-167. (doi:10.1016/j.ecocom.2016.12.005) Crossref, Web of Science, Google Scholar**32** - 4.
Baker RE, Peña JM, Jayamohan J, Jérusalem A . 2018 Mechanistic models versus machine learning, a fight worth fighting for the biological community?**Biol. Lett.**, 20170660. (doi:10.1098/rsbl.2017.0660) Link, Web of Science, Google Scholar**14** - 5.
Rosenzweig ML, MacArthur RH . 1963 Graphical representation and stability conditions of predator–prey interaction.**Am. Nat.**, 209-223. (doi:10.1086/282272) Crossref, Web of Science, Google Scholar**97** - 6.
Cordoleani F, Nerini D, Gauduchon M, Morozov A, Poggiale JC . 2011 Structural sensitivity of biological models revisited.**J. Theor. Biol.**, 82-91. (doi:10.1016/j.jtbi.2011.05.021) Crossref, PubMed, Web of Science, Google Scholar**283** - 7.
Myerscough MR, Darwen MJ, Hogarth WL . 1996 Stability, persistence and structural stability in a classical predator–prey model.**Ecol. Modell.**, 31-42. (doi:10.1016/0304-3800(95)00117-4) Crossref, Web of Science, Google Scholar**89** - 8.
Fussmann GF, Blasius B . 2005 Community response to enrichment is highly sensitive to model structure.**Biol. Lett.**, 9-12. (doi:10.1098/rsbl.2004.0246) Link, Web of Science, Google Scholar**1** - 9.
Seo G, Wolkowicz GSK . 2018 Sensitivity of the dynamics of the general Rosenzweig–MacArthur model to the mathematical form of the functional response: a bifurcation theory approach.**J. Math. Biol.**, 1873-1906. (doi:10.1007/s00285-017-1201-y) Crossref, PubMed, Web of Science, Google Scholar**76** - 10.
Jassby AD, Platt T . 1976 Mathematical formulation of the relationship between photosynthesis and light for phytoplankton.**Limnol. Oceanogr.**, 540-547. (doi:10.4319/lo.1976.21.4.0540) Crossref, Web of Science, Google Scholar**21** - 11.
Wood SN, Thomas MB . 1999 Super-sensitivity to structure in biological models.**Proc. R. Soc. Lond. B**, 565-570. (doi:10.1098/rspb.1999.0673) Link, Web of Science, Google Scholar**266** - 12.
Poggiale JC, Baklouti M, Queguiner B, Kooijman SALM . 2010 How far details are important in ecosystem modelling: the case of multi-limiting nutrients in phytoplankton–zooplankton interactions.**Phil. Trans. R. Soc. B**, 3495-3507. (doi:10.1098/rstb.2010.0165) Link, Web of Science, Google Scholar**365** - 13.
Gross T, Ebenhöh W, Feudel U . 2004 Enrichment and foodchain stability: the impact of different forms of predator–prey interaction.**J. Theor. Biol.**, 349-358. (doi:10.1016/j.jtbi.2003.09.020) Crossref, PubMed, Web of Science, Google Scholar**227** - 14.
Adamson MW, Morozov AY . 2012 When can we trust our model predictions? Unearthing structural sensitivity in biological systems.**Proc. R. Soc. A**, 20120500. (doi:10.1098/rspa.2012.0500) Link, Google Scholar**469** - 15.
Adamson MW, Morozov AY . 2014 Bifurcation analysis of models with uncertain function specification: how should we proceed?**Bull. Math. Biol.**, 1218-1240. (doi:10.1007/s11538-014-9951-9) Crossref, PubMed, Web of Science, Google Scholar**76** - 16.
Adamson MW, Morozov AY . 2014 Defining and detecting structural sensitivity in biological models: developing a new framework.**J. Math. Biol.**, 1815-1848. (doi:10.1007/s00285-014-0753-3) Crossref, PubMed, Web of Science, Google Scholar**69** - 17.
Aldebert C, Nerini D, Gauduchon M, Poggiale JC . 2016 Does structural sensitivity alter complexity–stability relationships?**Ecol. Complex.**, 104-112. (doi:10.1016/j.ecocom.2016.07.004) Crossref, Web of Science, Google Scholar**28** - 18.
Aldebert C, Nerini D, Gauduchon M, Poggiale JC . 2016 Structural sensitivity and resilience in a predator–prey model with density-dependent mortality.**Ecol. Complex.**, 163-173. (doi:10.1016/j.ecocom.2016.05.004) Crossref, Web of Science, Google Scholar**28** - 19.
Aldebert C, Kooi BW, Nerini D, Poggiale JC . 2018 Is structural sensitivity a problem of oversimplified biological models? Insights from nested dynamic energy budget models.**J. Theor. Biol.**, 1-8. (doi:10.1016/j.jtbi.2018.03.019) Crossref, PubMed, Web of Science, Google Scholar**448** - 20.
Aldebert C, Kooi BW, Nerini D, Gauduchon M, Poggiale JC . Submitted. Three-dimensional bifurcation analysis of a predator–prey model with uncertain formulation.**SIAM J. Appl. Math.**Web of Science, Google Scholar - 21.
Anderson TR, Gentleman WC, Sinha B . 2010 Influence of grazing formulations on the emergent properties of a complex ecosystem model in a global ocean general circulation model.**Prog. Oceanogr.**, 201-213. (doi:10.1016/j.pocean.2010.06.003) Crossref, Web of Science, Google Scholar**87** - 22.
Anderson TR, Hessen DO, Mitra A, Mayor DJ, Yool A . 2013 Sensitivity of secondary production and export flux to choice of trophic transfer formulation in marine ecosystem models.**J. Mar. Syst.**, 41-53. (doi:10.1016/j.jmarsys.2012.09.008) Crossref, Web of Science, Google Scholar**125** - 23.
Vallina SM, Ward BA, Dutkiewicz S, Follows MJ . 2014 Maximal feeding with active prey-switching: a kill-the-winner functional response and its effect on global diversity and biogeography.**Prog. Oceanogr.**, 93-109. (doi:10.1016/j.pocean.2013.08.001) Crossref, Web of Science, Google Scholar**120** - 24.
Guckenheimer J, Holmes P . 1983**Nonlinear oscillations, dynamical systems, and bifurcations of vector fields**. Berlin, Germany: Springer. Crossref, Google Scholar - 25.
Perko L . 1996**Differential equations and dynamical systems**, 2nd edn. New York, NY: Springer. Crossref, Google Scholar - 26.
Kuznetsov YuA . 2004**Elements of applied bifurcation theory**, 3rd edn. New York, NY: Springer. Crossref, Google Scholar - 27.
Claeskens G, Hjort NL . 2008**Model selection and model averaging**. Cambridge, UK: Cambridge University Press. Crossref, Google Scholar - 28.
Parent E, Rivot E . 2013**Introduction to hierarchical Bayesian modelling for ecological data**. Boca Raton, FL: CRC Press. Google Scholar - 29.
McElreath R . 2015**Statistical rethinking: a Bayesian course with examples in R and Stan**. Boca Raton, FL: CRC Press. Crossref, Google Scholar - 30.
Holling CS . 1959 Some characteristics of simple types of predation and parasitism.**Can. Entomol.**, 385-398. (doi:10.4039/Ent91385-7) Crossref, Google Scholar**91** - 31.
Holling CS . 1965 The functional response of predators to prey density and its role in mimicry and population regulation.**Mem. Entomol. Soc. Canada**, 3-60. (doi:10.4039/entm9745fv) Google Scholar**45** - 32.
Ivlev VS . 1955**Experimental ecology of the feeding of fishes**. Moscow: Pischepromizdat. (translated from Russian by D. Scott, Yale University Press, New Haven, 1961). Google Scholar - 33.
Wilmshurst JF, Fryxell JM, Colucci PE . 1999 What constraints daily intake in Thomson’s gazelles.**Ecology**, 2338-2347. (doi:10.2307/176914) Web of Science, Google Scholar**80** - 34.
Rosenzweig ML . 1971 Paradox of enrichment: destablilization of exploitation ecosystems in ecological time.**Science**, 385-387. (doi:10.1126/science.171.3969.385) Crossref, PubMed, Web of Science, Google Scholar**171** - 35.
Scheffer M *et al.*2012 Anticipating critical transitions.**Science**, 344-348. (doi:10.1126/science.1225244) Crossref, PubMed, Web of Science, Google Scholar**338** - 36.
Frost BW . 1972 Effects of size and concentration of food particles on the feeding behaviour of the marine planktonic copepod*Calanus pacificus*.**Limnol. Oceanogr.**, 805-815. (doi:10.4319/lo.1972.17.6.0805) Crossref, Web of Science, Google Scholar**17** - 37.
Aldebert C . 2016 Uncertainty in predictive ecology: consequence of choices in model construction. PhD thesis, Aix-Marseille University, Marseille, France. Google Scholar - 38.
Kooi BW, Kooijman SALM . 1994 The transient behaviour of food chains in chemostat.**J. Theor. Biol.**, 87-94. (doi:10.1006/jtbi.1994.1170) Crossref, Web of Science, Google Scholar**170** - 39.
Curtsdotter A, Banks HT, Banks JE, Jonsson M, Jonsson T, Laubmeier AN, Traugottand M, Bommarco R . 2018 Ecosystem function in predator–prey food webs—confronting dynamic models with empirical data.**J. Anim. Ecol.**1-15. (doi:10.1111/1365-2656.12892) PubMed, Web of Science, Google Scholar - 40.
Babtie AC, Kirk P, Stumpf MPH . 2014 Topological sensitivity analysis for systems biology.**Proc. Natl Acad. Sci. USA**, 18 507-18 512. (doi:10.1073/pnas.1414026112) Crossref, Web of Science, Google Scholar**111**