The effects of excess food nutrient content on a tritrophic food chain model in the aquatic ecosystem

The Effects of Excess Food Nutrient Content on a Tritrophic Food Chain Model in the Aquatic Ecosystem

Journal Pre-proof

The Effects of Excess Food Nutrient Content on a Tritrophic Food Chain Model in the Aquatic Ecosystem Lale Asik, Ming Chen, Angela Peace PII: DOI: Reference:

S0022-5193(20)30039-4 https://doi.org/10.1016/j.jtbi.2020.110183 YJTBI 110183

To appear in:

Journal of Theoretical Biology

Received date: Revised date: Accepted date:

8 May 2019 28 January 2020 31 January 2020

Please cite this article as: Lale Asik, Ming Chen, Angela Peace, The Effects of Excess Food Nutrient Content on a Tritrophic Food Chain Model in the Aquatic Ecosystem, Journal of Theoretical Biology (2020), doi: https://doi.org/10.1016/j.jtbi.2020.110183

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Ltd.

1

2

3 4

5

Highlights • Nutrient-rich food consumption causes (direct or indirect) adverse effects on species. • The tri-trophic stoichiometric knife-edge model exhibits periodic oscillations and chaotic behaviors. • The tri-trophic stoichiometric knife-edge model captures several types of species interactions.

1

7

The Effects of Excess Food Nutrient Content on a Tritrophic Food Chain Model in the Aquatic Ecosystem

8

Lale Asik1a , Ming Chenb , Angela Peacea

6

a

9 10

11

Department of Mathematics and Statistics, Texas Tech University, Lubbock, USA b School of Science, Dalian Maritime University, Dalian, Liaoning, China

Abstract Ecological stoichiometry is an approach that focuses on the balance of energy and elements in environmental interactions, and it leads to new insights and a better understanding of ecological processes and outcomes. Modeling under this framework enables us to investigate the effects of nutrient content (i.e., food quality) on organisms, whether the imbalance involves insufficient or excess nutrient content. In this paper, we develop and analyze a tritrophic food chain model that captures the phenomenon known as the “stoichiometric knife-edge”, where consumer growth is limited under conditions of excess nutrients. The model tracks two essential elements, carbon and phosphorus, in each species. The dynamics of the system such as boundedness and positivity of the solutions, existence and stability conditions of boundary and internal equilibria are analyzed. Through numerical simulations and bifurcation analyses, we observe the dynamics of the system switching between periodic oscillations and chaos. Our findings also show that nutrient-rich food consumption can cause adverse effects on species.

13

Keywords: Ecological stoichiometry, stoichiometric knife-edge, tritrophic food chain model, chaos, maximum Lyapunov exponent

14

1. Introduction

12

15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

Ecological stoichiometry is a branch of ecology that focuses on the effects of the balance of energy and elements on organisms and their interactions in ecosystems based on the laws of conservation of matter and energy (Sterner and Elser 2002). In the last two decades, multiple stoichiometrically-explicit models of producer–consumer interactions have been developed (Sterner, 1990; Andersen, 1997; Muller et al., 2001; Mulder and Bowden, 2007), which have provided insights into the coupling of population dynamics and nutrient recycling. Loladze et al. (2000) extended the stoichiometric concept by coupling chemical heterogeneity with Lotka-Volterra type equations to capture the effects of food quality and nutrient recycling feedbacks. Qualitative analysis of this model revealed a new phenomenon, the paradox of energy enrichment, where intense energy enrichment substantially elevates producer density but, despite such an abundant food supply, consumer decreases its growth rate and drives itself to deterministic extinction. While this model makes the assumption that nutrients in the dead and excreted matter is immediately recycled and acquired by the producer, Diehl (2007) showed that relaxing this assumption does not change the major features of the model outcomes. Recently, stoichiometric models are being applied to examine a phenomenon termed the “stoichiometric knife-edge”, where consumer dynamics are affected by both insufficient as well as excess food nutrient content (Elser et al., 2012; Peace et al., 2013, 2014; Asik and Peace, 2019). Common 1 corresponding author at: Tel: +1 806 742 1112; E-mail address: [email protected] Preprint submitted to Journal of Theoretical Biology

February 7, 2020

66

to these models is the (empirically supported-Plath and Boersma (2001)) assumption that the consumer ingests nutrients up to the rate required for its maximal growth but not more. This assumption leads to another paradox, the paradox of nutrient enrichment: excessive nutrient supply can inflate producer density but depress consumer density, eventually causing extinction. Similar to the findings of the knife-edge models, others have shown that nutrient-rich food can also cause consumer extinction (Branco et al., 2018). All existing models which incorporate the “stoichiometric knife-edge phenomenon” have focused on two species. They have been applied to autotrophic phytoplankton (algal) and herbivorous zooplankton (Daphnia) systems in aquatic environments. However, food web models involving only two species may miss important ecological behaviors. Non-stoichiometric continuous time models of a two species producer–consumer ecosystem have only two basic patterns: approach to an equilibrium or to a limit cycle (e.g., Rosenzweig and MacArthur 1963). The dynamic behavior of stoichiometric models (in continuous time) of two-species can show additional patterns: the occurrence of Hopf and saddle-node bifurcations (e.g., Loladze et al. 2000; Peace et al. 2013). However, in recent decades, it has been demonstrated that more complex dynamic behaviors can appear in continuous time models with three or more species. For example, Hastings and Powell (1991) investigated a three-species food chain model without stoichiometry and demonstrated that, for biologically reasonable choices of parameters, chaotic dynamics result. Chen et al. (2017) studied a three-dimensional model incorporating two resources to show that a resource competition model can generate oscillations and chaos when species compete for limiting resources. In this study, we extend the bitrophic stoichiometric knife-edge model (Peace et al., 2013) to the tritrophic case to obtain more realistic food web interactions. The third tropic level is composed of a larval gizzard shad (fish) which feeds only on Daphnia. The model tracks two essential elements, carbon (as a surrogate for biomass) and phosphorus (a limiting nutrient in many aquatic ecosystems), in each species. As far as we know, this is the first study that investigates the effects of excess nutrient content in tritrophic food chains using mathematical models. Our model exhibits nonintuitive paradoxes, (paradox of energy enrichment, paradox of nutrient enrichment). A closer look at the mechanisms of interactions reveals that variable stoichiometry of algal has many effects, direct and indirect, positive and negative, not only on Daphnia but also on lavral gizzard shad. Considering the third trophic level yields the existence of chaotic behaviour under biologically reasonable parameter ranges. In the next section, we provide the details of the model construction. In Section (3), we analyze it qualitatively. In Section (4), we construct numerical simulations and bifurcation diagrams of the system. The last section summarizes the main results and suggests directions for future research.

