<figure><img src="image_1.png" /></figure>

HAL
open science

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

Clement Aldebert, Daniel B Stouffer

► To cite this version:

Clement Aldebert, Daniel B Stouffer. Community dynamics and sensitivity to model structure: toward a 1 probabilistic view of process-based model predictions 2. Journal of the Royal Society Interface, 2018, 15 (149), 10.1098/rsif.2018.0741. hal-01950903

HAL Id: hal-01950903

https://hal.science/hal-01950903

Submitted on 11 Dec 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Community dynamics and sensitivity to model structure: toward a
probabilistic view of process-based model predictions

Clement Aldebert¹,*, Daniel B Stouffer²

November 1, 2018

¹ Mediterranean Institute of Oceanography, Aix-Marseille University, Toulon University, CNRS/INSU, IRD, MIO, UM 110, 13288, Marseille, Cedex 09, France.

² School of Biological Sciences, University of Canterbury, Christchurch 8140, New Zealand.

* corresponding author: clement.aldebert@mio.osupytheas.fr

Keywords

model predictions ; Bayesian statistics ; structural sensitivity ; community dynamics ; predation
11

Abstract

12

13

14

15

16

17

18

19

20

21

22

23

24

25

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 toward 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 dif-
ferent streams of research have been growing, statistical inference and mechanistic modelling. While the former
27
aims to make predictions based on uncovering statistical relationships in large data sets, mechanistic modelling
28
30
31
32
33
34
35
36
37
38
39

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-fertilisation. Specifically, we use Bayesian statistics
to present probabilistic predictions of a deterministic mechanistic model—built around empirical data—in way
that takes into account uncertainty both in model construction and model parameterization. We therefore im-
prove 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 and MacArthur [5] predator-prey
model, which is known to be “structurally sensitive” [as defined by 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].
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61

Structural sensitivity is a common phenomenon that emerges in mechanistic biological models, which usually
aim to summarise multi-level processes into equations after adopting simple assumptions regarding the com-
plexity 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 and 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.
52

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 and MacArthur [5] model (sec-
tion 2). In particular, we explore the consequences of alternative predator functional responses in changing the
behaviours predicted. We achieve the shift from determinism (section 3) to probabilism (section 4) by merging
bifurcation theory for model analysis [24–26] and Bayesian statistics [27–29] by 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 data sets– indicates that parametric uncertainty is regularly overcome by
model uncertainty (section 5), an observation that has broad implications across various other challenges faced
61
## 2 A simple predator-prey model

We study the Rosenzweig-MacArthur model [5]:

$$ egin{cases} 
\frac{dN}{d\tau} = rN \left( 1 - \frac{N}{K} \right) - f(N)Y \\ 
\frac{dY}{d\tau} = (ef(N) - \mu)Y. 
\end{cases} 
\quad (1) $$

Here, the prey population $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 $\mu$, 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 \in F = \{f^{(H)}, f^{(I)}, f^{(t)}\}$, all of which only depend upon two parameters:

$$ f^{(H)}(N) = \frac{a^{(H)}N}{1 + a^{(H)}h^{(H)}N}, \quad f^{(I)}(N) = \frac{1}{h^{(I)}}(1 - e^{-a^{(I)}h^{(I)}N}), \quad f^{(t)}(N) = \frac{1}{h^{(t)}}\tanh(a^{(t)}h^{(t)}N), \quad (2) $$

which are the Holling Type II [30, 31], later denoted Holling], Ivlev [32], and the Hyperbolic Tangent [10], respectively. The parameters in each have the same mathematical meaning: $1/h^{(·)}$ corresponds to the maximum uptake rate and $a^{(·)}$ corresponds to the function's slope in the limit of no prey, but they originate from different perspectives. For Holling, $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 (1): $y = Y/r$, $t = r\tau$, $m = \mu/e$, $\varepsilon = e/r$, and write $x = N$. The re-scaled model reads:

