Strong nutrient-plant interactions enhance the stability of ecosystems

Modular food web theory shows how weak energetic fluxes resulting from consumptive interactions plays a major role in stabilizing food webs in space and time. Despite the reliance on energetic fluxes, food web theory surprisingly remains poorly understood within an ecosystem context that naturally focuses on material fluxes. At the same time, while ecosystem theory has employed modular nutrient-limited ecosystem models to understand how limiting nutrients alter the structure and dynamics of food webs, ecosystem theory has overlooked the role of key ecosystem interactions and their strengths (e.g., plant-nutrient; R-N) in mediating the stability of nutrient-limited ecosystems. Here, towards integrating food web theory and ecosystem theory, we first briefly review consumer-resource interactions (C-R) highlighting the relationship between the structure of C-R interactions and the stability of food web modules. We then translate this framework to nutrient-based systems, showing that the nutrient-plant interaction behaves as a coherent extension of current modular food web theory; however, in contrast to the rule that weak C-R interactions tend to be stabilizing we show that strong nutrient-plant interactions are potent stabilizers in nutrient-limited ecosystem models.


A. Basic C-R Module
For the analysis of the basic C-R module, we use a classic Rosenzweig-MacArthur C-R system, with type II functional responses and logistic resource growth. The specific equations we use are as follows: where is resource biomass and is consumer biomass. For the logistic growth of the resource, we have intrinsic growth rate and carrying capacity . To model the interactions between the species we have the attack rate of the consumer on the resource , the assimilation rate , and half-saturation density 0 . The mortality rate of the consumer is .
Time series (Fig. 2) were conducted over 300 integration steps, with the resource initialized at * while the consumer was initialized at

B. P-C-R Food Chain Model
We extended the basic C-R module to a three-compartment food chain model, adding a top predator with a type II functional response. The specific equations were as follows: where is resource biomass, is consumer biomass, and is predator biomass. For the logistic growth of the resource, we have intrinsic growth rate and carrying capacity . To model the interaction between species we have two attack rates: consumer on the resource and predator on the consumer , with an assimilation rate , and two half-saturation densities: the consumer on the resource 0 and the predator on the consumer 0 . Finally, we have two mortality rates: the consumer , the predator .

C. P-C1-C2-R Specialist Model
We extended the P-C-R food chain model to a four-compartment model, adding a second consumer with a type II functional response. The specific equations were as follows: where is resource biomass, 1 is consumer 1 biomass, 2 is consumer 2 biomass and is predator biomass. For the logistic growth of the resource, we have intrinsic growth rate and carrying capacity . To model the interaction between species we have three attack rates: consumer 1 on the resource 1 , consumer 2 on the resource 2 and predator on consumer 1 1 , with an assimilation rate , and two half-saturation densities: consumer on the resource 0 and the predator on consumer 0 . Finally, we have three mortality rates: consumer and the predator .

D. P-C1-C2-R Generalist Model
We extended the P-C1-C2-R model to a diamond model, having a predator with a multi-species functional response. The specific equations were as follows: 1 + 1 ℎ 1 1 + 2 ℎ 2 2 = ( 1 1 + 2 2 ) 1 + 1 ℎ 1 1 + 2 ℎ 2 2 − where is resource biomass, 1 is consumer 1 biomass, 2 is consumer 2 biomass and is predator biomass. For the logistic growth of the resource, we have intrinsic growth rate and carrying capacity . To model the interaction between species we have four attack rates: consumer on the resource and the predator on consumer , with an assimilation rate , a handling time of consumer by the predator ℎ and a half-saturation density of consumer on the resource 0 . Finally, we have three mortality rates: consumer and the predator .

A. R-N Module
For the analysis of the R-N module, we use a limiting-nutrient pool and resource that are open to the environment, with Monod-like nutrient uptake by the resource. The specific equations we use are as follows: where is limiting-nutrient pool nutrients and is nutrients assimilated by the resource. The limiting-nutrient pool has external inputs . To model the uptake of nutrients we have the maximum rate of nutrient uptake by the resource and a half saturation density of . Nutrients are lost from compartment according to .

