The effects of climatic fluctuations and extreme events on running water ecosystems

Most research on the effects of environmental change in freshwaters has focused on incremental changes in average conditions, rather than fluctuations or extreme events such as heatwaves, cold snaps, droughts, floods or wildfires, which may have even more profound consequences. Such events are commonly predicted to increase in frequency, intensity and duration with global climate change, with many systems being exposed to conditions with no recent historical precedent. We propose a mechanistic framework for predicting potential impacts of environmental fluctuations on running-water ecosystems by scaling up effects of fluctuations from individuals to entire ecosystems. This framework requires integration of four key components: effects of the environment on individual metabolism, metabolic and biomechanical constraints on fluctuating species interactions, assembly dynamics of local food webs, and mapping the dynamics of the meta-community onto ecosystem function. We illustrate the framework by developing a mathematical model of environmental fluctuations on dynamically assembling food webs. We highlight (currently limited) empirical evidence for emerging insights and theoretical predictions. For example, widely supported predictions about the effects of environmental fluctuations are: high vulnerability of species with high per capita metabolic demands such as large-bodied ones at the top of food webs; simplification of food web network structure and impaired energetic transfer efficiency; and reduced resilience and top-down relative to bottom-up regulation of food web and ecosystem processes. We conclude by identifying key questions and challenges that need to be addressed to develop more accurate and predictive bio-assessments of the effects of fluctuations, and implications of fluctuations for management practices in an increasingly uncertain world.