67

2. Model Construction

33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

68

We start with the two dimensional stoichiometric knife-edge model studied by Peace et al. (2013)

69

( )   fˆθy dx x = bx 1 − − min f (x), y dt min{K, (P − θy y)/q} Q ( )   fˆθy dy Q = min eˆy , min f (x), y − dy y dt θy Q 70 71

(1a) (1b)

where Q = (P − θy y)/x describes the variable phosphorus to carbon (P:C) ratio of the producer. x(t) and y(t) are the biomass of the producer and consumer, respectively, measured in terms of C. 3

72 73 74 75 76 77

b is the maximum growth rate of producer. K is the producer carrying capacity (due to external constraints of C, such as light limitation). P is the total phosphorus in the system. θy is the consumer’s constant P:C ratio. q is the minimal P:C ratio in producer. eˆy is the maximal production efficiency of consumer and the second law of thermodynamics requires that eˆy < 1. dy is the the consumer loss rate. f (x) is the consumer’s ingestion rate, which is taken here as a Holling type II functional response. In general, the function f (x) is a bounded differentiable function that satisfies: f (0) = 0, f 0 (x) > 0, f 0 (0) < ∞ and f 00 (x) < 0 for x ≥ 0.

78

79 80

81 82

83 84

85 86

87 88 89 90 91 92 93 94

95 96 97 98 99 100

101 102 103

104 105

(2)