$$ egin{cases} 
\frac{dx}{dt} = x \left( 1 - \frac{x}{K} \right) - f(x)y \\ 
\frac{dy}{dt} = \varepsilon (f(x) - m)y. 
\end{cases} 
\quad (3) $$

Upon re-scaling, the remaining parameters are (i) the prey carrying capacity $K$, (ii) the scaled predator mortality $m$, (iii) the time-scale factor $\varepsilon$, and (iv) the two parameters $a^{(·)}$ and $h^{(·)}$ 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
84

(e.g. environmental degradation, additional predator mortality due to harvesting).

85

# 3 Probabilistic predictions

86

## 3.1 Fitting functional responses on data

87

88

89

90

91

92

93

94

95

96

97

98

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 section 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, figure 1]. This consumer-resource system can be modeled 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 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.

99

100

101

102

103

104

105

106

107

108

109

110

111

112

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(\theta_f)$ that the set of functional-response parameters $\theta_f = (a^{(f)}, h^{(f)})$ of function $f$ predict observed data (figure 1a-c), assuming a Gaussian noise around $f$ [29, details in Appendix]. The probability $P(\theta_f)$ is proportional to the likelihood of $\theta_f$, and the set of probabilities corresponding to the inferred parameter values is called the posterior distribution in Bayesian statistics. This posterior distribution of $\theta_f$ is estimated by the HMCMC algorithm by learning from available data (figure 1a-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 1d-f).

113

114

115

116

117

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 \in 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 data set, the Hyperbolic Tangent has the highest probability of
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Results of the Markov Chain Monte Carlo estimate of parameters probability densities (posterior distributions) for each function. (a-c) Bivariate posterior probability distributions (grey levels) for each functional response parameters, i.e. the likelihood of these parameter values based on available data. Points '+' correspond to parameter values sampled to perform the probabilistic model analysis. Black bullets (one per panel) indicate the parameter values with the maximum likelihood, i.e. parameter values giving the best fit to data. (d-f) Functional responses fit to experimental data (points '+'). The parametric uncertainty from the HMCMC estimation gives a confidence interval (95%, shaded area) around the best-fitted functions (curves). Model uncertainty is derived through the relative likelihood that one function fits new data better than the others, knowing their respective parametric uncertainty.</figcaption></figure>
<table><caption>Table 1. Comparison of the three functional responses fit to data on *Gazella thomsoni* feeding on grass swards, based on WAIC (Widely-Applicable-Information-Criterion). WAIC estimates are given plus or minus standard error.</caption><thead><tr><th>functional response</th><th>WAIC</th><th>difference with best WAIC</th><th>weights P(f)</th><th>ranking</th></tr></thead><tbody><tr><td>Holling</td><td>-31.5 ± 9.2</td><td>4.0 ± 2.6</td><td>0.079</td><td>3</td></tr><tr><td>Ivlev</td><td>-30.5 ± 9.1</td><td>1.0 ± 0.9</td><td>0.350</td><td>2</td></tr><tr><td>Hyp. Tangent</td><td>-27.5 ± 9.1</td><td>0.0</td><td>0.572</td><td>1</td></tr></tbody></table>

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.

## 3.2 Overview of model predictions

Let use first present an overview of model predictions for a given functional response $f$ and its plausible parameter values $\theta_f$. We performed a bifurcation analysis based on analytical and numerical results (mathematical details in 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 \in [0, x_{max}]$, where $x_{max} \approx 197$ is the maximal prey abundance in functional response data, and $m \in [0, m_{max}]$, where $m_{max} := \sum_{f \in F} P(f) \int_{\mathbb{R}_+^2} P(\theta_f, D) f(x_{max}) d\theta_f \approx 5.44$ 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}$ restrain 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 since all models necessarily give identical predictions.

Figure 2a-c presents model predictions based on each functional response and including parametric uncertainty. The effect of parametric uncertainty will be introduced in the 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 $\theta_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 destabilised 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 time scale $\epsilon$, except for the bistability area that we computed numerically for $\epsilon = 1$; this area starts at lower (higher) K for higher (lower) $\epsilon$, which does not affect our conclusions.