B. Local Stability Analyses
The local stability analyses for the model was based around the interior equilibrium: * = ( − * )( + * ) * * = − We impose biological feasibility on our system with the following restrictions on parameters: (1) parameter values must be positive and (2) parameterizations must confer positive abundances of * and * , such that the following inequalities are satisfied: With this restriction, we now consider stability criteria outlined in Strogatz (2018) to determine local stability (i.e., the sign of the real parts of the eigenvalues) about the interior equilibrium.
We first linearize the system about the interior equilibrium, giving the Jacobian matrix: The trace of matrix is given by: which satisfies the inequality < 0 for all biologically feasible parameterizations.
The determinant of matrix A is given by: which can be shown to satisfy the inequality det( ) > 0 as follows: Given the inequality > , ( − ) must be negative, and the inequality det( ) > 0 can be satisfied if which rearranges to: which has been satisfied by our restriction for biological feasibility (2.03). Therefore, given that the trace and determinant of the Jacobian for the interior equilibrium satisfy the criteria < 0 and det( ) > 0 for biologically feasible parametrizations, we can conclude that the dominant eigenvalue has real parts that are only negative and that the R-N system is locally stable.

C. Interaction Strength and Local Stability
We determine how local stability (i.e., the magnitude of the dominant eigenvalue) changes as is increased (Fig. 4). Again, we restrict our parameters for biologically feasibility as outlined in the previous section, satisfying inequalities (2.02) and (2.03). As is increased from the point at which the resource can persist, three qualitatively distinct patterns of stabilization occur depending on which of the following inequalities/equality are satisfied: > , = , < ( Fig. 4). We are able to prove the case of = as general with the following proof.
We first expand the determinant det( − ) = 0 which yields the characteristic equation: with the solutions to 1,2 given by: Substituting and for in equation (2.11) and rearranging gives: simplifying to (2.10) (2.11) (2.12) where 1 evaluates the positive root and 2 evaluates the negative root. Evaluating 2 as approaches ∞ yields: where 1 evaluates the positive root and 2 evaluates the negative root. Therefore, must tend from 0 to − as is increased.
Numerical analysis reveals that the horizontal asymptote can be approached from below or above, with approaching from above for > and from below for < (Fig. 4).  1,1). In all cases, has only real parts.

D. Nutrient Loading and Local Stability
We extended the analysis of interaction strength and stability by implicitly strengthening the R-N interaction through nutrient loading (i.e., forcing ) and asking how stability changes.
Starting with the special case of = , equation (2.12) gives the following: where 1 evaluates the positive root and 2 evaluates the negative root of equation ( where 1 evaluates the positive root and 2 evaluates the negative root. Therefore, must tend from 0 to − as is increased. Numerical analysis reveals that the horizontal asymptote can be approached from below or above, with approaching from above for > and from below for < ( Supplementary   Fig. 1). Parameters for Supplementary Fig. 1 1,1). In all cases, has only real parts.

A. C-R-N Model
We extend the basic R-N module to a three-compartment nutrient-limited model, adding a consumer with a type II functional response. The specific equations are as follows: where is limiting-nutrient pool nutrients, is nutrients assimilated by the resource, and is nutrients assimilated by the consumer. The limiting-nutrient pool has external inputs . To model the uptake of nutrients we have two maximum uptake rates: uptake from the limiting-nutrient pool by the resource and uptake from resource by the consumer , with an assimilation efficiency , and two half saturation densities: the resource on the limiting nutrient pool and the consumer on the resource 0 . Nutrients are lost from compartment according to .

B. P-C-R-N Model
We extended the C-R-N model to a four-compartment nutrient-limited model, adding a top predator with a type II functional response. The specific equations are as follows: where is limiting-nutrient pool nutrients, is nutrients assimilated by the resource, is nutrients assimilated by the consumer, and is nutrients assimilated by the predator. The limiting-nutrient pool has external inputs . To model the uptake of nutrients we have three maximum uptake rates: uptake from the limiting-nutrient pool by the resource , uptake from resource by the consumer , and uptake from consumer by the predator , with an assimilation efficiency , and three half saturation densities: the resource on the limiting nutrient pool , the consumer on the resource 0 , and the predator on the consumer 0 . Nutrients are lost from compartment according to .

C. P-C1-C2-R-N Model
We extended the P-C-R-N model to a five-compartment nutrient-limited model, adding a second consumer with a type II functional response. The specific equations are as follows: where is limiting-nutrient pool nutrients, is nutrients assimilated by the resource, 1 is nutrients assimilated by the first consumer, 2 is nutrients assimilated by the second consumer and is nutrients assimilated by the predator. The limiting-nutrient pool has external inputs . To model the uptake of nutrients we have four maximum uptake rates: uptake from the limiting-(3.2) (3.2) nutrient pool by the resource , uptake from resource by the first consumer 1 , uptake from resource by the second consumer 2 , and uptake from the first consumer by the predator 1 , with an assimilation efficiency , and three half saturation densities: the resource on the limiting nutrient pool , consumer on the resource 0 , and the predator on the first consumer 0 . Nutrients are lost from compartment according to .

