Sergio M Vallina
email: vallina@mit.edu

Corinne Le Quéré

Stability of complex food webs: Resilience, resistance and the average interaction strength

In the face of stochastic climatic perturbations, the overall stability of an ecosystem will be determined by the balance between its resilience and its resistance, but their relative importance is still unknown.

Using aquatic food web models we study ecosystem stability as a function of food web complexity.

We measured three dynamical stability properties: resilience, resistance, and variability. Specifically, we evaluate how a decrease in the strength of predator-prey interactions with food web complexity, reflecting a decrease in predation efficiency with the number of prey per predator, affects the overall stability of the ecosystem. We find that in mass conservative ecosystems, a lower interaction strength slows down the mass cycling rate in the system and this increases its resistance to perturbations of the growth rate of primary producers. Furthermore, we show that the overall stability of the food webs is mostly given by their resistance, and not by their resilience. Resilience and resistance display opposite trends, although they are shown not to be simply opposite concepts but rather independent properties. The ecological implication is that weaker predator-prey interactions in closed ecosystems can stabilize food web dynamics by increasing its resistance to climatic perturbations.

Introduction

The ability of an ecosystem to return to its reference state after a perturbation stress is given by its resilience [START_REF] May | Stability and complexity in model ecosystems[END_REF][START_REF] Harwell | Nutrient recycling and stability: A reevaluation[END_REF][START_REF] Pimm | Food Webs[END_REF]. A high resilient ecosystem is one that is able to recover fast after inbalances in the populations densities induced by climatic fluctuations. The ability of an ecosystem to resist displacement from its reference state during a perturbation stress is given by its resistance [START_REF] Webster | Nutrient cycling and the stability of ecosystems[END_REF][START_REF] Harwell | Nutrient recycling and stability: A reevaluation[END_REF][START_REF] Harrison | Stability under environmental stress: resistance, resilience, persistence, and variability[END_REF]. A high resistant ecosystem is one which is displaced slightly after imbalances in the populations rates induced by climatic fluctuations [START_REF] Harrison | Stability under environmental stress: resistance, resilience, persistence, and variability[END_REF][START_REF] Harrison | Resistance of nutrient cycling systems to perturbations of the flow rates[END_REF]. Both properties will influence how close an ecosystem remains to its reference equilibrium state and the ecosystem's variability in the face of stochastic climatic perturbations. In other words, resilience and resistance are complementary properties that will determine the overall stability of the ecosystem [START_REF] Ives | Stability and diversity of ecosystems[END_REF]. 27

Ecosystem resilience has been extensively studied and it is probably the most widely used metric for food web stability [START_REF] May | Will a large complex system be stable?[END_REF][START_REF] Pimm | Food Webs[END_REF][START_REF] Pimm | The complexity and stability of ecosystems[END_REF][START_REF] Loreau | A new look at the relationship between diversity and stability[END_REF]. The theory of Lyapunov stability predicts that ecosystems of higher complexity should be less stable in the face of perturbations [START_REF] May | Stability and complexity in model ecosystems[END_REF]. In this regard, systems are deemed L-stable when the density of all species returns to equilibrium following a perturbation [START_REF] Pimm | Food Webs[END_REF][START_REF] Pimm | The complexity and stability of ecosystems[END_REF]. Both the probability of displaying L-stable dynamics and the rate of recovering from a perturbation in population densities have been shown to decrease with ecosystem complexity in models [START_REF] May | Stability and complexity in model ecosystems[END_REF][START_REF] Fussmann | Food web complexity and chaotic population dynamics[END_REF]Chen and Cohen, 2001b;[START_REF] Ives | Stability and diversity of ecosystems[END_REF]. 35

