Strong nutrient-plant interactions enhance the stability of ecosystems
Review of C–R stability theoryTo set the context for how the R–N module will be used to understand the dynamics of nutrient-limited ecosystem models, we first briefly review stability results from modular food web theory. We do this by laying out a set of examples that serve to illustrate that in general, strong C–R interactions promote oscillatory dynamics while carefully placed weak C–R interactions dampen them5. We begin with the Rosenzweig–MacArthur C–R system as our base C–R module (Fig. 1a). It is biologically supported and produces a range of biologically plausible dynamics5, making it an appropriate system for this analysis. It exhibits three different dynamical phases over a gradient of interaction strengths (energetically defined sensu Nilsson et al. 2018) such that increasing the attack rate (({a}_{{CR}})) increases interaction strength15 (Fig. 2). We use the return time after a small perturbation (i.e., eigenvalues) to highlight the natural stability trade-off that occurs as interaction strength is changed, (i.e., the “checkmark” stability pattern)5,6. Equations and parameters can be found in Supplementary Results 1A. We draw your attention to three notable dynamical phases of the C–R module. At low interaction strengths the dominant eigenvalue (({lambda }_{{max }})) is negative and real and the C–R module follows a monotonic return to a stable equilibrium (Fig. 2a). During this phase ({lambda }_{{max }}) decreases from 0 (i.e., where ({a}_{{{CR}}}) allows the consumer to persist) to more negative values and thus stronger interactions tend to increase stability (Fig. 2d, i). At moderate interaction strengths, there is a sudden shift to eigenvalues with a non-zero complex part and population dynamics overshoot the equilibrium (Fig. 2b). Increases in interaction strength then further excite population dynamics and we observe less stable dynamics across this phase (Fig. 2d, ii). Last, the system reaches a Hopf bifurcation where the dominant eigenvalue becomes positive, yielding sustained cycles or oscillations (Fig. 2c, d, iii). As interaction strength increases across this phase, it is difficult to determine stability from the magnitude of a positive dominant eigenvalue; however, destabilization with increased interaction strength is readily observed in that the cycles become increasingly larger oscillations with a high coefficient of variation (CV)5. Note that while the Rosenzweig–MacArthur C–R system is shown here under a single set of parameters, analysis of the Jacobian shows the qualitative results to be general5. Moreover, the qualitative stability pattern remains for a type I and type III functional response5.Fig. 1: C–R and R–N base modules.a Rosenzweig–MacArthur C–R module modelled with Holling type II functional response and logistic resource growth, where (R) is resource biomass and (C) is consumer biomass. Parameters: (r) is the intrinsic growth rate of (R), (K) is the carrying capacity of (R), ({a}_{{mathrm {CR}}}) is the attack rate of (C) on (R), (e) is the assimilation rate of (C), ({R}_{0}) is the half-saturation density of (C), ({m}_{R}) and ({m}_{C}) are the mortality rates of (R) and (C), respectively. b R–N module modelled with a Monod nutrient uptake equation and external nutrient input, where (N) is a limiting-nutrient pool and (R) is the resource biomass. Parameters: ({I}_{N}) is external nutrient input to (N), ({a}_{{RN}}) is nutrient uptake rate by (R), (k) is the half-saturation density of (R), ({l}_{N}) and ({l}_{R}) are nutrient loss rates from (N) and (R), respectively.Full size imageFig. 2: C–R checkmark stability response.d Local stability (real and complex parts of the dominant eigenvalue; ({lambda }_{{max }})) as a function of interaction strength (({a}_{{{mathrm {CR}}}})) for the Rosenzweig–MacArthur C–R module. Time series reflect dynamics associated with region i, ii, and iii, respectively, following a perturbation that removes 50% of consumer biomass: a Stable equilibrium; monotonic dynamics. b Stable equilibrium; overshoot dynamics. c Unstable equilibrium; limit cycle. Boldness of arrows indicates the strength of interaction (({a}_{{CR}})).Full size imageWe now couple C–R modules into higher order food web modules to demonstrate how the addition of weak and/or strong interactions to a system can be used to predict dynamics at steady state (Fig. 3), constituting the “algebra” of C–R modules. Equations and parameters can be found in Supplementary Results 1B–D. We start with the three trophic level food chain (Fig. 3a), consisting of two coupled C–R modules (i.e., C1-R and P–C1). Theory has tended to find two weakly interacting C–R modules to generally produce locally stable equilibria16 (Fig. 3a). Increasing the strength of the C1–R interaction causes it to act like an oscillator (see Fig. 2c, above), and with enough increase this underlying oscillation is reflected in the limit cycles of the entire food chain (Fig. 3b). If the P–C1 interaction is strengthened as well, we end up with two coupled oscillators—the recipe for chaos17,18 (Fig. 3c). As such, coupled strong interactions are not surprisingly the recipe for complex and highly unstable dynamics.Fig. 3: Algebra of C–R modules.Time series showing the general dynamical outcomes for the food chain and diamond module at steady state with varied combinations of C–R interaction strengths. a Weak–weak interaction; point attractor. b Strong–weak interaction; limit cycle. c Strong–strong interaction; chaos. d Strong–strong, weak interaction; limit cycle. e Strong–strong, weak–weak interaction; point attractor.Full size imageFollowing McCann et al.19, we now add a weakly coupled consumer C2 to the food chain system of Fig. 3c. This weak consumer essentially draws energy away from the strong P–C1–R pathway and in doing so partially mutes the coupled oscillators, bringing the dynamics back to a more even limit cycle (Fig. 3d) and under certain conditions can drive equilibrium dynamics19. Last, the predator is weakly coupled to C2, creating a strong and weak pathway. The second weak interaction further draws energy away from the strong pathway, muting the oscillators entirely and bringing the system in this example to a point attractor (Fig. 3e). These examples show that well placed weak interactions (i.e., non-oscillatory phases, Fig. 2a, b) can be used to draw energy away from strong pathways and act as potent stabilizers of potentially oscillatory pathways. Note that weak interactions play a similarly stabilizing role in the omnivory module20 and further, weak interactions have been shown to stabilize large food web networks4,6 suggesting the principles derived from modular theory scale up to whole systems. Taken altogether, the oscillatory nature of strong C–R interactions generally promotes oscillatory dynamics in higher order systems, while the careful placement of weak C–R interactions—which are monotonic in nature—act to dampen oscillations. Although not discussed to our knowledge, we conjecture that if a subsystem exists such that strong interactions lead to monotonic dynamics (i.e., without oscillatory decay), strong interactions in this case would serve as a potent stabilizer. Below, we show the R–N module appears to be such a case.R–N module and stabilityTowards understanding how the R–N subsystem may interact in a higher order system, we first briefly consider the stability of the R–N module alone (akin to what we discussed for the C–R module above). The R–N module consists of a resource that takes up nutrients according to a Monod-like growth term, is open to flows from the external environment as a result of geochemical processes, and nutrients are lost to the external environment according to a linear term11 (Fig. 1b). Performing a local stability analysis about the interior equilibrium reveals the R–N module to be locally stable for all biologically feasible parameterizations, as determined by the signs of the trace and determinant of the Jacobian matrix (see Supplementary Results 2B). We now perform further numerical and analytical analyses to understand how stability is influenced by interaction strength.As the maximum rate of nutrient uptake (({a}_{{RN}})) is increased (i.e., R–N interaction strength), stability is generally increased (Fig. 4d), with the real part of the dominant eigenvalue (({lambda }_{{max }})) tending from 0 (i.e., where ({a}_{{RN}}) allows the resource to persist) towards an asymptote of ({-l}_{R}) (see Supplementary Results 2C). Numerical analysis reveals that the asymptote at ({-l}_{R}) can be approached from above or below depending on the relative leakiness of the R and N compartments (i.e., the rate at which nutrients are lost to the external environment from compartment R (({l}_{R})) and N (({l}_{N}))). For ({l}_{N} , > , {l}_{R}) (Fig. 4d), the R–N module only follows a monotonic return to equilibrium as interaction strength is increased, with increased interaction strength only tending to increased stability (i.e., reduce return time). For ({l}_{N} < {l}_{R}) (Fig. 4d), the R–N module follows a monotonic return to equilibrium for weak (Fig. 4a) and strong (Fig. 4c) interaction strength, but modest overshoot dynamics are observed for intermediate interaction strength (Fig. 4b). Stability tended to increase with interaction strength for weak to intermediate interaction strength (i.e., dominant eigenvalue becomes more negative), then slightly decrease as interaction strength became strong. A special case exists when ({l}_{R}={l}_{N}) (Fig. 4d), where stability increases with interaction strength until ({lambda }_{{max }}) becomes locked in at ({-l}_{R}), indicating stability does not change regardless of any further increase in interaction strength. Overall, the R–N interaction tends to generally stabilize in all cases (dominant eigenvalue goes from zero to a more negative saturating value with monotonic dynamics), although there are some intermediate cases that produce complex eigenvalues that suggest population dynamic overshoot potential (Fig. 4b). Note that we obtain qualitatively similar results when implicitly strengthening the R–N interaction by increasing nutrient loading (see Supplementary Results 2D and Supplementary Fig. 1). Now, given the above framework for coupled C–R modules—where weak C–R interactions with underlying monotonic dynamics dampen the oscillatory potential of strong C–R interactions—the underlying monotonic dynamics of the R–N module suggest that R–N interactions ought to be stabilizing when coupled to strong C–R interactions. Further, the underlying increase in stability (i.e., more rapid return to equilibrium) as R–N interaction strength is increased suggests the stabilizing potential of the R–N module ought to increase as the interaction becomes stronger.Fig. 4: R–N stability response to increasing interaction strength.Time series showing R density following a perturbation that lowered R density to 50% of equilibrium density for a low (({a}_{{RN}}=0.8)), b intermediate (({a}_{{RN}}=1)), and c high maximum rate of nutrient uptake (({a}_{{RN}}=2.8)). d Local stability (dominant eigenvalue; ({lambda }_{{max }})) of the R–N subsystem as ({a}_{{RN}}) is increased for ({l}_{N} , > , {l}_{R}), ({l}_{N}={l}_{R}), and ({l}_{N} < {l}_{R}), where ({l}_{R}) and ({l}_{N}) are the rate at which nutrients are lost to the external environment from compartment R and N, respectively. Solid lines are real parts and dashed lines are complex parts of ({lambda }_{{max }}).Full size imageTo look into this conjecture, we first coupled R–N to multiple configurations of strong and expectantly oscillatory C–R interactions and increased R–N interaction strength (({a}_{{RN}})). Following this, we added nutrient cycling and repeated the experiment to demonstrate that our results can be generalized to nutrient-limited ecosystem models. The full equations and parameter values for each model are listed in Supplementary Results 3A–D and 4A, B. We begin with the C–R–N system, where C–R and R–N are coupled through R (Fig. 5a). The initial increase in ({a}_{{RN}}) implicitly strengthens the C–R interaction and fuels the oscillatory potential of C–R and cycles emerge almost immediately after C is able to persist. As ({a}_{{RN}}) is increased further the cycles disappear and we obverse a steep stabilization phase, followed by a modest period of destabilization. Adding a weakly coupled predator gives a similar outcome, although the system continually stabilizes as ({a}_{{RN}}) is increased (Fig. 5b). If the P–C interaction is strengthened (i.e., both C1–R and P–C1 are strong, the recipe for chaos), R–N is unable to dampen oscillations even with a strong interaction strength, although a strong interaction gives tighter bound cycles than a weak interaction (Fig. 5c). We next add a weakly coupled consumer to the nutrient-limited food chain with strong P–C1 and C1–R interactions (Fig. 5d). As seen previously, this interaction draws energy out of the strong pathway, partially muting oscillatory potential. Thus, the ability for a strong R–N interaction to once again return the system to a stable equilibrium is not surprising. Finally, we add a detrital compartment to show that strong R–N interactions remain potent stabilizers in the context of nutrient cycling (Fig. 6b) when compared to a nutrient-limited food chain without nutrient cycling (Fig. 6a).Fig. 5: Nutrient-limited food chain stability.a–d Non-equilibrium dynamics (log10(C1,max/C1,min)) and equilibrium stability (real part of the dominant eigenvalue; ({lambda }_{{max }})) of the C–R–N, P–C–R–N with a single oscillator, P–C–R–N with coupled oscillators, and P–C1–C2–R–N modules, respectively, as ({a}_{{RN}}) is varied.Full size imageFig. 6: Nutrient-limited ecosystem module stability.a, b Non-equilibrium dynamics (log10(Cmax/Cmin)) and equilibrium stability (real part of the dominant eigenvalue; ({lambda }_{{max }})) of the C–R–N nutrient-limited food chain model and the C–R–N–D nutrient-limited ecosystem model, respectively, as ({a}_{{RN}}) is varied.Full size imageNote that we repeat our analysis of higher order modules by implicitly increasing R–N interaction strength through nutrient loading (see Supplementary Results 3E and 4C and Supplementary Figs. 2 and 3). In all cases, increased nutrient loading led to less stable dynamics, consistent with DeAngelis’ (1992) paradox of enrichment finding where increased nutrient loading lead to destabilizing autotroph–herbivore oscillations. More