D. Stability Analysis
Global stability was determined over a range of values for (Fig. 5). The lower endpoint of the range for is defined as the point where all compartments persist at densities greater than zero, with the upper endpoint for set to 5. The analyses were run for 10 000 integration steps, with the local maxima and minima of the consumer density collected over the last 100 integration steps. The analyses were only performed on the last 100 integration steps to ensure steady state.
The local stability of the interior equilibrium was analyzed over a range of values for identical to global stability analyses (Fig. 5). The interior equilibrium, Jacobian matrix, and eigenvalues were calculated numerically using Mathematica 11.
The parameter values for To determine the generality of the stability relationship between the R-N and C-R module we perform the experiment of increasing 100 times for the C-R-N system. We randomized all parameters within uniform ranges and checked for feasibility across the range of values (0.1-10). The upper limit of the range was increased from 5 to 10 as changes to the parameterization of the system are relative and can shift the quantitative points of stabilization. Thus, a wider range of tested values allows us to confirm that a transition to increasing stability (as observed in Fig. 5a where an initial increase is destabilizing, before the system is stabilized) that is shifted right is still consistent with the qualitative pattern of stabilization (i.e., as the interaction strength is increased, any instability that initially emerges is dampened by stronger interaction strengths). We repeated the randomization process until 100 feasible parameter sets were obtained. We then performed the experiment for each parameter set, using visual confirmation to ensure that the qualitative pattern described in Fig. 5 (i.e., strong R-N interactions correspond to cycles of decreased amplitude or a stable equilibrium is reached, when compared to weak R-N interactions) remained for all cases.
Note that was held constant at the same value, as the point of the experiment is to determine the ability of a strong R-N interaction to dampen out the oscillatory potential of a strong C-R interaction. Parameters were randomized within the range of the following distributions: In all cases, strong R-N interactions corresponded to cycles of decreased amplitude when compared to weak R-N interactions, or cycles were damped out to a stable equilibrium as interaction strength was increased.

E. Nutrient Loading Stability Analysis
Global stability was determined over a range of values for ( Supplementary Fig. 2). The lower endpoint of the range for is defined as the point where all compartments persist at densities greater than zero, with the upper endpoint for set to 3. The analyses were run for 10 000 integration steps, with the local maxima and minima of the consumer density collected over the last 100 integration steps. The analyses were only performed on the last 100 integration steps to ensure steady state.
The local stability of the interior equilibrium was analyzed over a range of values for identical to global stability analyses. The interior equilibrium, Jacobian matrix, and eigenvalues were calculated numerically using Mathematica 11.
The parameter values for Supplementary Fig. 2

A. C-R-N-D Model
We extended the C-R-N model to a four-compartment nutrient-limited ecosystem model, adding a detrital compartment that mineralizes nutrients back to the limiting-nutrient pool. The specific equations are as follows: where is limiting-nutrient pool nutrients, is nutrients assimilated by the resource, is nutrients assimilated by the consumer, and is nutrients present in detritus. The limiting-nutrient pool has an external input and an internal recycling input where is the rate of mineralization of nutrients in the detrital compartment. To model the uptake of nutrients we have two maximum uptake rates: uptake from the limiting-nutrient pool by the resource and uptake from resource by the consumer , with an assimilation efficiency , and two half saturation densities: the resource on the limiting nutrient pool and the consumer on the resource 0 . Nutrients are lost from compartment according to . Nutrients from compartment C and R are conserved in the system and contribute to detritus according to the mortality term . Inefficiencies in nutrient assimilation and sloppy feeding by consumers also contribute towards the detrital pool.

B. Stability Analysis
Global stability was determined over a range of values for (Fig. 6). The lower endpoint of the range for is defined as the point where all compartments persist at densities greater than zero, with the upper endpoint for set to 5. The analyses were run for 10 000 integration steps, with the local maxima and minima of the consumer density collected over the last 100 integration steps. The analyses were only performed on the last 100 integration steps to ensure steady state. (4.0) The local stability of the interior equilibrium was analyzed over a range of values for identical to global stability analyses. The interior equilibrium, Jacobian matrix, and eigenvalues were calculated numerically using Mathematica 11.

C. Nutrient Loading Stability Analysis
Global stability was determined over a range of values for ( Supplementary Fig. 3). The lower endpoint of the range for is defined as the point where all compartments persist at densities greater than zero, with the upper endpoint for set to 3. The analyses were run for 10 000 integration steps, with the local maxima and minima of the consumer density collected over the last 100 integration steps. The analyses were only performed on the last 100 integration steps to ensure steady state.
The local stability of the interior equilibrium was analyzed over a range of values for identical to global stability analyses. The interior equilibrium, Jacobian matrix, and eigenvalues were calculated numerically using Mathematica 11.