As a general conclusion, using Holling and Ivlev functional responses lead 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 (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 $\alpha = (K, m, \epsilon)$ that do not relate to $f$, the probability that functional-response $f$ leads to the prediction $X$ is:

$$P(X|f) = \int_{\mathbb{R}_+^2} M(X, \alpha, f, \theta_f)P(\theta_f)d\theta_f, \quad (4)$$

where $M(X, \alpha, f, \theta_f) = 1$ if the model based on function $f$ and with parameter values $(\theta_f, \alpha)$ predicts dynamics $X$ and 0 otherwise, and the integral is over the range of the posterior distributions of the two parameters in $\theta_f$. In practice, this probability (4) is estimated by performing model analysis for a finite sample of $\theta_f$ values (here 1000 parameter sets) drawn randomly from the posterior distribution $P(\theta_f)$ [29]. Computing the probability (4) for different model parameters $\alpha$ 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 2a-c), showing how parametric uncertainty propagates into predictions.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Probabilistic predictions of Rosenzweig-MacArthur model, made with each alternative functional response (a-c) and averaged (d). Probabilistic predictions include both parametric uncertainty for each function (a-c), and model uncertainty for the average prediction (d). Each panel presents a probabilistic bifurcation diagram, where the colours indicate the qualitative system dynamics depending on the predator mortality rate and the prey carrying capacity: predator extinction (white), prey-predator coexistence at equilibrium (blue), prey-predator oscillations (green), and bistability with equilibrium coexistence or oscillations depending on initial population sizes (red). Colour gradients indicate the probability (Red-Green-Blue levels) of each model predictions. Thus, blurred areas (e.g. top left of panel (d)) indicate uncertain predictions. Calculations at point “+” are detailed in table 2.</figcaption></figure>
<table><caption>**Table 2.** Example to detail probabilities related to uncertainty, for each functional response (knowing *f*) and then averaged together (*f*-independent) based on WAIC weights from table 1. This example corresponds to $\alpha = (K^* = 42.6, m^* = 1.08)$, the “+” point in figures 2-4.</caption><thead><tr><th>description</th><th>probability</th><th>Holling</th><th>Ivlev</th><th>Hyp. Tangent</th><th>average</th></tr></thead><tbody><tr><td>function weight</td><td>P(f)</td><td>0.258</td><td>0.326</td><td>0.416</td><td></td></tr><tr><td rowspan="4">predicting qualitative dynamics</td><td>P("extinction"|f)</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>P("equilibrium"|f)</td><td>0.004</td><td>0.502</td><td>1.000</td><td>0.748</td></tr><tr><td>P("oscillations"|f)</td><td>0.996</td><td>0.498</td><td>0.000</td><td>0.252</td></tr><tr><td>P("bistability"|f)</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>source of predictive uncertainty:</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>- parametric uncertainty</td><td>U<sub>param</sub>(f)</td><td>0.001</td><td>0.175</td><td>0.000</td><td>0.061</td></tr><tr><td>- model uncertainty</td><td>U<sub>model</sub>(f)</td><td>0.373</td><td>0.162</td><td>0.126</td><td>0.158</td></tr><tr><td>- total</td><td>U<sub>tot</sub>(f)</td><td>0.373</td><td>0.337</td><td>0.126</td><td>0.219</td></tr><tr><td>sensitivity index</td><td>S(f)</td><td>0.997</td><td>-0.038</td><td>1.000</td><td>0.637</td></tr></tbody></table>

## 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$:

$$P(X) = \sum_{f \in F} P(f)P(X|f). \qquad (5)$$