where ε is a random variable with mean = 0 and variance σ [5,6]. This introduces density-dependent environmental stochasticity in population sizes -larger populations experience a proportionally larger change in density. This would be the case if the fluctuation affected a fundamental rate such as rmax, such as through a temperature-driven impact on physiology [7,8]. The model can easily be extended to directly include fluctuations in temperature-dependent life history [8,9] or interaction [10,11] parameters. These fluctuations may be correlated between populations if populations tend to have a similar thermal sensitivity, for example, which we model using a multivariate Gaussian process having a characteristic correlation between each population's perturbations, but uncorrelated with respect to time (no temporal auto-correlation). In other words, all populations may experience "good" and "bad" in a perfectly correlated manner (correlation = 1), or completely independently (correation = 0), but for any two time steps t and k, ε t is independent of ε k . A low correlation implies a poor high physiological mismatch between species.
Extreme events: Populations may also experience catastrophic declines due to sudden changes in conditions (extreme events), such as the advent of droughts or temperature spikes. The impact of such events can be both density-independent as well as dependent. For example, in the case of temperature change, metabolic stress would lead to density-independent mortality, whereas concurrent reduction in trophic resources would increase competition in larger populations, resulting in density dependent mortality. For simplicity, here we model only the density-dependent scenario by adding a probabilistic catastrophe term as an instantaneous reduction by a fraction δ in all populations (Lande 1993 Where , is given by the right hand size of eqn (S3), and pδ is the per-time step probability of an extreme event. This allows us to model extreme events occurring at different characteristic frequencies.
We now include metabolic constraints into the food web model (eqn (S3)) as follows. For simplicity, we will focus only on the size (and not temperature) components of these metabolic constraints, leaving the stochastic fluctuation and extreme events parameters in eqns (S3) & (S4) to simulate, albeit indirectly, the effect of thermal fluctuations.

Size-scaling parameterizations
For intrinsic growth rate we use where r0 is a normalization constant that includes the effect of temperature (main text eqn 2.1), m is the species' average adult body mass, and β = 0.75. These relationships are empirically well supported [8,9,12,13], and have been used previously in similar contexts [1,[14][15][16][17]. For the search rate coefficient between the i th resource and j th consumer species we use where, a 0 is a normalization constant, and φ ij ∈ [0, 1] is a dimensionless function that embodies attack success probability. Thus from eqn (S2), Note that in eqns (S6)-(S7), the commonly used quarter-power exponent that has been used for the scaling of search rate, as in numerous previous studies [15][16][17][18]. According to the results of Pawar et al [1,11,19,20], this is approximately the scaling exponent for 2D (two spatial dimensions; e.g., benthic) interactions only. For simplicity, we am thus assuming here that all food webs have 2D interactions only [e.g., 21]. Future work should consider the effect of 3D or a mixture of 2D-3D interactions (main text Fig.  1). Next, we assume that attack success probability is unimodal ( [22]) with respect to average resource mass (consumer-resource body mass ratio, or size-ratio): Here s determines how rapidly the function reaches its peak k, the size-ratio of maximum consumption rate. This is a simplification of the potentially complex dynamics of consumer-resource encounter and consumption rates in nature, but captures an important feature of empirically observed attack and capture rates -these are unimodal because attack success decreases and handling time increases at extreme sizeratios [23][24][25]. The actual values of k and s are expected to vary with type of consumption (foraging) strategy as well as habitat type. For example, in the case of predator-prey interactions, because smaller organisms have greater mass-specific power relative to larger ones, and can hence handle a larger range of prey sizes, k may be closer 1 (or even <1) for small consumers, and s smaller (a more gradual decline of consumption rate at extreme ratios). Brose et al. [26] have shown that invertebrate consumers do indeed have a k closer to 1 than vertebrates across disparate habitat types, suggesting their superior ability to attack and capture prey closer to their own size. Our results are qualitatively insensitive to a wide range of choice of parameters k and s (Table S2).

Stochastic community assembly under environmental fluctuations
We simulate open, dynamically assembling food webs till the system reaches immigration-extinction equilibrium (IEE) [22,27] as follows:  Immigration. Beginning with the establishment at least one basal species, with some per-time-step probability pc, a species' population is introduced at an extinction threshold biomass abundance xe. Each immigrant species is generated by sampling a body size from a Beta(1,ω) distribution that captures the qualitative properties of empirical body size distributions [22]. Inter-and intra-specific interaction parameters are determined by body sizes (eqns (S5)-(S8)).  Trophic linking. Upon colonization, the j th immigrant establishes a trophic link to the i th pre-existing one with a connectance probability p ij . For each assembly simulation, conditional upon p ij , a "vulnerability probability" pv ranging between 0.5-1 is set: pv = 0.5 means that the j th immigrant is equally likely to be a resource or a consumer of the i th resident species (provided it was not basal), while pv = 1 meant that the j th immigrant can only be a consumer.  Interaction driven extinction. At every time step, a species is deleted from the system if its density drops below x e . This simulation is continued till the system reaches IEE. Simulations were performed in Python with scipy. Simulation parameters were chosen as follows [22,28] (Table S2):  e was fixed at 0.5 for all species, the approximate midpoint of the range reported from empirical data [9,13]. Values ranging from 0.1-1 don't change the simulation results qualitatively.  p v was set to 0.9 because this is the midpoint of the range [0.75-1] that yields communities with structural and dynamical characteristics similar to that of real ones, similar to the niche model [28]. The results do not change qualitatively over the full range of [0.5-1].  xe was set to 10 -20 ; the results do not change qualitative for values ranging from 10 -32 -10 -3 .
In addition, body size related simulation parameters were chosen as follows:  ω, which determines the shape of the immigrant pool size distribution was varied between 1 (uniform distribution; immigration rate independent of body size) and 2 (power law-like with slope = -2; immigration probability decreases with body mass). This range of ω was chosen because: (i) Empirical data show the distributions of sizes at large spatial scales are right-skewed [29,30], probably partly driven by speciation rate, which appears to follow a negative power law relationship with body size [31][32][33], possibly linked to metabolic scaling [34][35][36] and, (ii) dispersal ability should increase with body size [13,37]. Only considering (i) means that ω should be > 1; choosing ω = 2 sets a reasonably high upper limit to this bias. Considering (ii) means that the effect of (i) may be somewhat negated due to dispersal ability. However, because speciation within the local community also adds to the effective bias towards immigration by smaller species and because data on dispersal ability itself is biased towards larger organisms [13], it is unlikely that (ii) can overwhelm the effects of (i). Hence ω = 1, which yields the uniform distribution expected if the effects of (i) and (ii) exactly cancel out, is a reasonable lower limit to the immigration rate bias [22].  The log-body mass range [ , ] was chosen to be [-12,12] because this is approximately the range of species' log-body masses observed across empirical communities [26] (also see [22]).  k was chosen to vary randomly between 10 -3 and 10 3 with uniform probability (consumers 1000 times smaller to 1000 times larger than resources), which covers most of the range considered to be "optimal" (in the sense of viable, or evolutionarily stable strategy for the consumer) in previous studies on consumer-resource interactions [14,38], and accommodates potential differences in k across the most common trophic interaction types seen in food webs (i.e., predator-prey, herbivoreplant and parasitoid-host) [26].  The parameter s was set to 0.1; however, varying it between a wide range (0.05-0.5) does not change the results (results not shown; also see [22]).  The allometric constant r 0 was chosen to be 1.  The search-rate scaling constant a 0 was varied between 10 -3 -1 based upon recent empirical results [1] -the results shown here are for a 0 = 1. Other values in this range do not change the results qualitatively.  a ii was chosen according to the target mean n at IEE; larger values give larger feasible communities [6,21,22,39]. Table S1; Fig S1-4). For stability properties, we calculated the Jacobian matrix J, which is the S x S matrix of partial derivatives,

Food web properties. At IEE, we calculated various metrics, including stability properties (Main text Fig. 2 &
The local stability criterion of the system is that the real parts of all the S eigenvalues, λi = 1,…,S of C lie in the left half of the plane of complex numbers, i.e., max [λi = 1,…,S] = λmax < 0. We also calculated the system's "dynamic dimensionality", D -the number of eigenvalues required to account for 95% of the total sum of the eigenvalues. This is a measure of the number of dimensions that determine the system's return time and hence is indicative of the amount of dynamical redundancy, which we define as R = D -1 , i.e., greater the dimensionality, the lower its dynamical redundancy.

Food web characteristics:
We now identify and describe a set of topological and non-topological features of networks that were used to characterize emergent and evolving food webs.

Distribution of rmax's and body sizes:
The rmax of each population was calculated given its size using eqn (S5). We calculated two summary statistics about both size and growth rate distributions in the emergent food webs: the mean and the skewness.
Degree and degree distributions: We measured the degree of the directed, weighted food-webs by converting the corresponding adjacency matrix A to an unweighted, binary (but directed) one A UB with its n(n -1) off-diagonal elements given by, The degree of each node was then calculated from A UB as a sum of its in-and out-degree. From this, the average degree (k av ), was calculated.
Connectance: We measure connectance as, where and are directed edges between the i th and the j th nodes.
Trophic levels: Using the unweighted version of the interaction matrix, we calculated the shortest path length Gij between each i-j vertex pair ( [40]). The average shortest path length, Gav was then calculated as the average of all the G ij 's. Here because A UB represents a directed graph, we calculated G av as, ∑ , (S11 )  Table S1: Summary of qualitative theoretical results of the food web assembly modelling, comparison with previous theoretical work, and empirical support (for or against). An asterisk (*) indicates that the previous theoretical or empirical study supports or partly supports our results. See main text for references, except for the following which do not appear there:    Figure S2: Effects of mismatches between species on the consequences of environmental fluctuations on open, dynamically assembling food webs. Mismatch-level is defined as inverse of cross-correlation between species' sensitivity to fluctuations -so a correlation of 0 implies perfect mismatch between species across the whole community, while a correlation of 1 implies a perfect match such that all species are identically sensitive to fluctuations. A correlation of 1 (not shown) results in qualitatively similar results as those shown above. The bars represent 95% confidence intervals around the mean of 200 community assembly simulations. Each model community at the end of a simulation is at immigration-extinction equilibrium. Note that some food web features are significantly sensitive to mismatches, while others are not.    2) and (E) no stochastic fluctuations but strong and frequent extreme events (δ = .1, pδ = 1). Note that for all cases |λmax|≪ as expected from May's [6] analytical results. 1/|λmax| of a locally stable dynamical system captures the return time of the system to the equilibrium following perturbation, and is a measure of resilience. All the eigenvalues are negative because we only considered the locally-stable subset of the community -as such, all stochastically assembled communities have at least one species that is on a path to extinction. That is, no dynamically assembling system is locally stable, a result of the constant and inherent immigration-extinction dynamics of open food webs in dynamically changing environments such as running waters.