f (x) is saturating with limx→∞ f (x) = fˆ. The model makes the following four assumptions. Assumption 1. The total amount of phosphorus in the entire system is fixed, i.e., the system is closed for phosphorus with a total of P (mg P/L). Assumption 2. P : C ratio in the producer varies, but it never falls below a minimum value q (mg P/ mg C); the consumer maintains a constant P : C, θy (mg P/ mg C). Assumption 3. All phosphorus in the system is divided into two pools: phosphorus in the producer and phosphorus in the consumer. Assumption 4. The consumer ingests phosphorus up to the rate required for its maximal growth but not more. These assumptions introduce a new constraint in terms of C and P. Three minimum functions have been obtained, according to Liebig’s Law of Minimum. Here, the first minimum function, min{K, (P − θy y)/q}, is used to describe carrying capacity of the producer determined by light (or C) and P availability. The second minimum function, min{f (x), fˆθy /Q}, is used to describe the consumer’s ingestion rate. The first input, f (x), is the consumer’s ingestion rate when P is not in excess and the second input, fˆθy /Q, is the consumer’s ingestion rate when P is in excess. The third minimum function, min{ˆ ey , Q/θy }, is used to describe the production efficiency of consumer, which depends on both energetic and nutrient limitations. Notice that, since eˆy f (x) < fˆ, one obtains that ) (     fˆθy θy Q Q ˆ min eˆy , min f (x), = min eˆy f (x), f (x), eˆy f . θy Q θy Q Biologically, this translates into three cases growth is determined by energy  in which consumer   θy Q ˆ limitation (ˆ ey f (x)), P limitation θy f (x) , and P in excess eˆy f Q . In order to investigate the effects of excess nutrient on systems of three trophic levels, we expanded the the bitrophic model (System (1)) by following the procedure used by Peace (2015). A predator (z) with a fixed P:C ratio (θz ) is added into the system. As a result of this, the second and third assumptions should be modified: Assumption 5. The P : C ratio of the producer varies, but it never falls below a minimum q (mg P/ mg C); the consumer and the predator maintain constant P : C ratios, θy and θz (mg P/ mg C), respectively. Assumption 6. All phosphorus in the system is divided into three pools: phosphorus in the producer, phosphorus in the consumer, and phosphorus in the predator.

4

106 107

108 109 110 111 112

From these modified assumptions, it follows that P available for the producer at any given time is P − θy y − θz z (mg P/L). The tritrophic stoichiometric knife-edge model takes the following form: ( )   fˆθy dx x − min f (x), y (3a) = bx 1 − dt min{K, (P − θy y − θz z)/q} Q   θy dy Q ˆ y − g(y)z − dy y (3b) = min eˆy f (x), f (x), eˆy f dt θy Q   θy dz g(y)z − dz z (3c) = min eˆz , dt θz where Q = (P − θy y − θz z)/x is the variable P:C ratio of the producer. eˆz is the predator maximal production efficiency in the terms of C which is less than 1 due to the second law of thermodynamics. dz is the predator loss rate. θz is the predator’s constant P:C ratio. The predator’s ingestion rate, g(y), is assumed to have a similar form as the consumer’s ingestion rate, f (x), i.e., it satisfies the following assumptions: g(0) = 0, g 0 (y) > 0, g 0 (0) < ∞ and g 00 (y) < 0 for y ≥ 0.

(4)

116

g(y) is saturating with limy→∞ g(y) = gˆ. Here, minimum function, min{ˆ ez , θy /θz }, is used to describe the production efficiency of the predator, which depends on both energetic and nutrient limitations. It is important to note that this tritrophic model is simply an expansion of the two species stoichiometric knife-edge model (1) to include three species in a linear food chain.

117

3. Model Analysis

113 114 115

118 119

The following theorem shows that System (3) is biologically well defined and that solutions are unique.

120

Theorem 3.1. The model, given in System (3) is well defined as x → 0.

121

Proof. x0 (t) is well defined at x → 0, since

( )   fˆθy x dx = bx 1 − − min f (x), y dt min{K, (P − θy y − θz z)/q} Q ( )   fˆθy x x = bx 1 − − min f (x), y. min{K, (P − θy y − θz z)/q} P − θy y − θ z z

122

123

124

From condition (2) it follows that for f (x)/x the following holds:   f (x) f (x) 0 < 0 for x > 0 and lim = f 0 (0) < ∞. x→0 x x Thus, y 0 (t) is well defined at x → 0, since  θ  eˆy f (x)y − g(y)z − dy y, if eˆy f (x) < θQy f (x), eˆy f (x) < eˆy fˆ Qy   dy (P −θy y−θz z) θ = θQy f (x)y − g(y)z − dy y = f (x) − g(y)z − dy y, if θQy f (x) < eˆy f (x), θQy f (x) < eˆy fˆ Qy x θy  dt  θ θ eˆy fˆθy y − g(y)z − dy y = eˆy fˆ θy x if eˆy fˆ Qy < eˆy f (x), eˆy fˆ Qy < θQy f (x). Q P −θy y−θz z y − g(y)z − dy y, 5

126

The boundedness and positive invariance of the solutions of (3) are assured by the following theorem.

127

Theorem 3.2. Let k = min{K, P/q}. Then, solutions with initial conditions in

125

Ω = {(x, y, z) : 0 < x < k, 0 < y < P/θy , 0 < z < P/θz , qx + θy y + θz z < P } 128

129 130 131 132 133 134

remain there for all forward times. Proof. We consider a solution S(t) ≡ (x(t), y(t), z(t)) of system (3) with initial conditions in Ω. Hence, 0 < x(0) < k = min{K, P/q}, 0 < y(0) < P/θy , 0 < z(0) < P/θz , qx(0)+θy y(0)+θz z(0) < P . ¯ (closure of Assume that there exists a time t1 > 0 such that S(t) touches or crosses a boundary of Ω Ω) for the first time, and then (x(t), y(t), z(t)) ∈ Ω for 0 ≤ t < t1 . We have several cases to consider. P Case 1. x(t1 ) = 0. Let f = f 0 (0) = limx→0 f (x) x , y = maxt∈[0,t1 ] y(t) < θy , and z = maxt∈[0,t1 ] z(t) < P θz .

Then, for every t ∈ [0, t1 ], we have  dx = bx 1 − dt min{K, (P  = bx 1 − min{K, (P "  = b 1− min{K, (P "  ≥ b 1− min{K, (P "  ≥ b 1− min{K, (P

135 136 137

138 139 140

(5)

x − θy y − θz z)/q}



(

fˆθy − min f (x), Q ( 

)

y

) fˆθy x − min f (x), y P − θy y − θz z ( ) #  fˆθy f (x) x − min , y x − θy y − θz z)/q} x P − θy y − θz z ( ) #  ˆθy f x − min f¯, y x − θy y − θz z)/q} P − θy y − θz z ( ) #  ˆθy f k − min f¯, y¯ x ≡ α1 x − θy y¯ − θz z¯)/q} P − θy y¯ − θz z¯

x − θy y − θz z)/q}

where α1 is a constant. Thus x(t) ≥ x(0)eα1 t for t ∈ [0, t1 ], which implies that x(t1 ) ≥ x(0)eα1 t1 > 0, a contradiction. This proves that S(t1 ) does not reach this boundary. Case 2. x(t1 ) = k = min{K, P/q}. ( )   fˆθy dx x = bx 1 − − min f (x), y dt min{K, (P − θy y − θz z)/q} Q   x ≤ bx 1 − min{K, (P − θy y − θz z)/q}    x x ≤ bx 1 − = bx 1 − . min{K, P/q} k Then, the standard comparison argument yields that x(t) < k for all t ≥ 0. Thus S(t1 ) does not reach this boundary. P Case 3. y(t1 ) = 0. Let g = g 0 (0) = limy→0 g(y) y and z = maxt∈[0,t1 ] z(t) < θz . Then, for every

6

141

t ∈ [0, t1 ], we have

  θ dy Q y = min eˆy f (x), f (x), eˆy fˆ y − g(y)z − dy y dt θy Q     θy Q g(y) ˆ = min eˆy f (x), f (x), eˆy f − z − dy y θy Q y     θy Q ˆ − g¯z¯ − dy y ≥ min eˆy f (x), f (x), eˆy f θy Q ≥ −(¯ g z¯ + dy )y ≡ α2 y.

142 143 144

145 146 147 148 149

where α2 is constant. This implies that y(t1 ) ≥ y(0)eα2 t1 > 0, a contradiction. This proves that S(t1 ) does not reach this boundary. Case 4. z(t1 ) = 0. For every t ∈ [0, t1 ], we have   θy dz = min eˆz , g(y)z − dz z ≥ −dz z. dt θz This implies that z(t1 ) ≥ z(0)e−dz t1 > 0. This contradicts z(t1 ) = 0 and proves that S(t1 ) does not reach this boundary. Case 5. Assume that qx(t1 ) + θy y(t1 ) + θz z(t1 ) = P and qx(t) + θy y(t) + θz z(t) < P for t ∈ [0, t1 ). Then x(t1 ) = (P − θy y(t1 ) − θz z(t1 ))/q and Q(t1 ) = (P − θy y(t1 ) − θz z(t1 ))/x(t1 ) = q. It is easy to see that qx0 (t1 ) + θy y 0 (t1 ) + θz z 0 (t1 ) ≥ 0. ( )   fˆθy dx x(t1 ) − min f (x(t1 )), y(t1 ) = bx(t1 ) 1 − dt t = t min {K, (P − θy y(t1 ) − θz z(t1 ))/q} Q(t1 ) 1 ( )   fˆθy x(t1 ) x(t1 ) − min f (x(t1 )), = bx(t1 ) 1 − y(t1 ) min {K, (P − θy y(t1 ) − θz z(t1 ))/q} P − θy y(t1 ) − θz z(t1 ) ( )   fˆθy x(t1 ) x(t1 ) − min f (x(t1 )), y(t1 ) ≤ bx(t1 ) 1 − (P − θy y(t1 ) − θz z(t1 ))/q q (P − θy y(t1 ) − θz z(t1 ))/q ( )   fˆθy x(t1 ) x(t1 ) = bx(t1 ) 1 − − min f (x(t1 )), y(t1 ) x(t1 ) q x(t1 ) ( ) fˆθy = − min f (x(t1 )), y(t1 ) q   θy Q(t1 ) dy ˆ = min e ˆ , f (x(t )), f (x(t )), e ˆ f y(t1 ) − g(y(t1 ))z(t1 ) − dy y(t1 ) y 1 1 y dt t = t θy Q(t1 ) 1 ( )   fˆθy Q(t1 ) = min eˆy , min f (x(t1 )), y(t1 ) − g(y(t1 ))z(t1 ) − dy y(t1 ) θy Q(t1 ) ( )   fˆθy q = min eˆy , min f (x(t1 )), y(t1 ) − g(y(t1 ))z(t1 ) − dy y(t1 ) θy q ( ) fˆθy q < min f (x(t1 )), y(t1 ) − g(y(t1 ))z(t1 ) θy q 7

  θy dz = min eˆz , g(y(t1 ))z(t1 ) − dz z(t1 ) dt t = t θz 1 θy < g(y(t1 ))z(t1 ) θz (

fˆθy qx0 (t1 ) + θy y 0 (t1 ) + θz z 0 (t1 ) < − q min f (x(t1 )), q

)

(

fˆθy y(t1 ) + q min f (x(t1 )), q

)

y(t1 )

− θy g(y(t1 ))z(t1 ) + θy g(y(t1 ))z(t1 ) = 0 150

151 152 153 154 155

156

This contradicts the assumption qx0 (t1 ) + θy y 0 (t1 ) + θz z 0 (t1 ) ≥ 0. ˘ we denote the region Ω with its entire boundary except of the biologically unrealistic By Ω, edges γ1 = {(x, y, z) : x = 0, P = θy y + θz z} and γ2 = {(x, y, z) : y = 0, P = qx + θz z}. (At the first edge consumer and predator contain all P, thus there is no producer in the system; and at the second one producer and predator contain all P, hence there is no consumer in the system.) In other words, ˘ = (∂Ω \ (γ1 ∪ γ2 )) ∪ Ω. Ω To simplify our analysis, we rewrite System (3) in the following form: dx = xF (x, y, z), dt

157

dy = yG(x, y, z), dt

dz = zH(x, y, z) dt

(6)

where ( )  fˆθy x f (x) F (x, y, z) = b 1 − − min , y min{K, (P − θy y − θz z)/q} x P − θy y − θ z z   P − θy y − θz z θy x g(y) ˆ f (x), eˆy f − z − dy G(x, y, z) = min eˆy f (x), θy x P − θy y − θ z z y   θy H(x, y, z) = min eˆz , g(y) − dz . θz 

158 159

Conditions (2) and(4) ensure that all partial derivatives of F, G, and H exist almost everywhere on ˘ Ω:

8

Fx

Fy

Fz

Gx

  0 fˆθ f (x) b  − − y, if f (x) < Qy ∂F x min{K,(P −θy y−θz z)/q} = = fˆθ b − ∂x if f (x) > Qy min{K,(P −θy y−θz z)/q} < 0,  fˆθ P −θy y−θz z f (x)  if f (x) < Qy , K < − x < 0, q    fˆθy fˆθy (P −θz z) P −θy y−θz z  < 0, if f (x) > − , K < ∂F Q q (P −θy y−θz z)2 = = ˆθy f bxqθ P −θ y−θz z f (x) y y  ∂y if f (x) < Q , K > − (P −θy y−θz z)2 − x < 0,   q    bxqθy P −θy y−θz z fˆθy (P −θz z) fˆθ − (P −θy y−θ if f (x) > Qy , K > 2 − (P −θ y−θ z)2 < 0, q z z) y z  P −θy y−θz z fˆθy  0, if f (x) < Q , K <  q    fˆθy yθz fˆθy P −θy y−θz z  − < 0, if f (x) > , K < ∂F Q q (P −θy y−θz z)2 = = fˆθy P −θy y−θz z bxqθz  ∂z − (P −θy y−θz z)2 < 0, if f (x) < Q , K >   q    P −θy y−θz z fˆθy yθz fˆθy bxqθz − (P −θy y−θz z)2 − (P −θy y−θz z)2 < 0, if f (x) > Q , K > q  θ 0 if eˆy f (x) < θQy f (x), eˆy f (x) < eˆy fˆ Qy  eˆy f (x) > 0,    0 ∂G θ f (x) zz = = P −θyθy−θ < 0, if θQy f (x) < eˆy f (x), θQy f (x) < eˆy fˆ Qy x y  ∂x  θy θ θ  ˆ eˆy f > 0, if eˆy fˆ y < eˆy f (x), eˆy fˆ y < Q f (x) P −θy y−θz z

Gy =

∂G ∂y

Gz =

∂G ∂z

Q

Q

(7)

θy

0   θ g(y)  − z > 0, if eˆy f (x) < θQy f (x), eˆy f (x) < eˆy fˆ Qy  y    0 θ g(y) = − f (x) − z, if θQy f (x) < eˆy f (x), θQy f (x) < eˆy fˆ Qy x y  0    θ θ g(y)  eˆy fˆθy2 x − z > 0, if eˆy fˆ Qy < eˆy f (x), eˆy fˆ Qy < θQy f (x) 2 y (P −θy y−θz z)  g(y) θ  < 0, if eˆy f (x) < θQy f (x), eˆy f (x) < eˆy fˆ Qy −   y θy f (x) g(y) Q Q = − θθyz x − y < 0, if θy f (x) < eˆy f (x), θy f (x) < eˆy fˆ Q    eˆy fˆθy θz x − g(y) , if eˆ fˆθy < eˆ f (x), eˆ fˆθy < Q f (x) y y y 2 (P −θy y−θz z)

y

Q

Q

θy

∂H =0 ∂x   θy ∂H Hy = = min eˆz , g 0 (y) > 0 ∂y θz ∂H Hz = = 0. ∂z

Hx =

160 161

3.1. Boundary equilibria To find equilibria we solve xF (x, y, z) = 0,

162 163

yG(x, y, z) = 0,

zH(x, y, z) = 0.

The boundary equilibria are E0 = (0, 0, 0), E1 = (k, 0, 0) = (min{K, P/q}, 0, 0), and E2 = (¯ x, y¯, 0). To determine the local stability of these equilibria we calculate the Jacobian of (6),   F (x, y, z) + xFx (x, y, z) xFy (x, y, z) xFz (x, y, z) . yGx (x, y, z) G(x, y, z) + yGy (x, y, z) yGz (x, y, z) J=  zHx (x, y, z) zHy (x, y, z) H(x, y, z) + zHz (x, y, z) 9

164

165

166

167 168

169

170

171

Theorem 3.3. The equilibrium point at the origin, E0 , is a saddle (always unstable). Proof. At the origin, the Jacobian matrix takes the form   b 0 0 0 . J(E0 ) = 0 −dy 0 0 −dz The eigenvalues have different signs. Therefore, E0 is a saddle. Theorem 3.4. The consumer and predator extinction equilibrium, E1 , is a locally asymptotically stable node if   θy k P < dy ; min eˆy f (k), f (k), eˆy fˆ θy k P otherwise, it is unstable (being of the saddle type). Proof. At E1 = (k, 0, 0) = (min{K, P/q}, 0, 0), the Jacobian matrix is   −b kFy (k, 0, 0) kFz (k, 0, 0)  G(k, 0, 0) 0 J= 0 0 0 −dz where 

θy k P G(k, 0, 0) = min eˆy f (k), f (k), eˆy fˆ θy k P 172 173 174

175 176 177 178 179 180 181 182



The stability of E1 depends on the sign n of G(k, 0, 0). If G(k, 0, o 0) is positive, then E1 is a saddle. If θy k P ˆ G(k, 0, 0) is negative, i.e., dy > min eˆy f (k), θy k f (k), eˆy f P , then E1 is a locally asymptotically stable node. The boundary equilibrium E2 = (¯ x, y¯, 0) indicates the extinction of the predator. We define the conditions for the existence of E2 and check its stability by following closely the terminologies in Loladze et al. (2000), Peace et al. (2013), and Chen et al. (2017). E2 can be viewed as the internal equilibrium of the degenerating simpler two-dimensional system only with one consumer on one producer. It is, therefore, sufficient to study the x − y phase plane (see Fig. 1). To investigate the interior equilibrium the x − y phase plane is divided into three biologically significant regions by ˆ fˆθ the two lines eˆy = θQy and f (x) = Qy . Peace et al. (2013) proved that if eˆy > θQy , then f (x) < fQθ . This shows the line f (x) = fˆθy Q .

fˆθy Q

lies beneath the line eˆy =

Q θy .

Region I is defined by eˆy <

183

f (x) <

184

(energy or C) limits consumer growth. Region II is defined by eˆy >

185 186

− dy .

Q θy

and

This represents the cases where P is neither limiting nor in excess, i.e., food quantity

limited by a deficiency of P. Region III is defined by eˆy < and reduces consumer growth.

10

Q θy

Q θy ;

here, consumer growth is

and f (x) >

fˆθy Q ,

where P is in excess

Figure 1: x − y phase plane for low total phosphorus P = 0.02 mg P/L. The three different regions are depicted: P is not deficient or in excess in Region I, P is limiting in Region II, and P is in excess in Region III. The solid lines represent consumer nullclines and the dashed lines represent producer nullclines. The open circles denote unstable equilibria and the filled in circles denote stable equilibria. 187

188 189

190 191 192 193 194 195 196

197

The following two inequalities provide sufficient criteria for the existence of E2 :   θy k P ˆ G(k, 0, 0) = min eˆy f (k), f (k), eˆy f − dy > 0 θy k P which shows that the consumer growth rate is larger than its death rate, ensuring the survival of the consumer and   θy H(¯ x, y¯max , 0) = min eˆz , g(¯ ymax ) − dz < 0 θz   d where y¯ ≤ y¯max = θPy − eˆy f −1 eˆyy . This reveals that even with sufficiently large biomass of the consumer, the predator death rate is greater than its growth rate, and hence the predator cannot survive. To analyze the local stability of E2 we apply the method of nullclines in Loladze et al. (2000). The knife-edge shaped consumer nullcline creates a possibility for multiple positive equilibria as Fig. 1 shows. Suppose that E2 = (¯ x, y¯, 0) is one such equilibria, i.e., F (¯ x, y¯, 0) = 0 = G(¯ x, y¯, 0). Then, the Jacobian matrix at E2 is   Jsub (E2 ) J1 J(E2 ) = , (8) 0 H(¯ x, y¯, 0) where Jsub (E2 ) =



 x ¯Fx (¯ x, y¯, 0) x ¯Fy (¯ x, y¯, 0) y¯Gx (¯ x, y¯, 0) y¯Gy (¯ x, y¯, 0) 11

198

and J1 =

199



 x ¯Fz (¯ x, y¯, 0) . y¯Gz (¯ x, y¯, 0)

The trace and the determinant of Jsub (E2 ) are Tr(Jsub (E2 )) = x ¯Fx + y¯Gy

(9)

Det(Jsub (E2 )) = x ¯y¯(Fx Gy − Fy Gx ).

(10)

200

201 202 203 204 205

206 207 208 209

210 211 212 213

We split the analysis into three cases depending on whether E2 is in Region I, II, or III. The slope of the producer and consumer nullclines at(x, y) are defined −Fx /Fyby and −Gx /Gy respectively. fˆθ

Case 1. Suppose that E2 lies in Region I eˆy < θQy and f (x) < Qy . Then (7) yields that at E2 , Fy < 0, Gx > 0, and Gy = 0. Therefore, the sign of the determinant (10) is positive, while   Fx sign(Tr(Jsub )) = sign(Fx ) = sign − . Fy If the producer nullcline is increasing at E2 (i.e., Fx > 0), then E2 is a saddle. If the producer nullcline is decreasing at E2 (i.e., Fx < 0), then E2 is locally asymptotically stable.  fˆθ

Case 2. Suppose that E2 lies in Region II eˆy > θQy and f (x) < Qy . Then (7) yields that at E2 , Fy < 0, Gx < 0, and Gy < 0. An elementary derivation indicates that      Fx Gy − Fy Gx Gx Fx sign(Det(Jsub )) = sign(Fx Gy − Fy Gx ) = sign = sign − − − . Fy Gy Gy Fy Hence, if the slope of consumer nullcline at E2 is less than the slope of the producer nullcline, i.e., −Gx /Gy < −Fx /Fy , then the determinant (10) is negative, yielding E2 as a saddle. If the consumer nullcline has a larger slope, then the determinant (10) is positive which makes the eigenvalues of the same sign as the trace (9). Recalling that Gx and Gy are negative, we find that 0 > −Gx /Gy > −Fx /Fy .

214 215 216 217

218 219 220 221 222 223

(11)

Condition (11) together with Fy < 0 implies that Fx < 0. Therefore, the trace (9) is negative, so E2 is locally asymptotically stable if the consumer nullcline has larger   slope than the producer’s. fˆθy Q Case 3. Suppose that E2 lies in Region III eˆy < θy and f (x) > Q . Then (7) yields that at E2 , Fx < 0, Fy < 0, Gx > 0, and Gy > 0. An elementary derivation indicates that      Fx Gy − Fy Gx Fx Gx sign(Det(Jsub )) = sign(Fx Gy − Fy Gx ) = sign = sign − − − . −Fy Gy Fy Gy If the slope of producer nullcline at E2 is less than the slope of consumer nullcline, i.e., −Gx /Gy > −Fx /Fy , then the determinant (10) is negative, yielding E2 as a saddle. If producer nullcline has larger slope, then the determinant (10) is positive. The stability of E2 depends on the sign of trace (9). Hence, if x ¯Fx + y¯Gy > 0, then the eigenvalues of the Jacobian matrix (8) have positive real parts and E2 is a saddle. If x ¯Fx + y¯Gy < 0, then the eigenvalues of the Jacobian matrix (8) are negative real parts and E2 is locally asymptotically stable.

224 225 226

In the following theorem, we summarize the existence and stability of the predator-extinction boundary equilibrium. 12

θ

d

233

Theorem 3.5. If H(¯ x, y¯max , 0) = min{ˆ ez , θyz }g(¯ ymax ) − dz < 0 where y¯max = θPy − eˆy f −1 ( eˆyy ), then the predator-extinction equilibrium E2 = (¯ x, y¯, 0) exists. In Region I, E2 is locally asymptotically stable if the producer nullcline is decreasing and unstable if it is increasing. In Region II, if the slope of consumer nullcline at E2 is greater than the producer’s, then E2 is locally asymptotically stable; otherwise, if the consumer nullcline has a smaller slope, it is a saddle. In Region III, if the slope of producer nullcline at E2 is less than the consumer’s, then E2 is a saddle; if producer nullcline has a larger slope, then its stability depends on the sign of x ¯Fx + y¯Gy .

234

3.2. Internal equilibria

227 228 229 230 231 232

235

If the following system has a solution F (x, y, z) = G(x, y, z) = H(x, y, z) = 0,

236 237 238

239

240 241 242 243

244 245 246 247

248 249 250 251 252

then system (6) has an internal (positive) equilibrium, which we denote as (x∗ , y ∗ , z ∗ ). To determine the type of species interactions occurring at this equilibrium, it is convenient to examine the Jacobian of system (6) at (x∗ , y ∗ , z ∗ ) which takes the following form:  ∗  x 0 0 J(x∗ , y ∗ , z ∗ ) =  0 y ∗ 0  E(x∗ , y ∗ , z ∗ ), 0 0 z∗ where



Fx E(x∗ , y ∗ , z ∗ ) = Gx Hx

Fy Gy Hy

 Fz Gz  Hz

(12)

is the ecosystem matrix of our system. In this matrix, the ij-th term measures the effect of the j-th species on the i-th species per-capita growth rate. ∗ If P is neither limiting nor in excess, i.e., Qeˆy > θy and f (x∗ )Q∗ < fˆθy , then by substituting signs in (7) into (12) we obtain   +/− − −  + + − , 0 + 0

which means that conventional consumer-producer, i.e., (+, −) type, interaction exists between the consumer and the producer. ∗ If the consumer growth is limited by a deficiency of P, i.e., Qeˆy < θy and f (x∗ )Q∗ < fˆθy , then by substituting signs in (7) into (12) we obtain   +/− − −  − +/− − , 0 + 0

which implies that the interaction between the consumer and the producer is shifted from conventional (+, −) type to competitive (−, −) type. The competition lies in two species not only internally but also interspecifically for the limited available P. ∗ If P is in excess, i.e., Qeˆy > θy and f (x∗ )Q∗ > fˆθy , then by substituting signs in (7) into (12) we obtain

13



− + 0

 − − + +/− , + 0

253

which shows that competitive consumer–producer is shifted to conventional consumer–producer.

254

3.3. Numerical analysis of internal equilibria

255 256 257

258 259 260 261 262 263 264

Here we provide a numerical analysis on the interior equilibria for varying values of total phosphorus P . We fix all other parameters with values listed in Table 1. If the system (3) has an internal equilibrium (x∗ , y ∗ , z ∗ ), from H ∗ (x∗ , y ∗ , z ∗ ) = 0, it follows that   d nz o . y ∗ = g −1  (13) θ min eˆz , θyz

Then the internal equilibrium can be studied by degenerating the system onto the 2-dimensional x − z plane by setting y = y ∗ . The internal equilibrium can be found at the intersection of the two curves F (x, y ∗ , z) = 0 and G(x, y ∗ , z) = 0, which are denoted as F0 (x, z) = 0 and G0 (x, z) = 0, in the x − z plane, respectively. There is only one internal equilibrium if it exists (See Fig. 2a-c). For low P level, the internal equilibrium is stable (See Fig. 2a). At intermediate levels of P the internal equilibria are unstable (See Fig. 2b-c). In the case of excessive P level, there is no point of intersection which means that the internal equilibrium does not exist (See Fig. 2d).

14

(a) P = 0.03 mg P/L

(b) P = 0.05 mg P/L

(c) P = 0.08 mg P/L

(d) P = 0.15 mg P/L

Figure 2: Red dashed curves and black solid curve are defined F0 (x, y) and G0 (x, y) with different P values, respectively. A filled circle denotes a stable equilibrium while an open circle denotes unstable equilibrium.

265

266 267 268 269 270 271 272 273 274 275 276 277 278 279 280

4. Numerical Simulations The parameter P represents the total phosphorus into the system. The level of P can change the producer’s quality, and the analytical analyses in the previous section suggest that the change of producers P: C can profoundly affect the food web and predation between different trophic levels. Bifurcation diagrams provide clear and visual ways to understand how behaviors of a dynamical system depend on a parameter. Here we analyze the dynamics of System (3) by choosing P as the bifurcation parameter. The spectrum of the maximum Lyapunov exponent diagram (Fig. 3) and bifurcation diagrams (Fig. 4) provide further insight on how System (3) responds to nutrient enrichment as P increases from 0 to 0.15 mg P/L. We also present time series simulations for four different values of P to assist intuitive analysis (Fig. 5). The parameter values are listed in Table (1). These values are biologically realistic values obtained from Andersen (1997) and Urabe and Sterner (1996) and used by Loladze et al. (2000) and Peace et al. (2013). The value of eˆz is realistic for most planktivorous fish according to Leidy and Jenkins (1977). The parameterization of θz is in the range of juvenile gizzard shad P:C ratios presented in Dickman et al. (2008). The values of gˆ and az are in the ranges of parameter values compiled by MacKenzie et al. (1990) for larval fish ingestion rate. The parameterization of dz 15

281 282 283

includes both natural mortality and respiration rates reported by Leidy and Jenkins (1977). The initial conditions are x(0) = 0.5 mg C/L, y(0) = 0.25 mg C/L, and z(0) = 0.125 mg C/L (the predator population is half of the consumer population that is half of the producer population). Parameter P eˆy eˆz b dy dz θy θz q K f (x) fˆ ay g(y) gˆ az

Description Total phosphorus Maximal production efficiency in C terms for consumer Maximal production efficiency in C terms for predator Maximal growth rate of the producer Consumer loss rate Predator loss rate Consumer constant P:C Predator constant P:C Producer minimal P:C Producer carrying capacity limited by light Consumer ingestion rate Maximal ingestion rate of the consumer Half-saturation of the consumer ingestion response Predator ingestion rate Maximal ingestion rate of the predator Half-saturation of the predator ingestion response

Value 0–0.15 0.8 0.8 1.2 0.25 0.25 0.03 0.04 0.0038 1.5 fˆx/(ay + x) 0.81 0.25 gˆy/(az + y) 0.81 0.25

Units mg P/L – – /day /day /day mg P/mg C mg P/mg C mg P/mg C mg C/L /day /day mg C/L /day /day mg C/L

Table 1: Model paremeters 284 285 286 287 288 289 290 291 292

The usual test of whether a deterministic dynamical system is chaotic or nonchaotic involves the calculation of the maximal Lyapunov exponents (MLE) (Solomon et al. 2003; Wang and Zhao 2016). A positive MLE indicates chaos (Hern´andez and Rodr´ıguez, 2009). In this paper, the MLE is calculated following the procedure in Sprott and Sprott (2003, p. 116-117). The general idea behind this procedure is to follow two nearby orbits and to calculate their average logarithmic rate of separation. Whenever they get too far apart, one of the orbits has to be moved back to the vicinity of the other along the line of separation. Here, we take a conservative approach to this step and do this at each iteration. The estimated MLE values for varying total phosphorus levels are shown in Fig 3.

16

Figure 3: Estimation of the maximal Lyapunov exponent with P as bifurcation parameter. 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308

Fig. 4 displays bifurcation diagrams of System (3) with bifurcation parameter P . Plotted are the minima and maxima of the producer (x), consumer (y), and predator (z) densities from computer simulated time series solutions in the interval [5000, 10000] days. We waited 5,000 days to ensure the transient dynamics had stabilized out. For very low levels of P , the consumer and the predator cannot persist due to starvation and the boundary equilibrium (K, 0, 0) is stable. As P increases, the consumer starts to survive and (¯ x, y¯, 0) is the attractor; P at such level can support the producer and the consumer, but not the predator. As P increases further, all the populations survive and coexist at a stable internal equilibrium (x∗ , y ∗ , z ∗ ) (e.g., P = 0.03 in Fig. 5a). As P is increased through a threshold value, the internal equilibrium loses its stability through a Hopf bifurcation. As a result, a limit cycle is born, and as P increases, the amplitudes of the oscillations increase. When P increases to some intermediate level (e.g., P = 0.05 in Fig. 5b), the system exhibits chaotic behavior in all population densities. The calculation of the MLE effectively supports this claim (see Fig. 3). When P increases through an intermediate threshold value, the dynamics of the system change abruptly, the chaotic behavior is transformed into periodic oscillations (e.g., P = 0.08 in Fig. 5c). For excess values of P , the consumer is heading toward deterministic extinction due to low food quality while the predator goes extinct because of low food quantity (e.g., P = 0.15 in Fig. 5d).

17

(a) Producer vs. P

(b) Consumer vs.P

(c) Predator vs. P Figure 4: Bifurcation diagrams of the System (3) with bifurcation parameter total phosphorus for the (a) producer, (b) consumer, and (c) predator densities. Parameter values are listed in Table 1 where total phosphorus varies from 0 to 0.15 mg P/L. Data were generated using MATLAB.

18

(a) P = 0.03 mg P/L

(b) P = 0.05 mg P/L

(c) P = 0.08 mg P/L

(d) P = 0.15 mg P/L

Figure 5: Numerical simulations of the full model presented in System (3) performed using parameters found in Table 1 and varying values for P , (a) low total phosphorus P = 0.03 mg P/L, (b) P = 0.05 mg P/L, (c) P = 0.08 mg P/L, (d) excess total phosphorus P = 0.15 mg P/L. Producer, consumer, and predator densities (mg C/L) are given by solid black, blue dashed, and red dash-dot lines, respectively.

309

310 311 312 313 314 315 316 317 318 319 320 321

5. Discussion We constructed and analyzed a tritrophic stoichiometric knife-edge model. The model stresses the importance of incorporating the effects of both food quantity and food quality into food chains. Conditions for positivity and boundedness of solutions are obtained in Theorem 3.2. The existence and stability conditions of boundary and internal equilibria are analyzed in Theorem 3.3, Theorem 3.4, and Theorem 3.5. For each species, bifurcation analysis of the model with total P is performed. We observe that the system shows complex dynamics such as chaotic and periodic oscillations under stoichiometric constraints. We calculate the maximum Lyapunov exponent to measure the chaotic behavior of the system. Chaotic dynamics have been observed in classical food chains that neglect stoichiometric constraints. Hastings and Powell (1991) demonstrated that a simple linear food chain with three trophic levels can exhibit chaotic dynamics. Several aspects of our model are similar to that of Hastings and Powell (1991); each considers three trophic levels where the primary producer follows 19

322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353

logistic growth, the consumer predates solely on the primary producer, and the predator predates solely on the consumer, each following Holling type II saturating functional responses. However, our model incorporates additional complexities due to stoichiometric constraints. Hastings and Powell (1991) warned that the existence of chaotic behavior may depend very delicately on the exact form of the model equations. Here, we show that modifying the equations to consider stoichiometric constraints can also yield chaotic behaviors. Additionally, this enables us to explore how the environmentally relevant parameter value of the phosphorus load can lead to chaotic dynamics (Fig. 3). Our model made the simplifying assumptions that the producer is stoichiometrically variable, but the consumer and predator have constant nutrient contents, called “fixed stoichiometry” or “strict homeostasis.” Under these strict homeostatic assumptions for the consumer and predator, our model suggested new insights into how energy flow and nutrient cycling regulates tritrophic interactions. When the producer is nutrient limited, the consumer is limited by nutrient, and can even go extinct despite an abundance of food. As a natural result of this, the predator also goes extinct due to lack of food. However, empirical studies confirm that under natural conditions many consumers are not strictly homeostatic, as consumers consistently display some small changes in body elemental composition in response to stoichiometric variations in their diet (Laspoumaderes et al., 2010). In the past decade, some theoretical studies have developed a new generation of stoichiometric models with nonhomeostatic consumers (Grover 2003; Mulder and Bowden 2007; Mulder 2007; Wang et al. 2012, 2018). But all of them were restricted with the first two trophic levels of the food chain. It is, therefore, worth developing a stoichiometric model with a strictly homeostatic predator that allows variable nutrient contents in both the producer and the consumer. Our model is a stoichiometrically explicit tritrophic (producer, consumer, and natural enemy of consumer) food chain model where interactions between species are observed as either predation (+, −) or competition (−, −). However, interactions in natural communities are far more complex than these, for example, intraguild predation (IGP), defined as killing and eating among potential competitors, seems to be a ubiquitous interaction, differing from competition or predation (Polis et al., 1989). The study of IGP has led to reconsidering many classical topics such as the stability and diversity of communities, trophic cascades in food webs, and the effects of ecosystem productivity on species interactions (Polis et al. 1989; Holt and Polis 1997; Borer et al. 2003; Holt and Huxel 2007.) Therefore, construction of a stoichiometrically explicit model of IGP with the knife-edge phenomenon may provide valuable and useful insights into better understanding population dynamics.

354

Funding

355

This research is partially supported by NSF grant DMS-1615697.

356

Acknowledgements

358

The authors wish to express their sincere gratitude to the anonymous referee for his/her constructive suggestions and comments have improved the manuscript a great deal.

359

References

360

Andersen T (1997) Pelagic Nutrient Cycles: Herbivores as Sources and Sinks. New York: Springer

357

361 362

Asik L, Peace A (2019) Dynamics of a producer–grazer model incorporating the effects of phosphorus loading on grazers growth. Bulletin of mathematical biology 81(5):1352–1368 20

363 364 365

366 367

368 369

370 371 372

373 374

375 376

Borer ET, Briggs CJ, Murdoch WW, Swarbrick SL (2003) Testing intraguild predation theory in a field system: does numerical dominance shift along a gradient of productivity? Ecology Letters 6(10):929–935 Branco P, Egas M, Elser JJ, Huisman J (2018) Eco-evolutionary dynamics of ecological stoichiometry in plankton communities. The American Naturalist 192(1):E1–E20 Chen M, Fan M, Kuang Y (2017) Global dynamics in a stoichiometric food chain model with two limiting nutrients. Mathematical biosciences 289:9–19 Dickman EM, Newell JM, Gonz´ alez MJ, Vanni MJ (2008) Light, nutrients, and food-chain length constrain planktonic energy transfer efficiency across multiple trophic levels. Proceedings of the National Academy of Sciences 105(47):18,408–18,412 Diehl S (2007) Paradoxes of enrichment: effects of increased light versus nutrient supply on pelagic producer-grazer systems. The American Naturalist 169(6):E173–E191 Elser JJ, Loladze I, Peace AL, Kuang Y (2012) Lotka re-loaded: modeling trophic interactions under stoichiometric constraints. Ecological Modelling 245:3–11

378

Grover JP (2003) The impact of variable stoichiometry on predator-prey interactions: a multinutrient approach. The American Naturalist 162(1):29–43

379

Hastings A, Powell T (1991) Chaos in a three-species food chain. Ecology 72(3):896–903

377

380 381

382 383

384 385

386 387

388 389 390

391 392

393 394 395

396 397

398 399

Hern´andez JBA, Rodr´ıguez PH (2009) Nonlinear techniques for signals characterization. In: Encyclopedia of Artificial Intelligence, IGI Global, pp 1266–1272 Holt RD, Huxel GR (2007) Alternative prey and the dynamics of intraguild predation: theoretical perspectives. Ecology 88(11):2706–2712 Holt RD, Polis GA (1997) A theoretical framework for intraguild predation. The American Naturalist 149(4):745–764 Laspoumaderes C, Modenutti B, Balseiro E (2010) Herbivory versus omnivory: linking homeostasis and elemental imbalance in copepod development. Journal of Plankton Research 32(11):1573–1582 Leidy GR, Jenkins RM (1977) The development of fishery compartments and population rate coefficients for use in reservoir ecosystem modeling. Department of Defense, Department of the Army, Corps of Engineers, Waterways Loladze I, Kuang Y, Elser J (2000) Stoichiometry in producer-grazer systems: Linking energy flow with element cycling. Bull Math Biol 62:1137–1162 MacKenzie B, Leggett W, Peters R, et al. (1990) Estimating larval fish ingestion rates: Can laboratory derived values be reliably extrapolated to the wild?. Marine ecology progress series Oldendorf 67(3):209–225 Mulder K (2007) Modeling the dynamics of nutrient limited consumer populations using constant elasticity production functions. ecological modelling 207(2-4):319–326 Mulder K, Bowden WB (2007) Organismal stoichiometry and the adaptive advantage of variable nutrient use and production efficiency in daphnia. Ecological Modelling 202(3-4):427–440 21

400 401

402 403

404 405 406

407 408

409 410

411 412

413 414

Muller EB, Nisbet RM, Kooijman SA, Elser JJ, McCauley E (2001) Stoichiometric food quality and herbivore dynamics. Ecology Letters 4(6):519–529 Peace A (2015) Effects of light, nutrients, and food chain length on trophic efficiencies in simple stoichiometric aquatic food chain models. Ecol Model 312:125–135 Peace A, Zhao Y, Loladze I, Elser JJ, Kuang Y (2013) A stoichiometric producer-grazer model incorporating the effects of excess food-nutrient content on consumer dynamics. Math Biosci 244(2):107–115 Peace A, Wang H, Kuang Y (2014) Dynamics of a producer–grazer model incorporating the effects of excess food nutrient content on grazers growth. Bulletin of mathematical biology 76(9):2175–2197 Plath K, Boersma M (2001) Mineral limitation of zooplankton: stoichiometric constraints and optimal foraging. Ecology 82(5):1260–1269 Polis GA, Myers CA, Holt RD (1989) The ecology and evolution of intraguild predation: potential competitors that eat each other. Annual review of ecology and systematics 20(1):297–330 Rosenzweig ML, MacArthur RH (1963) Graphical representation and stability conditions of predatorprey interactions. The American Naturalist 97(895):209–223

416

Solomon TH, Wallace BR, Miller NS, Spohn CJ (2003) Lagrangian chaos and multiphase processes in vortex flows. Communications in Nonlinear Science and Numerical Simulation 8(3-4):239–252

417

Sprott JC, Sprott JC (2003) Chaos and time-series analysis, vol 69. Citeseer

415

418 419

420 421

422 423

424 425

426 427

428 429

Sterner RW (1990) The ratio of nitrogen to phosphorus resupplied by herbivores: zooplankton and the algal competitive arena. The American Naturalist 136(2):209–229 Sterner RW, Elser JJ (2002) Ecological stoichiometry: the biology of elements from molecules to the biosphere. Princeton University Press Urabe J, Sterner RW (1996) Regulation of herbivore growth by the balance of light and nutrients. Proc Nat Acad Sci 93(16):8465–8469 Wang H, Sterner RW, Elser JJ (2012) On the strict homeostasis assumption in ecological stoichiometry. Ecological Modelling 243:81–88 Wang H, Lu Z, Raghavan A (2018) Weak dynamical threshold for the strict homeostasis assumption in ecological stoichiometry. Ecological modelling 384:233–240 Wang Y, Zhao J (2016) Advances in Energy, Environment and Materials Science: Proceedings of the International Conference on Energy, Environment and Materials Science. CRC Press

22

Credit Author Statement 430

Lale Asik: Conceptualization, Methodology, Software, Formal analysis, Writing - Original Draft, Writing - Review & Editing. Ming Chen: Software. Angela Peace: Methodology, Software, Writing - Review & Editing, Funding acquisition.