Despite the importance of resistance for ecosystem stability [START_REF] Pimm | The complexity and stability of ecosystems[END_REF], theoretical works that address ecosystem resistance are much less abundant than in the case of resilience, probably due to the lack of a unified theory [START_REF] Harrison | Stability under environmental stress: resistance, resilience, persistence, and variability[END_REF]. Therefore, resistance has been a stability property more difficult to quantify [START_REF] Harwell | Nutrient recycling and stability: A reevaluation[END_REF][START_REF] Webster | Nutrient cycling and the stability of ecosystems[END_REF][START_REF] Harrison | Resistance of nutrient cycling systems to perturbations of the flow rates[END_REF][START_REF] Loreau | Phenotipic diversity and stability of ecosystem processes[END_REF]. Also, there is not yet a clear picture of how ecosystem resistance may be affected by food web diversity [START_REF] Loreau | Phenotipic diversity and stability of ecosystem processes[END_REF]. Nevertheless, there is some empirical evidence of increased resistance to environmental perturbations (e.g. drought) with biodiversity in natural grasslands (Tilman It is well know that resilience is a property that characterize the ecosystem as a whole [START_REF] Loreau | Phenotipic diversity and stability of ecosystem processes[END_REF]. However, contrary to previous suggestions that resistance is a property that characterizes each ecosystem component separately [START_REF] Harrison | Resistance of nutrient cycling systems to perturbations of the flow rates[END_REF][START_REF] Loreau | Phenotipic diversity and stability of ecosystem processes[END_REF], resistance is also a property that characterize the ecosystem as a whole. In a euclidean multi-dimensional space (e.g. multi-species ecosystem), both resilience and resistance can be measured in an analogous way from the time needed to go from one point to another point in that space during the event of single climatic perturbation. 50

Due to stochastic environmental fluctuations, populations in real ecosystems are almost never able to remain close to their equilibrium densities [START_REF] Ives | Measuring resilience in stochastic systems[END_REF]. However, disentangling resilience and resistance individually for stochastic ecosystems can be challenging. Resilience is a measure of the ecosystem stability to fluctuations in the population densities due to environmental perturbations [START_REF] Harrison | Stability under environmental stress: resistance, resilience, persistence, and variability[END_REF].

It does not give, however, any information about how hard might be for the environment to make these populations fluctuate. Resistance, on the other hand, is a measure of the ecosystem stability to fluctuations in the populations rates (e.g. specific growth or mortality) [START_REF] Harrison | Stability under environmental stress: resistance, resilience, persistence, and variability[END_REF]. Thus, it tells how difficult is for the environmental perturbations to make the populations fluctuate.

The main objective of this work is to evaluate whether equilibrium stability properties (i.e. resilience and resistance) can be related to stochastic stability measures (i.e. population variability) and what mechanisms are likely to explain the observed relationship between ecosystem complexity and overall stability. 62 Ecosystem resilience to perturbations has been shown to depend negatively on the number of species in a food web and the strength of competitive interactions between the species [START_REF] May | Will a large complex system be stable?[END_REF][START_REF] May | Stability and complexity in model ecosystems[END_REF] and positively on the speed of mass cycling through the system [START_REF] Deangelis | Energy flow, nutrient cycling, and ecosystem resilience[END_REF][START_REF] Deangelis | Effects of nutrient recycling and food-chain length on resilience[END_REF].

Equilibrium theory suggest that an increase in species richness and connectance should decrease the resilience of ecosystems because the inter-specific competition between species are destabilizing forces that will tend to push some species to extinction [START_REF] May | Stability and complexity in model ecosystems[END_REF][START_REF] Mccann | Weak trophic interactions and the balance of nature[END_REF]. However, the presence of weak interactions in complex food webs has been suggested to help sustain high ecosystem diversity [START_REF] Kokkoris | Variability in interaction strength and implications for biodiversity[END_REF][START_REF] Neutel | Reconciling complexity with stability in naturally assembling food webs[END_REF]. 70

Natural communities tend to display skewed distributions of the interaction strengths towards weak links, i.e. many weak interactions and few strong ones [START_REF] Paine | Food-web analysis through field measurements of per capita interaction strength[END_REF][START_REF] Mccann | Weak trophic interactions and the balance of nature[END_REF][START_REF] Neutel | Stability in real food webs: weak links in long loops[END_REF][START_REF] Neutel | Reconciling complexity with stability in naturally assembling food webs[END_REF], although the underlying mechanisms are not completely well understood [START_REF] Berlow | Interaction strengths in food webs: issues and opportunities[END_REF].

Among them, the number of prey per predator seems to correlate to the strength of the interactions; polyphages (i.e. predators with many prey) have weaker predator-prey interactions than monophages [START_REF] Mccann | Weak trophic interactions and the balance of nature[END_REF][START_REF] Montoya | Ecological networks and their fragility[END_REF]. 76

Little is known about the dependence of resistance on ecosystem complexity, species interaction strength and the rate of mass cycling through the food web. Therefore, it is unknown what determines the resistance of ecosystems to climatic perturbations. It is also unknown what is the relative weight of resilience and resistance on the overall stability of stochastic ecosystems. 80

We address these questions by means of multi-species food web models that are mass-conservative (i.e. closed ecosystems). Specifically, we evaluate how a decrease in the strength of predator-prey interactions with food web complexity will affect ecosystem resilience, resistance and overall stability. We will use as climatic perturbations drastic changes in solar radiation levels (i.e. dark/light pulses). They will therefore affect the specific growth rate of primary producers. 85 Food web resilience and resistance will be evaluated by applying a single perturbation (i.e. one dark pulse) to the ecosystems. Resistance will be measured as the time needed to reach a given non-equilibrium perturbed state from the reference equilibrium state. Resilience will be measured as the return rate from the non-equilibrium perturbed state back to the reference equilibrium state. Food web overall stability will be evaluated by applying a stochastic perturbation (i.e. many random dark pulses) to the ecosystems.

We will define a metric for overall stability that merges both the temporal variability of the system and its average distance to the reference state, and we will compare it to more classical measures of ecosystem variability such as the coefficient of variation of population-level, community-level and ecosystem-level properties. 94

Model description

96

We constructed a series of food web models that differ in their complexity. Food web complexity here refers to the number of plankton species or plankton functional types (PFT) (Le [START_REF] Quéré | Ecosystem dynamics based on plankton functional types for global ocean biogeochemistry models[END_REF]. A minimum of 3 and a maximum of 6 trophic levels is considered. Complexity increases along two axis in the model: 3 levels of phytoplankton complexity (p-cmpx) and 4 levels of zooplankton complexity (zcmpx). Therefore there are 12 food web configurations between the least and most complex ecosystem.

The most complex food web is based on the schematic food web described by [START_REF] Duffy | The 524 functional role of biodiversity in ecosystems: incorporating trophic complexity[END_REF]. The least complex model is a short linear food web. Model equations are given in Box 1 and Box 2. Model equation terms are listed in Table 1. Parameter values are listed in Table 2. The topology diagrams for the 12 food web complexities is given in Figure 1. Note that all food web complexities can be generated as substructures of the most complex food web. For phytoplankton we chose to go from large to small when increasing p-cmpx, but we obtain analogous results if we go from small to large phytoplankton (see Supp.Mat). For zooplankton we go from lower trophic to higher trophic levels when increasing z-cmpx.

Primary production is limited by the availability of dissolved inorganic nitrogen (DIN) following Michaelis-Menten nutrient uptake kinetics. Predation is assumed to follow a Holling-Type III (sigmoid) functional response, which is common in filter feeders and marine copepods [START_REF] Jeschke | Consumer-food systems: why type i functional responses are exclusive to filter feeders[END_REF][START_REF] Fussmann | A direct, experimental test of resource vs. consumer dependence[END_REF]. Type III functional responses for predation are known to stabilize model dynamics and to promote biodiversity through a relaxation of feeding at low prey densities, allowing for a much greater probability of co-existence of all the species in a food web [START_REF] Haydon | Pivotal assumptions determining the relationship between stability and complexity: An analytical synthesis of the stability-complexity debate[END_REF]; [START_REF] Gentleman | Functional responses for zooplankton feeding on multiple resources: a review of assumptions and biological dynamics[END_REF]; and references therein). The use of functional relationships including a relaxation of feeding at low prey densities is a way of implicitly taking into account several ecological mechanisms like prey refuge, prey switching and/or predator interference that are difficult to model explicitly. Alternatively, one can make use of Holling-Type II functional response with an explicit refuge for the prey. In any case, our results are not dependent on this choice (see Supp.Mat). Fish (top carnivore) is used in the model simply as a closure mortality term for zooplankton. All losses (i.e. mortality, non-assimilated food) are assumed to be instantaneously recycled to the nutrients DIN pool, making the ecosystem mass-conservative. The total amount of mass in the system is constant and it is they same for all food web configurations (N tot = 2 122 mmolN m -3 ). 123

The strength of zooplankton predation upon their prey is assumed to decrease with the number of prey per predator [START_REF] Mccann | Weak trophic interactions and the balance of nature[END_REF][START_REF] Montoya | Ecological networks and their fragility[END_REF]. Therefore, the model assumes that pressure on individual prey decreases when predators scavenge multiple prey. This can be understood as a way of implicitly reflecting an increase in the predators difficulty to attack their prey due to an increase in habitat complexity with species diversity [START_REF] Real | The kinetics of functional response[END_REF], such as heterogeneities in prey distribution (i.e. patches) or prey defense-strategies, which are not explicitly accounted for in the model. The decrease in predator-prey interaction strength comes from two complementary mechanisms. First, we assume that an increase in the number of prey per predator decreases the predator's efficiency in capturing any of its prey, which is parametrized as an increase in the predator's half-saturation constant for ingestion (eq (10) gives the predation efficiency as a function of the number of prey). This might be seen as a transition from a highly efficient specialist to a low efficiency generalist, or it can reflect an increase in the average transit time spent from going from one prey type patch to another [START_REF] Oaten | Transit time and density-dependent predation on a patchily distributed prey[END_REF]. Second, since each prey has its own implicit grazing refuge in the Type III functional response (or explicit if using Type II with refuge), the presence of more prey types implies a higher overall refuge for the prey as a whole.

Both mechanisms will tend to decrease the average interaction strength between predators and prey as food web complexity increases. 139

The modelled PFTs were allowed to differ in their phenotypic traits for the exploitation of resources.

Both the phytoplankton maximum growth rates and nutrient half-saturation constants, as well as the zooplankton maximum ingestion rates and half-saturation constants for ingestion, are assigned randomly from a normal distribution with mean values given in Table 2 and a standard deviation of 0.2 times the mean value. In practice this means a range for parameter values of approximately ± 50% the mean values. (Using a uniform distribution with a range ± 50% the mean values gave similar results). For each of the 12 food web complexities we performed 400 runs with randomly-assigned phenotypic traits. This amounts to a total of 12 x 400 = 4800 individual runs. We also performed an extra run in which all phytoplankton and all zooplankton share the same traits respectively (i.e. parameter values given in Table 2) which will be the control run for a parameters' sensitivity analysis (see Supp.Mat). For all runs the ecosystems have a single L-stable reference equilibrium under constant light levels, regardless of the initial conditions. Therefore, all food web complexities are globally L-stable [START_REF] Pimm | Food Webs[END_REF]Chen and Cohen, 2001a). 152

Dynamical Stability Indices 153

The only external forcing in the model is solar radiation. Therefore, we will use as climatic perturbations dark-pulses in solar radiation levels. The ecosystems will be allowed to be under only two possible environmental conditions: perturbed (full dark) or unperturbed (full light). Three dynamical stability indices (DSI) will be computed numerically: resilience, resistance, and overall stability. Resilience and resistance will be estimated by applying a single perturbation (i.e. one dark pulse) to the ecosystems. Overall stability will be evaluated by applying a stochastic perturbation (i.e. many random dark pulses) to the ecosystems. 160

Resilience Index 161

The reference steady-state point (SSP), or unperturbed equilibrium, is the point X * with coordinates 162

(X * 1 , ..., X * n )
in the phase space in which the concentration of the PFTs does not change over time. Since all our model food webs are globally Lyapunov stable, any departure from the reference state caused by a climatic perturbation will decay over time once the perturbation has ceased, and the ecosystems will asymptotically return back to their reference equilibrium. The resilience of the systems will be estimated by measuring the time needed for recovery, and then converting the return times into return rates. This has the advantage of giving a comparable measure to the dominant eigenvalue of the Jacobian matrix [START_REF] May | Stability and complexity in model ecosystems[END_REF]. 169

Model equations

•

Phytoplankton [mmolN m -3 ] ∂Pi ∂t = FP i -GP i Z -MP i (1) • Zoooplankton [mmolN m -3 ] ∂Zj ∂t = FZ j -GZ j Z -GZ j C -MZ j (2) • Top Carnivore Fish [mmolN m -3 ] ∂C ∂t = FC • (1 - C Kc ) -MC (3) • Dissolved Inorganic Nutrients [mmolN m -3 ] ∂DIN ∂t = - m i FP i + m i MP i + n j EZ j + n j MZ j + EC + FC • C Kc + MC (4)

Model terms

• Phytoplankton production (FP j )

FP i = μP i • Q DIN P i • Q P AR P i • Pi (5) Q DIN P i = DIN kP j + DIN (6) Q P AR P i = ( * ) (7)
(*) = 1 or 0 (unperturbed / perturbed)

• Zooplankton grazing / predation (GX i Z j )

GX i Z j = gZ j • Zj • X 2 i (kZ j /εZ j ) 2 + m i X 2 i ( 8
)
GX i Z = n j GX i Z j (9) εZ j = 1 m α (10)
• Zooplankton production (FZ j ) and excretion (EZ j )

FZ j = βZ • m i GX i Z j (11) EZ j = (1 -βZ) • m i GX i Z j (12)
• Top Carnivore Fish predation (GZ j C )

GZ j C = gC • C • Z 2 j k 2 C + n j Z 2 j ( 13
)
• Top Carnivore Fish production (FC) and excretion (EC)

FC = βC • n j GZ j C (14) EC = (1 -βC ) • n j GZ j C (15)
• Natural mortality of phytplankton (MP j ), zooplankton (MZ j ) and top carnivore fish (MC )

MP i = mP • (1 -e -γ P •P i ) • Pi (16) MZ j = mZ • (1 -e -γ Z •Z j ) • Zj (17) MC = mC • (1 -e -γ C •C) • C (18)
until the system reaches a perturbed non-equilibrium state X(t) = X * + x(t), where x(t) is a given 173 perturbation distance. The distance between the perturbed point and the steady-state point is given by the classical Euclidean norm [START_REF] Neubert | Alternatives to resilience for measuring the responses of ecological 584 systems to perturbations[END_REF]: 175

x t ≡ x(t) = x 2 1 (t) + x 2 2 (t) + ... + x 2 n (t) = x 2 i (t) (19)
The maximum potential distance is given by the distance between an extreme unperturbed point in which the total mass N tot in the system is allocated in a single PFT, and an extreme perturbed point in which all PFT concentrations tend to zero so that the nutrients pool in the system tends to N tot . Therefore: 178

x max 0 ≡ x(0) max = N 2 tot + 0 + ... + N 2 tot = 2 • N 2 tot ≈ N tot (20)
The dark-pulse perturbation will be in place until a given 10% of the maximum potential distance between the perturbed point and the steady-state point was attained (i.e. x 0 = δ 0 • x max 0

; with δ 0 = 0.1). 180

Then full light conditions will be back again and we will numerically compute the time needed by the system to recover (δ < δ 0 ) to a 99.99% level. From the return times we can calculate the return rates 182 [d -1 ] assuming an exponential decay of the perturbation distance: 183

x t = x 0 • e η•t (21)
where x 0 is the perturbation distance, x t is the return distance, t is the return time, and η is the return 184 rate. Making these distances a fraction of the maximum potential distance of a perturbation:

δ • x max 0 = δ 0 • x max 0 • e η•t (22) δ = δ 0 • e η•t (23)
Taking the logarithm on both sides and isolating η, we obtain the resilience index:

η = ln( δ δ 0 ) t (24)

Resistance Index 187

The resistance of the ecosystems to the dark-pulse perturbation was measured as the perturbation time (the length of the dark-pulse), i.e. the time needed by the system under the influence of the dark-pulse to reach the selected perturbation distance. Longer perturbation times imply a slower displacement per unit of time and thus a higher resistance to the perturbation [START_REF] Harrison | Stability under environmental stress: resistance, resilience, persistence, and variability[END_REF].

The ultimate displacement is the distance between the "all alive" reference equilibrium state and the "all died" equilibrium state that would occur if the disturbance lasted indefinitely [START_REF] Harrison | Resistance of nutrient cycling systems to perturbations of the flow rates[END_REF]. Although the ultimate displacement can be different for each ecosystem depending on the coordinates of its reference state, it can be shown that the differences are negligible and all the modelled food webs have similar ultimate displacements that are, in fact, very close to the maximum potential distance (i.e. N tot ). Therefore, the measured perturbation time is a dynamical resistance index related to the speed at which the ecosystems depart from the reference state towards its ultimate displacement.

Stability Index 199

To obtain an estimate of the overall stability of the food webs, we performed a simulation in which the ecosystems were submitted to a stochastic perturbation of light levels (i.e. dark/light pulses) during one year without seasonality. Each day had a probability of 50% to be assigned a zero light level. For those random dark-days, phytoplankton do not grow (perturbed states). For all the other days, phytoplankton growth was not light-limited (unperturbed states). The severity of every perturbation is given by the amount of consecutive dark-days. Prior to the start of the stochastic perturbation each food web was allowed to reach its reference equilibrium state after a year under no light-limited growth.

Ecosystems under stochastic perturbations are always in a non-equilibrium state [START_REF] Ives | Measuring resilience in stochastic systems[END_REF]. Every time that a perturbation comes in, the system moves away from its reference equilibrium state; once the perturbation is gone, the system moves back towards its reference equilibrium state. Since there is never enough time to fully recover before another perturbation comes in, the system is unable to reach the steady-state. Under those conditions the equilibrial state can be given by a stationary distribution characterized by its mean and standard deviation [START_REF] Ives | Estimating community stability and ecological interactions from time-series data[END_REF]. Specifically we will use the trajectories in the phase space, evaluated with respect to the reference equilibrium (i.e. displacements), as the stationary distribution of the stochastic system. The mean of the stationary distribution measures the average distance to the reference equilibrium. The standard deviation measures the temporal variability of the distances. 216

We will show that the distances between the stochastic non-equilibrium states and the reference equilibrium state do not follow a normal distribution. They are best described by a Weibull distribution, which is characterized by two parameters: the shape and the scale parameter [START_REF] Rinne | The Weibull distribution: a handbook[END_REF]. We fit a Weibull probability density function (PDF) to the stochastic distances, and then calculated the mean (μ) and standard deviation (σ) of the Weibull distribution [START_REF] Rinne | The Weibull distribution: a handbook[END_REF]. (The only restriction imposed on was that the shape parameter had to be ≥ 1.25 to avoid exponential distributions.) Low mean and low standard deviation of the stationary distribution means high stability of the stochastic ecosystem. We now define the parameter ϕ as the square root of the product between the mean and the standard deviation of the stationary Weibull distribution: 225

ϕ = √ μ • σ (25)
We use the inverse of parameter ϕ of the Weibull PDF as the overall stability index:

ξ = 1 ϕ (26)
To measure the strength of species interactions we obtained the Interaction matrix from the Jacobian matrix computed at the reference equilibrium [START_REF] Berlow | Interaction strengths in food webs: issues and opportunities[END_REF]. The elements c ij in the Interaction 229 matrix give the linear effect of a small change in the concentration of species j on the specific (i.e. per-capita) rate of change of species i: 231

c ij = ∂( Ẋi /X i ) ∂X j X * (27)
where Ẋi denotes the rate of change of species i and X i its concentration. The ecosystem average 232 interaction strength (EAIS) was obtained from the absolute value of the off-diagonal elements in the Interaction matrix: 234

EAIS = n i =j | c ij | n(n -1) (28)
where n is the total number of species in the food web. Note that the interactions with nutrients were not included in the calculations. 236

We chose to use the elements of the Interaction matrix [m -3 mmolN -1 d -1 ] ( [START_REF] Kokkoris | Patterns of species interaction strength in assembled theoretical competition communities[END_REF] over 237 the elements in the Jacobian matrix [d -1 ] ( [START_REF] Neutel | Stability in real food webs: weak links in long loops[END_REF] because they give per-capita interaction 238 strengths, i.e. equivalent to the interaction coefficients in the Lotka-Volterra model [START_REF] Wootton | Estimates and test of per capita interaction strength: diet, abundance, and impact of intertidally foraging birs[END_REF], and therefore measure a level of interaction between species that is independent of the equilibrium concentrations, which may change from one food web to another. This provides a more objective measure of the imposed decrease in species interaction with complexity, since variations in interaction strength resulting from differences in species' densities are circumvented [START_REF] Laska | Theoretical concepts and empirical approaches to measuring interaction strength[END_REF]. This metric has been extensively used to study the stability of both real food webs [START_REF] Emmerson | Predator-prey body size, interaction strength and the stability of a 529 real food web[END_REF] and ecosystem models [START_REF] Kokkoris | Patterns of species interaction strength in assembled theoretical competition communities[END_REF][START_REF] Jansen | Complexity and stability revisited[END_REF]. In any case, using instead the elements in the Jacobian to measure the strength of species interactions did not change significantly the results since the correlation between average Jacobian Matrix strength and average Interaction Matrix strength is very high (see Supp.Mat). 248

In this section we will show the relationship between the ecosystems' state, productivity and stability properties as a function of the food web complexity, which is defined along two axis as the number of phytoplankton and zooplankton species present in each food web (p-cmpx and z-cmpx, respectively; see methods). For each of the 12 food web complexities, we take the average of 400 runs with randomlyassigned phenotypic traits. 254

Figure 2 shows the steady-state concentrations as a function of phytoplankton and zooplankton complexity. Total phytoplankton biomass increases (about a 70%) along its own phytoplankton complexity (p-cmpx) axis from ≈0.7 to ≈1.2 mmolN m -3 , but remains almost unchanged along the zooplankton 257 complexity (z-cmpx) axis. Total zooplankton biomass decreases (about a 30%) along the phytoplankton complexity axis from ≈1.0 to ≈0.70 mmolN m -3 , but slightly increases (about 15%) along its own 259 zooplankton complexity axis. Top carnivore fish biomass displays a moderate decrease both along the phytoplankton and zooplankton complexity axes (about a 30-40%), going from ≈0.25 to ≈0.10 mmolN m -3 between the simplest and most complex food web. Finally, dissolved inorganic nutrients (DIN) also decreases both along the phytoplankton (about a 80%) and zooplankton (about a 40%) complexity axes, going from ≈0.2 to ≈0.01 mmolN m -3 between the simplest and most complex food web. This implies 264 a higher degree of retention of nutrients by the organisms with food web complexity.

The increase of total phytoplankton along the p-cmpx axis can be mainly attributed to the imposed lower predation pressure upon phytoplankton when more species are added. For the same reason, this increase is mirrored by a small decrease in total zooplankton along the p-cmpx axis. Along the z-cmpx, however, since zooplankton can be at the same time both predator and prey, the lower food intake is balanced by a lower predation mortality and therefore the total zooplankton increase is much smaller.

Figure 3 gives the rate of mass input fluxes between compartments at the reference steady-state as a function of phytoplankton and zooplankton complexity. Between the simplest and most complex food web, primary production decreases by more than 50% (from ≈0.4 to ≈0.2); zooplankton production decreases by more than 60% (from ≈0.12 to ≈0.04); fish production decreases by more than 80% (from ≈0.03 to ≈0.005); and the recycling production of DIN decreases about 70% (from ≈0.9 to ≈0.3).

Thus, although the degree of retention of nutrients increases with food web complexity (Figure 2d), the rate of instantaneous production in the system decreases (Figure 3) due to the imposed weaker predatorprey interactions when more species are added. Weaker predator-prey interactions means slower recycling of DIN from zooplankton excretion, slowing down the pace of phytoplankton primary production.

Thus, assuming that predator-prey interactions are weaker with more prey choices results in a decrease with food web complexity of the speed of mass cycling in the ecosystem.

Figure 4a gives the ecosystems' average interaction strength as a function of phytoplankton and zooplankton complexity. The EAIS decreases strongly along both axes, which is to be expected since the model assumes a general decrease in the predator-prey interactions with food web complexity. The average interaction goes from ≈0.25 to ≈0.1 between the simplest and most complex food web, about a 60% decrease. Figures 4b,c,d give the food webs' DSI: resilience, resistance, and overall stability; respectively. We found that resilience decreases (about 75%) with food web complexity (Figure 4b), with the rate of return (i.e. 99.99% recovery) going from ≈0.2 to ≈0.08 [d -1 ] between the simplest and most 288 complex food web. On the other hand, resistance increases (about a 200%) with food web complexity (Figure 4c), with the perturbation time going from ≈0.4 to ≈1.2 [d] between the simplest and most complex food web. This means that in order to displace the ecosystems the same perturbation distance (see section 2.2.1), the most complex food web can resist the climatic perturbation 3 times longer than the simplest one. The overall stability also increases (about a 200%) with food web complexity, going from ≈3.0 to ≈9.0 [m 3 mmolN -1 ] between the simplest and most complex food web. The similarity between how the resistance index and the stability index change with food web complexity is apparent, which already indicates that the overall stability of the food webs seems more related to their resistance than to their resilience. 297

Figure 5 shows the probability density functions of the distances between the stochastic ecosystem and the reference equilibrium for each level of phytoplankton (top-down) and zooplankton complexity (left-right). We see that as food web complexity increases, both the mean and the standard deviation of the Weibull PDF decrease. This means that the stochastic ecosystem remains closer to the reference equilibrium and display lower temporal variability when food web complexity is higher. In other words, the ecosystems' stability increases with food web complexity. 303

Figure 6 shows the dynamics of the stochastic ecosystems just for the simplest and most complex food web. In this case the results are from the single control run (see methods), not the average values from the random runs, since the figure is illustrative only. (Note that all the stability properties of the control run are almost identical to the averaged stability properties of the random runs, see Supp.Mat).

Figure 6a shows the trajectories in the euclidean phase space of these two stochastic ecosystems (P.Z and PPP.ZZZZ). Figure 6b show the corresponding time-series. The PPP.ZZZZ stochastic ecosystem is 309 clearly more stable than the P.Z stochastic ecosystem; it displays both lower temporal variability and its trajectories remain closer to the reference equilibrium. 311

Figure 7 shows the relationship between the dominant eigenvalue of the Jacobian matrix versus the average interaction strength and the food webs' DSI. Of these four metrics, only the return rate (i.e. the resilience index) correlates well with the dominant eigenvalue, which is to be expected since in fact they are both measuring the same property (the resilience of the ecosystems). In all other cases, the relationship between the resilience measured by the dominant eigenvalue versus EAIS, resistance, and overall stability, shows negative but not significant trends. One important conclusion is that resilience and resistance are not simply inverse concepts but rather appear to be independent properties, as previously suggested [START_REF] Harrison | Stability under environmental stress: resistance, resilience, persistence, and variability[END_REF]. 319

Likewise, figure 8 shows the relationship between the total primary production in the ecosystem versus the average interaction strength, food webs' DSI resistance and overall stability, and the average rate of mass cycling in the system (excluding primary production). All four metrics correlate strongly with primary production. The correlation is positive, which means that a decrease in primary production is related to a decrease in EAIS, a decrease in the speed of mass cycling, and an increase in resistance and overall stability (note that the plots are against the inverse of these two stability indices).

Finally, figure 9 shows the relationship between the parameter ϕ of the Weibull distribution (i.e. the inverse of the stability index) and the coefficient of variation (CV) for the individual populations of phytoplankton and zooplankton species (figure 9a); the aggregated populations of total phytoplankton and total zooplankton (figure 9b); the average rate of mass cycling in the system (excluding primary production) (figure 9c); plus versus the parameter ϕ calculated using the mean and standard deviation of the raw data of the stochastic displacements, i.e. without fitting a Weibull distribution (figure 9d). All four metrics correlate well with parameter ϕ of the Weibull distribution. Again, the correlation is positive, which means that an increase in the stability index (note the Weibull ϕ is the inverse of the stability index) is related to an increase in the stability of properties at the population-level, the community-level and the ecosystem-level (i.e. time-series of individual populations, aggregated populations, and mass cycling rate; respectively). Also it shows that our definition of stability index is not sensitive to the choice of fitting a Weibull distribution to the raw data. 337

Discussion

338

The importance of the strength of species interactions for the stability of complex ecosystem has been studied extensively [START_REF] Yodzis | The stability of real ecosystems[END_REF][START_REF] Mccann | Weak trophic interactions and the balance of nature[END_REF][START_REF] Ives | Complex dynamics in stochastic tritrophic models[END_REF][START_REF] Kokkoris | Variability in interaction strength and implications for biodiversity[END_REF][START_REF] Jansen | Complexity and stability revisited[END_REF][START_REF] Neutel | Stability in real food webs: weak links in long loops[END_REF][START_REF] Neutel | Reconciling complexity with stability in naturally assembling food webs[END_REF][START_REF] Berlow | Interaction strengths in food webs: issues and opportunities[END_REF][START_REF] Emmerson | Predator-prey body size, interaction strength and the stability of a 529 real food web[END_REF]. Many works have evaluated the effect of interaction strength on the resilience of ecosystems [START_REF] Ives | Complex dynamics in stochastic tritrophic models[END_REF][START_REF] Neutel | Stability in real food webs: weak links in long loops[END_REF][START_REF] Emmerson | Predator-prey body size, interaction strength and the stability of a 529 real food web[END_REF][START_REF] Neutel | Reconciling complexity with stability in naturally assembling food webs[END_REF], the inhibition of chaotic oscillations [START_REF] Mccann | Weak trophic interactions and the balance of nature[END_REF] or the vulnerability of communities to species invasions [START_REF] Kokkoris | Patterns of species interaction strength in assembled theoretical competition communities[END_REF]. Weaker interactions were shown to stabilize ecosystem dynamics and this has been suggested to help sustain high levels of species diversity [START_REF] Kokkoris | Variability in interaction strength and implications for biodiversity[END_REF][START_REF] Neutel | Reconciling complexity with stability in naturally assembling food webs[END_REF]. We do not find, however, a clear robust relationship between EAIS and the resilience index, and if anything the trend is negative. That is, decreasing the species interaction strength with complexity does not preclude the ecosystem to become less resilient. This supports previous findings that ecosystem resilience decreases with food web complexity [START_REF] May | Will a large complex system be stable?[END_REF][START_REF] Ives | Stability and diversity of ecosystems[END_REF].

On the other hand, we are not aware of any work addressing the effect of interaction strength on the resistance of complex ecosystems to climatic perturbations. Since the overall stability of an ecosystem is determined by the balance between its resilience and its resistance [START_REF] Harrison | Stability under environmental stress: resistance, resilience, persistence, and variability[END_REF][START_REF] Ives | Stability and diversity of ecosystems[END_REF], it seems clear that there is a gap of knowledge about an important component of ecosystem stability. We do find that a decrease of species interaction strength with complexity produces an increase of ecosystem resistance, and that this promotes a higher overall stability of the system. The relationship between EAIS and the resistance index, and between resistance and the overall stability, are both robust and positive. That is, decreasing the species interaction strength with complexity increases the overall stability of the ecosystem. 359

The mechanism that explains the increase in overall stability with ecosystem complexity is the decrease in the speed of mass cycling due to the weakening of the species interaction strength with food web complexity, which causes an increase in the ecosystem resistance. In contrast, resilience has been previously suggested to increase with the rate of mass cycling through the food web [START_REF] Deangelis | Energy flow, nutrient cycling, and ecosystem resilience[END_REF][START_REF] Deangelis | Effects of nutrient recycling and food-chain length on resilience[END_REF]. While our results still support this view, they show that resilience and resistance are not just exact opposite properties. We find no clear positive relationship between the rate of mass cycling through the food webs and their resilience, at least not as clear as the negative relationship we find between mass cycling and resistance. It turns out that other effects, like competitive interactions among species, may be more important to determine the level of resilience than the mass turnover rate.

The reason why mass cycling rate and resistance are mechanistically connected is simple. At low cycling rates the ecosystem dynamical state is one of slow motion. Therefore, it will react slowly to imbalances in population parameters induced by climatic perturbations. In fact, affecting some population parameters amounts to a change in the fluxes between some compartments in the system [START_REF] Harrison | Resistance of nutrient cycling systems to perturbations of the flow rates[END_REF]. In our dark perturbations, the primary production is halted altogether (the mass flux between nutrients and phytoplankton is zero). The subsequent decline in population densities is going to be mainly governed by the speed of mortality due to predation, and because at the reference state all rates are slow, mortality by predation is slow as well. Thus, it will take longer for the ecosystem to reach the "all dead" ultimate equilibrium, which means it is dynamically more resistant.

The resistance of a nutrient cycling system to a disturbance that alter population rates, like a climatic perturbation, depends both in the ultimate displacement that would occur if the disturbance lasted indefinitely and the rate at which the system approaches the ultimate displacement [START_REF] Harrison | Resistance of nutrient cycling systems to perturbations of the flow rates[END_REF]. This may be an issue when comparing the resistance to perturbations for different ecosystems. For example, one ecosystem may have a slow approaching rate (high dynamical resistance) but with a very large ultimate displacement. If the perturbation last long enough, this ecosystem may thus be classified as less resistant than another ecosystem displaying fast approaching rates (low dynamical resistance) but with a small ultimate displacement, even if dynamically the former is more resistant than later. Variations in the ultimate displacement can partly explain the very complex behaviour of resistance with diversity previously reported for a mass-conservative ecosystem [START_REF] Loreau | Phenotipic diversity and stability of ecosystem processes[END_REF].

Our results, however, do not have such confounding effects. The dark-pulse perturbations ensures that the ultimate displacement is fairly similar for all the ecosystems. For long enough periods of darkness, all the mass in the system will be recovered in the dissolved inorganic nutrients pool. The euclidean distance to this "all dead" system equilibrium state from the "all alive" unperturbed reference equilibrium state (i.e. the ultimate displacement distance), although depending slightly on the particular position of the reference equilibrium for each food web configuration, can be shown to be approximately the total mass in the system for all of the food webs. Therefore, the resistance index we measure is a dynamical resistance index; it is directly related to the speed at which the ecosystems approach the same ultimate displacement. 396

The ecosystem resistance at the unperturbed reference equilibrium has been suggested to be related to the resilience at the perturbed equilibrium [START_REF] Harrison | Resistance of nutrient cycling systems to perturbations of the flow rates[END_REF]. Furthermore, for very short ultimate displacements, the unperturbed equilibrium and the perturbed equilibrium are so close that the ecosystem dynamical resistance may be well characterized by its resilience at the unperturbed equilibrium alone. In such case the dominant eigenvalue of Jacobian matrix at the unperturbed reference equilibrium gives the dominant "response time" of the ecosystem to perturbations, and resilience and resistance can be seen as opposite symmetric properties. Quick response times would imply high resilience and low resistance, i.e. an ecosystem that will closely track climatic fluctuations [START_REF] Harrison | Resistance of nutrient cycling systems to perturbations of the flow rates[END_REF].

However, for longer ultimate displacements like in our simulations, the information given by the Jacobian matrix at the reference state may not serve to evaluate the ecosystem resistance at this point. In these cases the dominant eigenvalue at the reference state cannot tell much about the speed of displacement towards the ultimate displacement, because the unperturbed and perturbed equilibriums are so far apart that the might well have very different resilience. Furthermore, even information about the local resilience at the perturbed equilibrium might not be related to the resistance at the unperturbed reference equilibrium because the Jacobian's linear approximation around the perturbed equilibrium is likely to break down for long distances between these two points. Nevertheless, we tried to calculate the Jacobian matrix at the "all dead" no-light perturbed equilibrium but we obtained an array of zeros for all food web complexities. Therefore we were unable to relate the ecosystem resistance at the unperturbed reference equilibrium to the ecosystem resilience at the perturbed equilibrium.

Yet, we find that information about the mass-cycling dynamical state of the reference equilibrium explains the observed variability in the resistance index. Our results show that the speed of displacement towards the ultimate displacement is regulated by the rate of mass cycling through the food web at the unperturbed equilibrium. Slow cycling of mass implies a high resistance to climatic disturbances. And because the ecosystem is closed and in mass balance, total primary production alone gives a good measure of the average flux between compartments. This is the same to say that primary production gives a measure of the mass turnover rate of the system at equilibrium. Thus, primary production might well be used as a proxy of ecosystem resistance to climatic perturbations affecting population rates.

Recycling of nutrients has been argumented to be a positive feedback [START_REF] Webster | Nutrient cycling and the stability of ecosystems[END_REF]. Less intense nutrient recycling from predation losses (e.g. unassimilated mass) means that less nutrient is immediately available for uptake by phytoplankton and thus instantaneous primary production is decreased [START_REF] Loreau | Coexistence of multiple food chains in a heterogeneous environment: interactions among community structure, ecosystem functioning, and nutrient dynamics[END_REF], which in turn slows down the whole cycling of mass through the system and therefore the mass turnover rate. Our results show that reducing predation strength with food web complexity leads to a slow nutrient recycling, slow primary production, slow mass turnover rate, and a high resistance. All this while maximizing abiotic resource (nutrients) retention by the ecosystem as a whole, i.e. a high ratio of total biotic to abiotic mass. 431 ronmental grain theory [START_REF] Nunney | The stability of complex model ecosystems[END_REF]. In coarse-grained environments (i.e. species tending to form patches), predators must partition their search time between different prey species, lowering the functional response of a predator to each of its prey [START_REF] Nunney | The stability of complex model ecosystems[END_REF]. However, whether or not higher number of species in natural aquatic communities may lead to a decrease in the average predation strength still remains speculative at this point, although this seems to be the case in terrestrial ecosystems [START_REF] Montoya | Ecological networks and their fragility[END_REF]. Nevertheless, the equations describing our aquatic food webs are analogues to those used for terrestrial ecosystems [START_REF] Fussmann | Food web complexity and chaotic population dynamics[END_REF]. Therefore our results could be generalized to any mass-conservative ecosystem, i.e. food webs for which primary production is mostly fuelled by the recycling of nutrients within the system. 441

In any case, increasing prey chances of survival to predation with more species diversity seems a simple mechanism by which an increase in food web complexity can lead to both an efficient nutrient retention by the organisms in the food web and to a high resistance of the ecosystem itself to climatic fluctuations.

It should be stressed, however, that the dynamical mechanistic link is only between mass cycling rate and resistance. The relationship between resistance and the average interaction strength just reflects that the mass cycling rate can be affected by the level of predator-prey interactions; there is no dynamical mechanistic link between resistance and the average interaction strength as such.

Furthermore, these results are not restricted to changes in predation efficiency with food web complexity.

The relationship between EAIS, primary production, mass turnover, and ecosystem resistance still holds when we remove any influence of food web complexity and the only source of variability comes from the randomly assigned parameter values (see Supp.Mat). This implies that any kind of source of variability affecting the mass turnover rate in the ecosystem will affect the resistance of the food web. Thus, it could be argued that ecosystems with slow primary production, regardless of its cause, should be more resistant to climatic changes. Ongoing research is addressing this hypothesis by means of a global marine ecosystem model where the number of plankton species is an emergent property [START_REF] Follows | Emergent biogeography of microbial communities in a model ocean[END_REF][START_REF] Dutkiewicz | Modeling the coupling of ocean ecology and biogeochemistry[END_REF][START_REF] Barton | Patterns of diversity in marine phytoplankton[END_REF]. 457

Theoretical studies predict that community-level variability should increase with species richness, while population-level variability should decrease [START_REF] Ives | Stability and species richness in complex communities[END_REF][START_REF] Cottingham | Biodiversity may regulate the temporal variability of ecological systems[END_REF]. The stability index we defined is based on the stationary distribution of the stochastic ecosystem [START_REF] Ives | Estimating community stability and ecological interactions from time-series data[END_REF].

Therefore, it gives a measure of population-level stability that characterize the system as a whole. Nevertheless, when we compared it to other measures of ecosystem variability such as the coefficient of variation at the population-level, community-level and the ecosystem-level, we obtain that all these measures of stability give similar results. This implies that in our simulations there is no significant difference between at which level (population, community, ecosystem) we evaluate the food web stability. That means that a higher stability with food web complexity at the population-level causes a higher stability at the community-and ecosystem-level as well. 467

What is more important for the overall stability of an ecosystem to climatic perturbations, its resilience or its resistance? Our results suggest that resistance dominates. The variability of the food webs analyzed is so buffered when they are resistant to climatic perturbations, that their resilience becomes unimportant.

In fact, the highest overall stability was found for the lowest resilient ecosystem. Furthermore, the positive relationship between the resistance index and the overall stability index is almost as clear as the one between mass turnover and resistance. There is no clear relationship between the resilience index and the stability index but the trend is negative. We can then conclude that, regardless of their resilience, persistent ecosystems with slow mass cycling through the food web are very stable thanks to their increased resistance. 476

While our conclusions apply to perturbations of population rates (i.e. model parameters), they do not necessarily apply to perturbations of population densities (i.e. model state-variables). The stability of an ecosystem to inputs/outputs of its nutrients or organisms will not be related to the dynamical resistance described here. For these type of perturbations resilience should dominate the overall stability of the system. Also, the dynamical stability properties analyzed here are not related to structural stability properties like robustness to species extinctions or vulnerability to species invasions. Being dynamically stable to climatic disturbances does not necessarily mean being also structurally stable to species losses or invasions. 484

We have analysed three dynamical stability properties of mass-conservative food web models: resilience, resistance, and variability. The goal was to evaluate what determines the overall stability of complex ecosystems to stochastic climatic perturbations affecting population rates. Although classical theory of ecosystem resilience predicts that the stability of food webs should decrease with its complexity, our results demonstrate that this is not necessarily the case. We do find that resilience decreases with food web complexity. However, resilience only measures measures the speed of recovery from fluctuations in population densities, it does not provide information about how difficult it may be to make the populations fluctuate in the first place. This ecosystem sensitivity to climatic perturbations is given by the stability property resistance. We find that if an increase of food web complexity is combined with a decrease in the strength of predator-prey interactions, the overall stability of the ecosystem does in fact increase with complexity. A lower average interaction strength causes a lower turnover rate of mass in the ecosystem, and this increases its resistance to climatic perturbations. We evaluated the relative weight of resilience and resistance on the overall stability of the food webs and show that resistance dominates. These results suggest that more research should address the resistance of natural ecosystems to climatic perturbations, since this stability property appears to be key for the overall stability of complex food webs and seems to have been overlooked in previous theoretical studies. 501 Grazing on phytoplankton i from zooplankton 9 mmolN m -3 d -1 GZ j Z Predation on zooplankton j from zooplankton 9 mmolN m -3 d -1 GZ j C Predation on zooplankton j from top carnivore fish 13 mmolN m

-3 d -1 MP i Natural mortality of phytoplankton i 16 mmolN m -3 d -1 MZ j Natural mortality of zooplankton j 17 mmolN m -3 d -1 MC Natural mortality of top carnivore 18 mmolN m -3 d -1

Figures

  Figures

Table 1 :

 1 List of model equations' terms

	Parameter	Symbol Value Units
	Phy. max. specific growth rate Phy. half-sat. for DIN uptake Phy. natural mortality specific rate	μP i kP i mP	1.0 0.1 0.05	d -1 mmolN m -3 d -1
	Phy. natural mortality density-dependence γP	2.0	m 3 mmolN -1
	Zoo. max. specific ingestion rate Zoo. half-sat. const. for ingestion	gZ j kZ j	1.0 0.7	d -1 mmolN m -3
	Zoo. predation efficiency const.	α	1/3	adim.
	Zoo. assim. efficiency	βZ	1/3	adim.
	Zoo. natural mortality specific rate	mZ	0.05	d -1
	Zoo. natural mortality density-dependence γZ	2.0	m 3 mmolN -1
	Fish max. specific ingestion rate	gC	1.0	d -1
	Fish half-sat. const. for ingestion	kC	1.0	mmolN m -3
	Fish assim. efficiency	βC	1/3	adim.
	Fish natural mortality rate	mC	0.2	d -1
	Fish natural mortality density-dependence γC	2.0	m 3 mmolN -1
	Fish carrying capacity	KC	0.5	mmolN m -3

Table 2 :

 2 List of model parameters.

Mean value μ from a Gaussian distribution with standard deviation σ = 0.2 μ

Acknowledgments 502 Our special thanks to Anje M. Neutel, Robert A. Armstrong, and many other participants of the "7 th 503 Green Ocean Workshop, 2008" held at Villefranche-sur-Mer (France) for their useful comments and fruitful discussions. We also would like to thank two anonymous reviewers for their very relevant critics and suggestions which have highly improved the quality of this work. This research was supported by the Natural Environment Research Council (NERC grant number R8/H10/19) through a postdoctoral fellowship and by the European Comission through a Marie Curie OIF fellowship (both to S. M. Vallina).

Figure legends

• Fig. 1 Topology diagrams for the 12 food web complexities. Note that only predator-prey interactions are shown (i.e. nutrients and recycling arrows are not included for simplicity). Color nodes represent food web plankton functional types (PFT): green for phytoplankton (small, medium, large; from left to right); blue for micro-zooplankton; black for meso-zooplankton; yellow for macro-zooplankton; red for mega-zooplankton; and gray for top carnivore fish. The size of the nodes indicates the size of the PFT (not true scale). The color of the lines indicate predation on a given prey, except for gray dashed lines which represent predation by top predator fish on all zooplankton (closure term). 620 • Fig. 9 Relationship between the inverse of the stability index and: a) coefficient of variation for