The resulting probabilistic predictions are shown in figure 2d. 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 a red-green and a 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 $m^* = 1.08$ and $K^* = 42.6$ (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 $\theta_f$ for function $f$ drawn independently from their posterior distribution. We know from (4) 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 follows 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:

$$P_{\text{param}}(f) = \sum_{X \in Q} P(X|f) [1 - P(X|f)], \qquad (6)$$

where the sum is across possible qualitative dynamics. Similarly, we can define the probability that different dynamics are predicted by one parameter set $\theta_f$ of function f and one parameter set of any of the other alternative functions:

$$P_{\text{model}}(f) = \sum_{g \in F, g \neq f} \frac{P(g)}{1 - P(f)} \sum_{X \in Q} P(X|f) [1 - P(X|g)]. \qquad (7)$$

The sum over alternative functions g is weighted in order to take into account the respective weight of each alternative to function f. The probabilities (6-7) 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:

$$U_{\text{param}}(f) = P(f)P_{\text{param}}(f), \qquad (8)$$

which equals the parametric uncertainty that is proportional to the probability that f is the best function to describe new data, and:

$$U_{\text{model}}(f) = \frac{1 - P(f)}{|F| - 1} P_{\text{model}}(f), \qquad (9)$$

which equals the model uncertainty that is proportional to the probability that another function than 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 (8) is small in comparison to model uncertainty defined by (9), as the Holling function is very unlikely to be the best function (probability of 0.079) in comparison to the two alternative functions (figure 3a-b). According to this analysis, the total uncertainty associated to the choice of function f is:

$$U_{\text{tot}}(f) = U_{\text{param}}(f) + U_{\text{model}}(f), \qquad (10)$$

with the example for the Holling function shown in figure 3c.

To have a global overview, one can look at the average uncertainty across all the alternative functions that
<figure><img src="image_5.png" /></figure>

(c) Total uncertainty
with Holling

<figure><img src="image_6.png" /></figure>

always same
prediction

<figure><img src="image_7.png" /></figure>

(d) Average total uncertainty

<figure><img src="image_8.png" /></figure>

always different
predictions

<caption>**Figure 3.** Uncertainty in predictions made with Rosenzweig-MacArthur model. Details of the calculations if one chooses to use Holling functional response: uncertainty due to functional response's parameter values (a), to the choice of the functional response (b), and their sum (c). (d) Total uncertainty (like in (c)) averaged over the three alternative functional responses, weighted by their respective likelihood. Calculations at point "+" are detailed in table 2.</caption>
224 are considered (figure 3d):

$$ \bar{U}_{\text{tot}} = \sum_{f \in F} P(f)U_{\text{tot}}(f). \qquad (11) $$

225 To go deeper, we propose to quantify the respective importance of parametric and model uncertainty in the predictions made with a function f with:

226

$$ S(f) = 2 \left( \frac{U_{\text{model}}(f)}{U_{\text{tot}}(f)} - 0.5 \right). \qquad (12) $$

227 Note that this index can be computed only if there is uncertainty in predictions ($U_{\text{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_{\text{param}}(f) = U_{\text{model}}(f)$. 228 229 230 231 232 233 234

235 Based on the index (12) 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: 236

$$ \bar{S} = \sum_{f \in F} P(f)S(f). \qquad (13) $$

237 Figure 4a-c shows the index for the three alternative functional responses and our first data set. 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). 238 239 240 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 241 242 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 243 244 close to the predator maximum growth rate that itself is proportional to $1/h$. Thus, this parameter value has a strong impact on the predicted dynamics, explaining why parametric uncertainty is higher. As a representative 245 246 example, we show detailed calculations of these probabilities for ($m^* = 1.08, K^* = 42.6$) in table 2. 247 248
<figure><img src="image_9.png" /><figcaption>**Figure 4.** Source of uncertainty in predictions made with Rosenzweig-MacArthur model. Relative importance of parametric (negative value, blue) and model (positive value, red) uncertainty in the resulting total predictive uncertainty (grey area: total prediction uncertainty lower than 0.01). (a,b,c) Source of uncertainty if one of the alternative functional responses is chosen. (d) Average over the three alternative functions, weighted by their respective likelihood. Calculations at point "+" are detailed in table 2.</figcaption></figure>
<table><caption>Table 3. Comparison of the three functional responses fit to data on copepods (*Calanus pacificus*) feeding on diatoms (centric sp. and *Thalassiosira fluviatilis*), based on WAIC (Widely-Applicable-Information-Criterion). WAIC estimates are given plus or minus standard error.</caption><thead><tr><th>prey</th><th>functional response</th><th>WAIC</th><th>difference with best WAIC</th><th>weights</th><th>ranking</th></tr></thead><tbody><tr><td rowspan="3">centric sp.</td><td>Holling</td><td>-39.4 ± 7.7</td><td>12.4 ± 2.9</td><td>0.002</td><td>3</td></tr><tr><td>Ivlev</td><td>-48.2 ± 7.7</td><td>3.5 ± 1.2</td><td>0.145</td><td>2</td></tr><tr><td>Hyp. Tangent</td><td>-51.7 ± 8.2</td><td>0.0</td><td>0.853</td><td>1</td></tr><tr><td rowspan="3"><i>Thalassiosira fluviatilis</i></td><td>Holling</td><td>-22.6 ± 4.8</td><td>0.4 ± 2.2</td><td>0.299</td><td>3</td></tr><tr><td>Ivlev</td><td>-22.8 ± 4.5</td><td>2.0 ± 1.6</td><td>0.332</td><td>2</td></tr><tr><td>Hyp. Tangent</td><td>-23.0 ± 3.5</td><td>0.0</td><td>0.370</td><td>1</td></tr></tbody></table>

## 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, figure 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 5b). Also, the areas of highest predictive uncertainty in figure 5c are due to the parametric uncertainty of the dominant function, as indicated in blue in figure 5d. Conversely, predictive uncertainty is low in areas where it is mostly due to model uncertainty (in red in figure 5d). 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.

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, figure 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. There, 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
<figure><img src="image_10.png" /><figcaption>**Figure 5.** Example where one equation is almost certainly the best one among candidates. Overview of the analysis with data on copepods (*Calanus pacificus*) feeding on diatoms (centric sp.). All panels are drawn similarly as in earlier figures: (a) functional responses are fitted to data by a HMCMC algorithm; (b) average qualitative predictions (probabilistic bifurcation diagram); (c) average total uncertainty in predictions; (d) source of uncertainty in predictions.</figcaption></figure>
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>Figure 6. Example where all alternative equations are equally likely. Overview of the analysis with data on starved copepods (*Calanus pacificus*) feeding on diatoms (*Thalassiosira fluviatilis*). All panels are drawn similarly as in earlier figures: (a) functional responses are fitted to data by a HMCMC algorithm; (b) average qualitative predictions (probabilistic bifurcation diagram); (c) average total uncertainty in predictions; (d) source of uncertainty in predictions. Note that here we took $K_{max}$ as 5 times the maximum prey density in data, in order to show all the possible qualitative dynamics predicted by the model, without altering our conclusions.</figcaption></figure>
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 given 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 time-scale parameter $\varepsilon$ 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 6a), 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 1d-f,5a), 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 maximises 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 (upper rows), where all functions are equally likely based solely on the data. Conversely, Holling functional response might remain the worst model in the example in table 3 (lower 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 vs. 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 examples 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 optimised 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-fertilisation 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 to make probabilistic predictions 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 toward a better knowledge of the uncertainty and forecast horizon [as defined in 2] of model predictions in environmental sciences.

## Author Contributions

Original idea by CA and DBS. CA performed the research and wrote a first version of the article. Both authors contributed to discussion, wrote the final manuscript and contributed to revisions.

## 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.

## Funding Statement

CA received funding from European FEDER Fund under project 1166-39417. DBS received support from a Rutherford Discovery Fellowship and the Marsden Fund Council from New Zealand Government funding, both
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 CA.

Competing Interests

We have no competing interests.

References

[1] TR Anderson. Progress in marine ecosystem modelling and the "unreasonable effectiveness of mathematics". *Journal of Marine Systems*, 81:4-11, 2010.

[2] OL Petchey, M Pontarp, TM Massie, S Kéfi, A Ozgul, M Weilenmann, GM Palamara, F Altermatt, B Matthews, JM Levine, and et al. The ecological forecast horizon, and examples of its uses and determinants. *Ecology Letters*, 18:597-611, 2015.

[3] F Pennekamp, MW Adamson, OL Petchey, JC Poggiale, M Aguiar, BW Kooi, DB Botkin, and DL DeAngelis. The practice of prediction: What can ecologists learn from applied, ecology-related fields? *Ecological Complexity*, 32:156-167, 2017.

[4] RE Baker, JM Peña, J Jayamohan, and A Jérusalem. Mechanistic models versus machine learning, a fight worth fighting for the biological community? *Biology Letters*, 14:20170660, 2018.

[5] ML Rosenzweig and RH MacArthur. Graphical representation and stability conditions of predator-prey interaction. *The American Naturalist*, 97(895):209-223, 1963.

[6] F Cordoleani, D Nerini, M Gauduchon, A Morozov, and JC Poggiale. Structural sensitivity of biological models revisited. *Journal of Theoretical Biology*, 283:82-91, 2011.

[7] MR Myerscough, MJ Darwen, and WL Hogarth. Stability, persistence and structural stability in a classical predator-prey model. *Ecological Modelling*, 89:31-42, 1996.

[8] GF Fussmann and B Blasius. Community response to enrichment is highly sensitive to model structure. *Biology Letters*, 1(1):9-12, 2005.

[9] G Seo and GSK Wolkowicz. Sensitivity of the dynamics of the general rosenzweigmacarthur model to the mathematical form of the functional response: a bifurcation theory approach. *Journal of Mathematical Biology*, 76(7):1873-1906, 2018.
[10] AD Jassby and T Platt. Mathematical formulation of the relationship between photosynthesis and light for phytoplankton. *Limnology and Oceanography*, 21(4):540-547, 1976.

[11] SN Wood and MB Thomas. Super-sensitivity to structure in biological models. *Proceedings of the Royal Society of London B*, 266:565-570, 1999.

[12] JC Poggiale, M Baklouti, B Queguiner, and SALM Kooijman. How far details are important in ecosystem modelling: the case of multi-limiting nutrients in phytoplankton-zooplankton interactions. *Philosophical Transactions of the Royal Society of London B*, 365:3495-3507, 2010.

[13] T Gross, W Ebenhöh, and U Feudel. Enrichment and foodchain stability: the impact of different forms of predator-prey interaction. *Journal of Theoretical Biology*, 227:349-358, 2004.

[14] MW Adamson and AY Morozov. When can we trust our model predictions? unearthing structural sensitivity in biological systems. *Proceedings of the Royal Society of London A*, 469(2149), 2012.

[15] MW Adamson and AY Morozov. Bifurcation analysis of models with uncertain function specification: how should we proceed? *Bulletin of Mathematical Biology*, 76:1218-1240, 2014.

[16] MW Adamson and AY Morozov. Defining and detecting structural sensitivity in biological models: developing a new framework. *Journal of Mathematical Biology*, 69:1815-1848, 2014.

[17] C Aldebert, D Nerini, M Gauduchon, and JC Poggiale. Does structural sensitivity alter complexity-stability relationships? *Ecological Complexity*, 28:104-112, 2016.

[18] C Aldebert, D Nerini, M Gauduchon, and JC Poggiale. Structural sensitivity and resilience in a predator-prey model with density-dependent mortality. *Ecological Complexity*, 28:163-173, 2016.

[19] C Aldebert, BW Kooi, D Nerini, and JC Poggiale. Is structural sensitivity a problem of oversimplified biological models? insights from nested dynamic energy budget models. *Journal of Theoretical Biology*, 448:1-8, 2018.

[20] C Aldebert, BW Kooi, D Nerini, M Gauduchon, and JC Poggiale. Three-dimensional bifurcation analysis of a predator-prey model with uncertain formulation. *SIAM Journal on Applied Mathematics*, submitted revision.

[21] TR Anderson, WC Gentleman, and B Sinha. Influence of grazing formulations on the emergent properties of a complex ecosystem model in a global ocean general circulation model. *Progress in Oceanography*, 87: 201-213, 2010.

[22] TR Anderson, DO Hessen, A Mitra, DJ Mayor, and A Yool. Sensitivity of secondary production and export flux to choice of trophic transfer formulation in marine ecosystem models. *Journal of Marine Systems*, 125: 41-53, 2013.
439

[23] SM Vallina, BA Ward, S Dutkiewicz, and MJ Follows. Maximal feeding with active prey-switching: a kill-the-winner functional response and its effect on global diversity and biogeography. *Progress in Oceanography*, 120:93-109, 2014.

440

441

442

[24] J Guckenheimer and P Holmes. *Nonlinear oscillations, dynamical systems, and bifurcations of vector fields*. Springer, 1983. 459 pp.

443

444

[25] L Perko. *Differential equations and dynamical systems*. Springer, New York, 2nd edition, 1996. 525 pp.

445

[26] YuA Kuznetsov. *Elements of applied bifurcation theory*. Springer, New York, 3rd edition, 2004. 361 pp.

446

[27] G Claeskens and NL Hjort. *Model selection and model averaging*. Cambridge University Press, 2008. 311 pp.

447

448

[28] E Parent and E Rivot. *Introduction to Hierarchical Bayesian Modelling for ecological data*. CRC Press, Boca Raton, 2013. 405 pp.

449

450

[29] R McElreath. *Statistical rethinking: a Bayesian course with examples in R and Stan*. CRC Press, 2015. 483 pp.

451

452

[30] CS Holling. Some characteristics of simple types of predation and parasitism. *The Canadian Entomologist*, 91(7):385-398, 1959.

453

454

[31] CS Holling. The functional response of predators to prey density and its role in mimicry and population regulation. *Memoirs of the Entomological Society of Canada*, 45:3-60, 1965.

455

456

[32] VS Ivlev. *Experimental ecology of the feeding of fishes*. Pischepromizdat, Moscow, 1955. 302 pp. (translated from Russian by D. Scott, Yale University Press, New Haven, 1961).

457

458

[33] JF Wilmshurst, JM Fryxell, and PE Colucci. What constraints daily intake in thomson's gazelles. *Ecology*, 80(7):2338-2347, 1999.

459

460

[34] ML Rosenzweig. Paradox of enrichment: destabilization of exploitation ecosystems in ecological time. *Science*, 171(3969):385-387, 1971.

461

462

[35] M Scheffer, SR Carpenter, TM Lenton, J Bascompte, WA Brock, V Dakos, J van de Koppel, IA van de Leemput, SA Levin, EH van Nes, M Pascual, and J Vandermeer. Anticipating critical transitions. *Science*, 338:344-348, 2012.

463

464

465

[36] BW Frost. Effects of size and concentration of food particles on the feeding behaviour of the marine planktonic copepod *Calanus pacificus*. *Limnology and Oceanography*, 17(6):805-815, 1972.

466

467

[37] C Aldebert. *Uncertainty in predictive ecology: consequence of choices in model construction*. PhD thesis, 2016. Aix-Marseille University, 276 pp.

468
469

[38] BW Kooi and SALM Kooijman. The transient behaviour of food chains in chemostat. *Journal of Theoretical Biology*, 170:87–94, 1994.

470

471

[39] A Curtsdotter, HT Banks, JE Banks, M Jonsson, T Jonsson, AN Laubmeier, M Traugottand, and R Bommarco. Ecosystem function in predator-prey food webs–confronting dynamic models with empirical data. *Journal of Animal Ecology*, pages 1–15, 2018. doi: 10.1111/1365-2656.12892.

472

473

474

[40] AC Babtie, P Kirk, and MPH Stumpf. Topological sensitivity analysis for systems biology. *Proceedings of the National Academy of Sciences of the United States of America*, 111(52):18507–18512, 2014.

475