Mathematics in Applied Sciences and Engineering https://ojs.lib.uwo.ca/mase Volume 3, Number 1, March 2022, pp.24-49 https://doi.org/10.5206/mase/14332 DYNAMICS OF A STOICHIOMETRIC PRODUCER-GRAZER MODEL WITH MATURATION DELAY HUA ZHANG, HAO WANG, AND BEN NIU Abstract. Ecological stoichiometry provides a multi-scale approach to study macroscopic phenomena via microscopic lens. A stoichiometric producer-grazer model with maturation delay is proposed and studied in this paper. The interaction between stoichiometry and delay is novel and leads to more interesting insights beyond classical delay-driven periodic solutions. For example, the period doubling bifurcation route to chaos can occur as the minimal phosphorous:carbon ratio in producer decreases. Mathematically, we establish the conditions for the existence and stability of positive equilibria, and study the occurrence of Hopf bifurcation at positive equilibria. Analytic results show that delay can change the number and stability of positive equilibria through transcritical bifurcation, saddle-node bifurcation and Hopf bifurcation, and it further determines the grazer’s extinction. Our model with a small delay behaves like LKE model in terms of light intensity, and Rosenzweig’s paradox of enrichment exists in a suitable light intensity. We plot bifurcation diagrams and show rich dynamics driven by delay, light intensity, phosphorous availability, and conversion efficiency, including that a large delay can drive the grazer to go extinct in an intermediate light intensity that is favorable for the survival of the grazer when there is no delay; a limit cycle can appear, then disappear as the delay increases; given the same initial condition, solutions with different delay values can tend to different attractors. 1. Introduction There is increasing evidence that elemental imbalances between producer and grazer can significantly influence their growth, reproduction and survival. For example, the experiment studying zooplankton- phytoplankton interactions [32] showed zooplankton growth may suffer at high algal density instead of always being positively correlated with algal density. This inspires many researchers to explain the producer-grazer interactions from the stoichiometric point of view, which focuses on the relations between multiple key elements in organisms and the abiotic environment [22]. Carbon (C, supplying energy) and phosphorus (P, measuring nutrient) are two vital elements for a cell that is the basic unit of living organisms. One classical stoichiometric producer-grazer model known as LKE model [19] tracks how energy flow and nutrient cycling affect the grazer’s dynamics [17, 37]. The LKE model allows P:C (phosphorous:carbon ratio) in producer to vary above a minimum structural P:C while to keep constant in grazer under the ”strict homeostasis” assumption [36, 35]. If P:C in producer is greater than that of grazer, then producer becomes low-quality food. Thus, this model incorporates the effect of both producer’s quantity (light dependent) and quality (nutrient dependent) on the grazer’s dynamics. Received by the editors 15 October 2021; accepted 2 February 2022; published online 6 February 2022. 2020 Mathematics Subject Classification. 34K18, 37D45, 37N25, 92D40. Key words and phrases. stoichiometric producer-grazer model, maturation delay, stability, bifurcation analysis, light intensity, phosphorus. Hua Zhang and Ben Niu were both supported by Shandong Provincial Natural Science Foundation (ZR2019QA020). Hao Wang was supported by Natural Sciences and Engineering Research Council of Canada (Nos. RGPIN-2020-03911 and RGPAS-2020-00090). 24 https://ojs.lib.uwo.ca/mase https://dx.doi.org/https://doi.org/10.5206/mase/14332 STOICHIOMETRIC PRODUCER-GRAZER MODEL WITH MATURATION DELAY 25 The LKE model reads dx dt = bx ( 1 − x min{K, (P −θy)/q} ) −f(x)y, dy dt = ê min { 1, P −θy θx } f(x)y −dy, (1.1) where x, y are the densities of producer and grazer, respectively. b is the intrinsic growth rate; K is the carrying capacity of producer, which is assumed to be determined positively by light intensity; ê is the maximal conversion efficiency of grazer; d is the loss rate of grazer involving metabolic losses and death; q is the minimal P:C in producer; θ is P:C in grazer keeping a constant; P is the total mass of phosphorus in the entire ecosystem system. All parameters are positive constants, and ê < 1, q < θ. f(x) reflects grazer ingestion rate, and it is a bounded smooth function with f(0) = 0, f ′ (0) < ∞, f ′ (x) > 0, for x ≥ 0. Note that limx→0 f(x) x = f ′ (0) < ∞, which implies that system (1.1) is meaningful as x → 0 although there is the term min { 1, P−θy θx } . Reasonable types of f(x) includes Holling I, Holling II and Holling III functions. In [19], it was revealed that model (1.1) with Holling II type of f(x) admitted multiple positive equilibria, limit cycles and various bifurcation phenomena. The authors also observed the paradox of energy enrichment, where intense energy enrichment substantially elevated the producer density but decreased the grazer growth rate and may drive the grazer to extinction. An explanation on such phenomenon is that the grazer becomes poor-quality food in an extremely high light intensity due to amounts of carbon element is fixed by photosynthesis. Complete analysis for model (1.1) on local and global stability of all equilibria and existence of limit cycles was provided by [17] and [37]. The authors in [17] dealt with two cases of f(x): Holling type I and Holling type II. For the former case, theoretical analysis showed the unique internal equilibrium was always globally asymptotically stable if it existed. For the later case, limit cycles, bistability and several bifurcation types were exhibited when all parameters were fixed at realistic values except K varied. This work was enriched by [37], where authors gave a comprehensive dynamics analysis for Holling II type of f(x) with all flexible parameters. Three new types of bistability comparing with [17] were found: between two positive equilibria, between one positive equilibrium and one boundary equilibria, between the limit cycle and one boundary equilibrium. For model (1.1), Yuan et al. [39] accounted for the effect of environmental noises on the switch between two stochastic attractors in the bistable situation. WKL model in [34] mechanistically incorporates free nutrient in media. Other extensions and applications of LKE model can be found in [33, 24, 38, 16, 3] and the references therein. These studies make contributions to apply stoichiometric models for the improvements in the predictive power of population ecology and cancer treatment. Notably, all existing stoichiometric models implicitly assume an instant process for grazer to be capable of preying on producer for its own growth. Nevertheless, grazer often needs some time to become mature so that it has the ability to prey producer. Inspired by [12], we assume grazer has two stage groups: immature grazer (yj) and mature grazer (y), and only mature grazer lives on producer. Therefore, we propose a stoichiometric producer-grazer model with stage-structure for the grazer as 26 H. ZHANG, H. WANG, AND B. NIU follows: dx(t) dt = bx(t) ( 1 − x(t) min{K, (P −θy(t))/q} ) −f(x(t))y(t), dy(t) dt = ê min { 1, P −θy(t− τ) θx(t− τ) } f(x(t− τ))y(t− τ)e−µτ −dy(t), dyj(t) dt = ê min { 1, P −θy(t) θx(t) } f(x(t))y(t) −µyj(t) − ê min { 1, P −θy(t− τ) θx(t− τ) } f(x(t− τ))y(t− τ)e−µτ, (1.2) where τ is the maturation delay, µ is the mortality rate of immature grazer, and e−µτ is the survival rate of immature grazer. System (1.2) can be derived from the standard age-structured population model, for example, see [28, 2]. Here, we include it for the completeness. Let Y (t,a) be the density of grazer of age a at time t, τ is the maturation period, µ is the mortality rate of immature grazer, and d is the loss rate of mature grazer. Assume that Y (t,a) satisfies the following age-structured population model ∂Y (t,a) ∂t + ∂Y (t,a) ∂a = −µY (t,a), t > 0, 0 < a < τ, (1.3) and the mature grazer density y(t) := ∫∞ τ Y (t,a)da safisfies dy(t) dt = Y (t,τ) −dy(t), t > 0, (1.4) with Y (t, 0) = y(t). Fix s > 0 and let ws(t) := Y (t,t−s) for s ≤ t ≤ s + τ. Together with model (1.3), we have dws(t) dt = −µws(t), 0 ≤ s ≤ t ≤ s + τ, with ws(s) = y(s). This leads to ws(t) = e−µ(t−s)y(s). Therefore, Y (t,τ) = wt−τ (t) = e−µτy(t− τ). Substituting Y (t,τ) into Eq. (1.4), and assuming that the prey activity of the mature grazer follows from that in system (1.1), we can obtain the second equation of system (1.2). It is reasonable to assumed that the density of immature grazer is small compared to that of mature grazer, so we only consider the dynamics of the first two equations in (1.2), that is, dx(t) dt = bx(t) ( 1 − x(t) min{K, (P −θy(t))/q} ) −f(x(t))y(t), dy(t) dt = ê min { 1, P −θy(t− τ) θx(t− τ) } f(x(t− τ))y(t− τ)e−µτ −dy(t). (1.5) The maturation delay in produce-grazer models has been widely studied, and it usually changes the stability/instability of equilibrium, for example, see [21, 4, 25, 18], but the interaction between stoi- chiometry and delay is novel. In this paper, we mainly study how delay affects the existence of positive equilibria and different types of bifurcation phenomena for system (1.5). As mentioned above, light in- tensity K and phosphorus availability P are two primary limiting factors for determining the persistence or extinction of the grazer. We also explore the change of their roles with the introduction of delay. The existence of delay and minimum function makes it challenging to give deeper theoretical analysis on the dynamics of the stoichiometric system, but numerical simulations provide another effective way to understand the dependence of underlying dynamics on delay and some key parameters in the stoi- chiometric system. Some interesting results relative to stoichiometry and delay can be summarized as follows, and the implications of these results are provided in Section 4. STOICHIOMETRIC PRODUCER-GRAZER MODEL WITH MATURATION DELAY 27 (1) There are two types of coexistence: a stable periodic solution and a locally asymptotically stable (LAS) positive equilibrium; two LAS positive equilibira, see Fig 1. (2) Our model with a small delay behaves like LKE model in terms of light intensity, and Rosen- zweig‘s paradox of enrichment can occur for a suitable light intensity and small delay, but a large delay annihilates the oscillations, see Figs 3 and 7. (3) The delay can not only produce a periodic solution (see Appendix), but it can annihilate a periodic solution, see Figs 4 and 5, where the amplitude of the periodic solution increases over delay, and an increasing delay drives the periodic solution to collide with the critical line x + y = p, which can make the solution change its convergent state. (4) There is a period doubling bifurcation route to chaos when both stoichiometry and delay are incorporated into the system, see Figs 9 and 10. The rest of the paper is organized as follows. In Section 2, we present the basic property of solutions to system (1.5) and the stability of boundary equilibria. With Holling II type functional response, the existence and stability of positive equilibria are discussed, and some conditions for the occurrence of Hopf bifurcation are obtained. In Section 3, using biologically meaningful parameter values, we depict some bifurcation diagrams to illustrate the conclusions obtained in Section 2. Further, we exhibit rich dynamics, and comprehensively show how the grazer’s dynamics depends on key parameters. In Section 4, we relate our analytic results to some important biological phenomena. Finally, we conclude the main results and suggest directions for future research. 2. Basic analysis In this section, we first establish the nonnegativity and boundedness of solutions, then we consider the existence and stability of non-trivial equilibria. For simplicity, let p = P θ and s = q θ , and then system (1.5) becomes dx(t) dt = bx(t) ( 1 − x(t) min{K, (p−y(t))/s} ) −f(x(t))y(t), dy(t) dt = ê min { 1, p−y(t− τ) x(t− τ) } f(x(t− τ))y(t− τ)e−µτ −dy(t). (2.1) 2.1. Nonnegativity and boundedness. In the biological perspective, the initial conditions are given as x(η) ≥ 0, q > y(η) ≥ 0,η ∈ [−τ, 0]. (2.2) Denote k = min{K,P/q}, the basic property of solutions with initial values (2.2) is stated in the following theorem. Theorem 2.1. Let (x(t),y(t)) be any solution of system (2.1) subject to initial conditions (2.2). Then, x(t) ≥ 0 for t ∈ (0,∞), and lim sup t→∞ x(t) ≤ k. Moreover, if M := maxx∈[0,k] f(x) < dê , then 0 ≤ y(t) ≤ p for t ∈ (0,∞). Proof. Solving x(t) from the first equation of (2.1) gives x(t) = x(0)e ∫ t 0 [ b ( 1− x(u) min{K,(p−y(u))/s} ) −f(x(u)) x(u) y(u) ] du . Thus, x(t) ≥ 0 for all t > 0, and x(t) > 0 if x(0) > 0. Moreover, we see that dx(t) dt ≤ bx(t) ( 1 − x(t) min{K,p/s} ) = bx(t) ( 1 − x(t) k ) , From the standard comparison argument, we have lim sup t→∞ x(t) ≤ k. 28 H. ZHANG, H. WANG, AND B. NIU Using the variation-of-constant formula to the second equation in (2.1), we have y(t) = y(0)e−dt + ∫ t 0 ê min { 1, p−y(u− τ) x(u− τ) } f(x(u− τ))y(u− τ)e−µτe−d(t−u)du (2.3) which indicates y(t) ≥ y(0)e−dt for t ∈ [0,τ] provided that p > y(η) ≥ 0, η ∈ [−τ, 0]. Moreover, if y(η) 6≡ 0 for η ∈ [−τ, 0], then y(t) > 0, t ∈ [0,τ]. For t ∈ [τ, 2τ], (I) if 0 ≤ y(t) < p, t ∈ [0,τ], then by a similar way, we have y(t) ≥ 0, and y(t) > 0 when y(η) 6≡ 0 for η ∈ [−τ, 0]. (II) If there exists t1 ∈ [0,τ] such that y(t1) tends to p from below, and 0 ≤ y(t) < p for t ∈ [0, t1), we claim that under the condition M < d ê , it holds that 0 < y(t) ≤ p for t ∈ [τ, 2τ]. It follows from the second equation of (2.1) that dy(t1) dt = ê min { 1, p−y(t1 − τ) x(t1 − τ) } f(x(t1 − τ))y(t1 − τ)e−µτ −dy(t1) ≤ êMy(t1) −dy(t1). Thus, dy(t1) dt < 0 when M < d ê , which implies that for any small ε > 0, y(t1 +ε) < y(t1). Thus, y(t) ≤ p, t ∈ [0, t1 + ε]. We also see that y(t) ≤ p for t ∈ [t1,τ] with the similar argument. Using (2.3), we have y(t) > y(0)e−dt ≥ 0 for t ∈ [t1 + ε,t1 + ε + τ]. As a consequence, 0 < y(t) ≤ p for t ∈ [0,τ]. Similar to (I), we prove the claim. We can repeat the process for any [nτ, (n + 1)τ], n ≥ 2. � Remark 2.1. Due to the existence of a minimum function in the second equation of system (2.1), it is difficult to determine the global existence and the nonnegativity of the component y(t). When there is a strong restriction on f(x), Theorem 2.1 gives the result. Actually, when f(x) does not satisfy the restriction given in Theorem 2.1, 0 ≤ y(t) < p can also hold (see section 3), but we can not prove it theoretically. Such a minimum function can induce complex dynamics that is difficult to provide rigorous proofs, see subsection 3.3. 2.2. Stability of boundary equilibria. System (2.1) always has equilibria: E0 = (0, 0) and E1 = (k, 0), and their local stability is obtained by the distribution of eigenvalues corresponding to the linearized system. Theorem 2.2. (i) E0 is always unstable. (ii) If k < p, E1 is LAS (locally asymptotically stable) when êf(k) < d for all τ ≥ 0. (iii) If k > p, E1 is LAS when f(k) k êp < d for all τ ≥ 0. Proof. The characteristic equation at E0 is (λ− b)(λ + d) = 0, (2.4) this suggests that two eigenvalues have different signs, so E0 is always unstable. At E1, if k < p, then the characteristic equation takes the following form (λ + b)[λ + d− êf(k)e−µτe−λτ ] = 0. (2.5) Note that the stability of E1 depends on the real parts of the zeros of λ + d− êf(k)e−µτe−λτ . When τ = 0, λ = êf(k) −d. Clearly, E1 is LAS if êf(k) < d and is unstable as êf(k) > d. When τ > 0, it is easy to see that λ = 0 is not a root of λ + d− êf(k)e−µτe−λτ = 0 provided êf(k)e−µτ 6= d. Moreover, λ = ±iω (ω > 0) are a pair of purely imaginary roots of (2.5) if and only if ω satisfies cos ωτ = d êf(k)e−µτ , sin ωτ = −ω êf(k)e−µτ , and ω2 = [êf(k)e−µτ ]2 −d2. STOICHIOMETRIC PRODUCER-GRAZER MODEL WITH MATURATION DELAY 29 Hence, (2.5) has no roots on the imaginary axis if êf(k)e−µτ < d. Then we claim that all roots of (2.5) have negative real parts if êf(k) < d. For k > p, the characteristic equation is (λ + b) [ λ + d− f(k) k êpe−µτe−λτ ] = 0. (2.6) Similarly, we obtain E1 is LAS when f(k) k êp < d. The proof is completed. � 2.3. Existence and stability of positive equilibria. In this part, the existence of positive equilibria and their local stability are investigated with f(x) = cx a+x . Motivated by [37], in the rest paper, we always assume (A) K ≤ { P q , θad q(cêe−µτ −d) } , ad + pd cp < êe−µτ. As a consequence, the positive equilibria are determined by{ 0 = b(1 − x K ) − cy a+x , 0 = cêx a+x e−µτ −d, or { 0 = b(1 − x K ) − cy a+x , 0 = cê(p−y) a+x e−µτ −d. Denote h(x) = b c ( 1 − x K ) (a + x), l1 : x ∗ = ad cêe−µτ −d , l2 : y = ( p− da êc eµτ ) − d êc eµτx, x ∈ [ x∗, min { K, cêp d e−µτ −a }] . (2.7) When the parabola y = h(x) and the line l1 intersect in the region {(x,y) ∈ R2 : 0 < x < K, 0 < y,x + y < p}, there exists a positive equilibrium, denoted by E2. As the parabola y = h(x) enters the region {(x,y) ∈ R2 : 0 < x < K, 0 < y,x + y > p}, when it is tangent with the line l2, a newly positive equilibrium exists; when it intersects with the line l2 at two different points, system (2.1) has two newly positive equilibria, denoted by E3 and E4. Thus, system (2.1) may have zero, one, two or three positive equilibria. To explore the existence of positive equilibria, it is convenient to define some critical curves, following the following steps: (i) Let K1(τ) = ad/(cêe −µτ −d). (ii) If the line l1, the parabola y = h(x) and the critical line y = p−x intersect, then b c ( 1 − x∗ K ) (a + x∗) = p−x∗. Denote K2(τ) = x ∗/[1 − c(p−x∗)/(b(a + x∗)], where x∗ is given in (2.7). (iii) If the parabola y = h(x) is tangent to the line l2 at (x̃, ỹ), then we have  ỹ = b c ( 1 − x̃ K ) (a + x̃), b cK (K −a− 2x̃) = − d êc eµτ, ỹ = ( p− da êc eµτ ) − d êc eµτx̃, x̃ ∈ [ x∗, min { K, cêp d e−µτ −a }] . Vanishing x̃ and ỹ yields( b + deµτ ê )2 K2 + [ 2ab ( b + deµτ ê ) − 4bcp ] K + (ab)2 = 0. (2.8) Therefore, if (2.8) has positive roots K+3 (τ) or K − 3 (τ) (it can hold that K + 3 (τ) = K − 3 (τ)), then the parabola y = h(x) is tangent to the line l2 in the region {(x,y) ∈ R2 : x ≥ x∗,x + y ≥ p} when K = K+3 (τ). 30 H. ZHANG, H. WANG, AND B. NIU (iv) Denoting the intersection point of the line l2 and x−coordinate as K4, we have K4(τ) = cêp d e−µτ − a. One can check that if τ satisfies assumption (A), then 0 < K1(τ) < p < K4(τ) for each τ. If cp < ab or cp > ab, êe−µτ < ad+pd cp−ab is further true, then 0 < K1(τ) < K2(τ) for each τ. K±3 (τ) exist if and only if cp ≥ a ( b + d êe−µτ ) , and K±3 (τ) = 2cpb−ab ( b + d êe−µτ ) ± 2 √ b2cp [ cp−a ( b + d êe−µτ )] ( b + d êe−µτ )2 . (2.9) Note that if K2(τ) > 0 and K + 3 (τ) exists, then it must hold K + 3 (τ) < K2(τ). Moreover, K+3 (τ) < K4(τ). Based on [37, Theorems 3.1-3.7], we state the existence of positive equilibria in the following theorem. Theorem 2.3. Assume that (A) is satisfied. (i) Suppose that K ∈ (0,K1(τ)) or K ∈ (max{K2(τ),K4(τ)},∞). Then system (2.1) has no positive equilibria. (ii) Suppose that cp < ab or cp > ab but êe−µτ < ad+pd cp−ab holds. If K ∈ (K1(τ),K2(τ)), then system (2.1) has the positive equilibrium E2. (iii) Suppose that cp > ab and ad cp−ab < êe −µτ hold. If K ∈ (K+3 (τ),K2(τ)), then system (2.1) has the positive equilibrium E3; If K ∈ (K+3 (τ),K4(τ)), then system (2.1) has the positive equilibrium E4. (iv) Suppose that cp > ab, and ad cp−ab < êe −µτ < ad+pd cp−ab hold. If K belongs to sets (K1(τ),K2(τ)) ∩ (K+3 (τ),K2(τ)), (K1(τ),K2(τ)) ∩ (K + 3 (τ),K4(τ)), and (K + 3 (τ),K2(τ)) ∩ (K + 3 (τ),K4(τ)), respectively, then system (2.1) has two positive equilibria: E2 and E3, E2 and E4, or E3 and E4, respectively. If (K1(τ),K2(τ)) ∩ (K+3 (τ),K2(τ)) ∩ (K + 3 (τ),K4(τ)) 6= ∅, then in this interval, system (2.1) has three positive equilibria: E2, E3 and E4. Remark 2.2. For the critical cases, we have (i) when K = K1(τ), that is, τ = 1 µ ln cêK d(a+K) , E2 collides E1 and disappears, system (2.1) undergoes a transcritical bifurcation. Recalling Theorem 2.2, we see that the boundary equilibrium E1 is LAS as E2 disappears, i.e., τ > 1 µ ln cêK d(a+K) . (ii) When K = K2(τ), E2 and E3 merge into one positive equilibrium. System (2.1) undergoes a saddle-node bifurcation, see Fig 2(a). (iii) When K = K+3 (τ), E3 and E4 collide into one positive equilibrium. System (2.1) undergoes a saddle-node bifurcation, see Fig 2(b). (iv) When K = K4(τ), namely, τ = 1 µ ln cêp d(a+K) , E4 collides E1 leaving one boundary equilibrium. System (2.1) undergoes a transcritical bifurcation. Moreover, Theorem 2.2 shows that E1 is LAS if τ > 1 µ ln cêp d(a+K) . The local stability of positives equilibria that are not on the critical line y = p−x can be analyzed using the method in [5], see Appendix for details. We summarize the existence and stability of positive equilibria of system (2.1) in Table 1, where I1, I2 and Sn(τ) are defined in Appendix. 3. Numerical simulations In this section, numerical simulations are carried out to illustrate theoretical results and reveal some biological mechanisms intuitively. Referring to [19], we take P = 0.025, ê = 0.6, b = 1.2, d = 0.25, θ = 0.03, q = 0.0038, c = 0.81, a = 0.25, µ = 0.003. (3.1) STOICHIOMETRIC PRODUCER-GRAZER MODEL WITH MATURATION DELAY 31 Table 1. The existence and stability of positive equilibria for system (2.1) c êe−µτ K Stability c < ab p K ∈ (K1(τ),K2(τ)) E2 Hopf bifurcation occurs if Sn(τ), n ∈ N0 has zeros in I1 c > ab p ad+pd cp < êe−µτ < ad cp−ab K ∈ (K1(τ),K2(τ)) E2 Hopf bifurcation occurs if Sn(τ), n ∈ N0 has zeros in I1 ad cp−ab < êe −µτ < ad+pd cp−ab K ∈ (K1(τ),K2(τ)) E2 Hopf bifurcation occurs if Sn(τ), n ∈ N0 has zeros in I1 K ∈ (K+3 (τ),K2(τ)) E3 unstable K ∈ (K+3 (τ),K4(τ)) E4 Hopf bifurcation occurs if Sn(τ), n ∈ N0 has zeros in I2 ad+pd cp−ab < êe −µτ K ∈ (K+3 (τ),K4(τ)) E4 Hopf bifurcation occurs if Sn(τ), n ∈ N0 has zeros in I2 3.1. Biological mechanisms related to τ. It can be seen that c > ab p . Due to assumption (A), we restrict 0 ≤ τ < 134.1278 and K ≤ min{6.5789, 2.0908}. 3.1.1. Joint effect of τ and K on asymptotic dynamics. We present the existence of positive equilibria on τ − K plane, see Fig 1(a). Bifurcation diagrams provide direct understanding about Theorem 2.3 and Remark 2.2, which reflect the species persistence/extinction regulated by τ and K. The pair (τ,K) satisfying K > max{K4(τ),K2(τ)} is in region D1; (τ,K) satisfying K < K1(τ) is in region D7. Theorem 2.3 indicates that system (2.1) has no positive equilibria in the two regions. The pair (τ,K) satisfying K1(τ) < K < K2(τ) is in regions D3, D4, D5 and D6, thus E2 exists in these regions. We further know that in regions D3 and D4, E2 is unstable induced by Hopf bifurcation, and it remains LAS in regions D5 and D6. When K+3 (τ) < K < K2(τ), (τ,K) locates in regions D3 and D5, and E3 exists in these regions. When K+3 (τ) < K < K4(τ), (τ,K) belongs to regions D2, D3 and D5, E4 exists in these regions. All bifurcation curves described in Remark 2.2 are also exhibited in Fig 1(a). Observe that three curves: K = K1(τ), K = K2(τ) and K = K4(τ) intersect at P3, at which two positive equilibria with one merged by E2 and E3 and the other E4 collide with E1 from its two sides, then all positive equilibria disappear and only one boundary equilibrium stays. In fact, the position of P3 on τ − K plane is (τ̂,K1(τ̂)), where τ̂ = 1 µ ln cêp d(a+p) . It is also seen that the curve K = K+3 (τ) is tangent with the curve K = K2(τ) at P2, thus E2, E3 and E4 collide into one positive equilibrium. Meanwhile, it can be calculated that there is a zero eigenvalue for the corresponding characteristic equation. Therefore, we assert P2 is a cusp bifurcation point. The Hopf bifurcation curve near E2 intersects K-axis at (0, 0.7797) and K = K2(τ) at P1, respectively. At point P4, one can see that E3 and E4 collide into one equilibrium, E2 goes to a boundary equilibrium and disappears. As a result, system (2.1) has exactly one positive equilibrium. Biologically, K = K1(τ) (the red curve) and K = K4(τ) (the magenta curve) are two critical curves determining the grazer’s persistence, which holds in the region between the two curves. We see that the survival region gradually declines as τ increases, which follows from the facts that K1(τ) is an increasing function and K4(τ) is a decreasing function of τ, respectively. Therefore, beneficial light intensity for the growth of grazer negatively correlates with τ, and less light is needed for the grazer’s persistence as τ increases. It has been recognized that the producer is low quality food at high light level due to the Conservation Law of Matter. The similar phenomenon is observed in Fig 1(a) that system (2.1) has no positive equilibria for large K. Of course, the extremely low light intensity leads to a very limited 32 H. ZHANG, H. WANG, AND B. NIU τ K 0 50 100 150 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 Assumption (A) is not satisfied E 2 E 4 P 1 E 2 D7 D2 D4 D1 P 2 P 3 D6 D3 E 2 ,E 3 ,E 4 D5 P 4 E 2 ,E 3 ,E 4 Transcritical bifurcation curve No positive equilibrium Saddle−node bifurcation curve Hopf bifurcation curve Transcritical bifurcation curve No positive equilibrium 0 50 100 150 200 0 0.5 1 1.5 2 τ K Assumption (A) is not satisfied D1 D5 D6 D7 P 3 P 4 D3 P 2 D2 P 1 D4 P 0 No positive EquilibriumTranscritical bifurcation curve Saddle−node bifurcation curve Hopf bifurcation curve Transcritical bifurcation curve No positive Equilibrium (a) (b) Figure 1. The bifurcation diagrams on τ − K plane. The red curve is determined by K = K1(τ); the blue curve stands for K = K2(τ); the magenta curve is given by K = K4(τ); the green one is K = K + 3 (τ); the black one is Hopf bifurcation curve at E2. (a) Parameter values are given in (3.1). (b) ê = 0.655, a = 0.3 and other parameter values are the same as those in (3.1). quantity of the producer which drives the grazer to go extinct as well. Moreover, the existence of P3 reflects that the grazer can go extinct at an intermediate light intensity when τ is too large. The above discussions show that system (2.1) can have sustainable oscillations without delay. Does the appearance of limit cycles is substantially impacted by τ? The simulation in Fig 1(b) illustrates it does, where ê = 0.655 and a = 0.3, and other parameters in (3.1) remain fixed. Dynamics exhibited here are almost identical as in Fig 1(a), except that Hopf bifurcation curve intersects the curve K = K2(τ) at two distinct points, denoted by P0 and P1. In this case, both producer and grazer densities keep stable for τ = 0, and they change periodically only when τ is greater than a certain value. Thus, we assert that τ plays a significant role in the oscillatory behavior of solutions. See Appendix for an example on how τ affects the stability of the positive equilibrium E2. In view of Fig 1(a) and (b), system (2.1) undergoes a saddle-node bifurcation when E2 and E3 (or E3 and E4) collide into one positive equilibrium. This phenomenon is presented in Fig 2. 3.1.2. The bifurcation diagrams of the grazer over τ. To reveal how the grazer’s dynamics depend on τ under different light intensities, we choose four levels of solar energy: K = 0.3, K = 0.7, K = 0.87, K = 1.3, and sketch the change of existence and stability of non-trivial equilibria in Fig 3. At low light (K = 0.3), the grazer declines at a stable equilibrium until it dies out. When light intensity is intermediate (K = 0.7), the variation of the grazer density is complex. The grazer density slowly reduces until τ = 1.278 where a supercritical Hopf bifurcation occurs at E2 such that E2 loses its stability, and the grazer density changes periodically for 1.278 < τ < 21.308. At τ = 21.308, a periodic solution disappears and E2 regains stability, the grazer density continues to decline. Two new equilibria merge through a saddle-node bifurcation at τ = 75.6630, then system has two stable equilibria simultaneously, both of which decrease as τ increases. As τ increases further, the grazer eventually becomes extinct. At high light (K = 0.87), whether the grazer density exhibits sustainable oscillations or declines at a stable equilibrium depends on the initial point. At τ = 0.0527, periodic solutions disappear through an infinite period bifurcation, after that the grazer density keeps at a stable equilibrium whatever its initial density is, and then gradually decreases when τ increases. When τ increases to τ = 123.0325, STOICHIOMETRIC PRODUCER-GRAZER MODEL WITH MATURATION DELAY 33 0.76 0.78 0.8 0.82 0.84 0.86 0.88 0.9 0.78 0.8 0.82 0.84 0.86 0.88 0.9 K x + y E 2 E 3 Bifurcation point 0.8 1 1.2 1.4 1.6 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 K x + y E 4 E 3 Bifurcation point (a) (b) Figure 2. The saddle-node bifurcation portraits. Solid and dot curves represent the existence and non-existence of equilibria, respectively. (a) The saddle-node bifurcation related to E2 and E3 at τ = 38.4432. (b) The saddle-node bifurcation related to E3 and E4 at τ = 10. the grazer becomes extinct. For extremely high light (K = 1.3), the trend of the grazer density is similar to that in the case of low light, except that the persistence interval of τ becomes small. Of course, we have known that the extinction mechanisms are opposite in low light and in very high light. The above observations show that both moderate and high light enrichment can produce sustainable producer-grazer oscillations. 3.1.3. The solution behavior over different delays and initial states. Note that unstable E2 and LAS E4 can coexist under certain conditions, such as in regions D3 and D4 in Fig 1. In this case, system (2.1) presents a bistable phenomenon: a stable periodic solution and a stable positive equilibrium, thus the eventual steady state depends on the initial point. Let K = 0.84 and other parameters be the same as those in (3.1). For initial points satisfying x + y < p, solutions converge to a periodic solution at τ = 0; as τ increases, the solution finally goes to E4, shown in Fig 4(a) and (b). Moreover, whatever the initial values are constant, periodic or monotone functions, when they satisfy x + y < p, solutions may converge to E4 if the periodic solutions collide with the critical line x + y = p. With initial values satisfying x + y > p, solutions directly converge to E4 for both τ = 0 and τ > 0 as shown in Fig 4(c) and (d). This tendency is kept when the initial values are constant, periodic or monotone functions satisfying x + y > p. It is pointed out by [17, Theorem 17] that under the following parameter values: P = 0.025, ê = 0.8, b = 1.2, d = 0.25, θ = 0.04, q = 0.004, c = 0.8, a = 0.25, µ = 0.003, (3.2) K = 0.585185 is an important threshold value, at which E2 is unstable and system (2.1) without delay has at least one limit cycle around E2, and E3 and E4 merge into a saddle-node equilibrium E3,4 on the critical line x + y = p. Besides, solutions always tend to a periodic solution. We are particularly interested in the asymptotic behavior of solutions with different initial values if τ is incorporated into the system. Simulation results imply that as τ increases, E3 and E4 leave the critical line x + y = p and separate, E3 is unstable and E4 gains stability. Both solutions initiating from the regions x+y < p and x + y > p converge to the periodic solution when τ = 0, but as τ increases, all solutions converge 34 H. ZHANG, H. WANG, AND B. NIU 0 20 40 60 80 100 120 0 0.02 0.04 0.06 0.08 0.1 τ y Transcritical bifurcation 0 20 40 60 80 100 120 0 0.1 0.2 0.3 0.4 0.5 τ y S ad d le − n o d e b if u rc at io n P er io d ic s o lu ti o n 0 5 0 0.2 0.4 τ y H o p f b if u rc at io n Transcritical bifurcation (a) (b) 0 20 40 60 80 100 120 0 0.1 0.2 0.3 0.4 0.5 τ y 0 0.05 0.45 0.55 τ y P er io d ic so lu ti o n 16 18 20 0.5337 0.534 τ y S ad d le − n o d e b if u rc at io n Transcritical bifurcation 0 20 40 60 80 100 120 −2 0 2 4 6 8 10 12 14 16 x 10 −3 τ y Transcritical bifurcation (c) (d) Figure 3. The bifurcation diagrams of the grazer density over τ. All other parameter values are given in (3.1), except that K = 0.3, K = 0.7, K = 0.87 and K = 1.3 that correspond to (a)-(d), respectively. Solid and dash curves stand for stable and unstable equilibria, respectively. Black line represents a boundary equilibrium, magenta curve is E2, green and brown curves are E3 and E4, respectively. Blue and red curves are the maxima and minima of amplitudes of periodic solutions, respectively. to E4, see Fig 5. When τ = 0, the periodic solution attracts solutions starting from anywhere of the phase plane; introducing τ into the model can increase the maxima of periodic solutions, and once the periodic solution arrives at the critical line x + y = p from the region x + y < p, the solution may converge to E4. 3.1.4. Joint effect of τ and P on asymptotic dynamics. By decreasing the ratio of P:C in the producer such that the chemical composition required and captured is unbalanced for the grazer, light enrichment harms the growth of the grazer and even leads to the extinction. Naturally, we are curious whether we can change phosphorous to balance the adverse effect of light and significantly facilitate the grazer’s persistence. Since p = P θ , where θ reflects the fixed P:C ratio in the grazer, the change of phosphorus availability can be described by varying p. We exhibit the dynamics of system (2.1) on τ −p plane in Fig 6. STOICHIOMETRIC PRODUCER-GRAZER MODEL WITH MATURATION DELAY 35 0 100 200 300 400 500 600 0.4 0.45 0.5 0.55 0.6 t y (t ) τ=1 τ=0 0.2 0.3 0.4 0.5 0.6 0.3 0.35 0.4 0.45 0.5 0.55 x y E 4 (a) (b) 0 50 100 150 200 250 300 0.4 0.5 0.6 t y (t ) τ=1 τ=0 0.2 0.3 0.4 0.5 0.6 0.3 0.4 0.5 0.6 x y E 4 (c) (d) Figure 4. Initial values satisfy x + y < p in (a) and (b). (a) The solution starting from (0.2715, 0.5224) converges to a stable periodic solution when τ = 0 and tends to E4 when τ = 1. (b) When τ = 1, under distinct initial values (0.2665, 0.5224), (0.2665 − 0.2 sin(2t), 0.5224 − 0.2 sin(2t)), (0.2665 + 0.3t, 0.5224 + 0.3t) and (0.2665 − 0.5t, 0.5224 − t), all solutions converge to E4. The blue, magenta, green and black curves correspond to the above initial values, respectively. Initial values are true for x + y > p in (c) and (d). (c) Both solutions initiating from (0.3648, 0.6522) tend to E4 for τ = 0 and τ = 1. (d) When τ = 1 and initial values are taken as (0.3648, 0.6522), (0.3648 + 0.3 sin(2t), 0.6522 + 0.5 sin(2t)), (0.3648 + 0.3t, 0.6522 + 0.3t) and (0.3648 − 0.5t, 0.6522 − t), respectively, E4 keeps an attractor. The blue, magenta, green and black curves are solutions with the above initial values, respectively. It can be observed that system (2.1) admits the unique positive equilibrium E2 as τ and p vary simultaneously. Moreover, when parameter values are in the region between two black vertical lines, E2 is unstable driven by Hopf bifurcation. The increase of phosphorous availability in an ecosystem can increase the chance of grazer’s survival. We further claim that increasing phosphorous can weaken the negative effect of poor-food quality caused by increasing light intensity. This is based on the fact that the survival region remarkably reduces as K increases to extreme as shown in Fig 1, however, 36 H. ZHANG, H. WANG, AND B. NIU 0 100 200 300 400 500 600 0.4 0.42 0.44 0.46 0.48 t y (t ) τ=3 τ=0 0 100 200 300 400 0.35 0.4 0.45 0.5 0.55 t y (t ) τ=3 τ=0 (a) (b) Figure 5. (a) The initial point satisfies x + y < p with (0.1636, 0.4469), the solution converges to a stable periodic solution for τ = 0 and tends to E4 for τ = 3. (b) The initial point satisfies x + y > p with (0.2626, 0.5469). For τ = 0, the solution converges to a stable periodic solution; for τ = 3, the solution eventually tends to E4. comparing (b) and (c) in Fig 6, we observe this tendency stops when p rises. Phosphorous availability has little effect on the oscillatory behavior which mainly depends on delay. 3.2. The bifurcation diagrams over K or ê. In this part, we study the relationship between dy- namical behavior and two other parameters in system (2.1): light-dependent carrying capacity K and conversion efficiency ê. 3.2.1. The bifurcation diagrams of the grazer over K. The work of [19] and the above discussions have shown that the light-dependent carrying capacity K is a key factor on the grazer’s fate. We sketch the bifurcation diagrams of the grazer with respect to K in Fig 7, where τ = 5 in (a) and τ = 100 in (b). In Fig 7 (a), when 0 < K < 0.2732 and K > 1.3459, system (2.1) has no positive equilibria, which suggests the grazer cannot survive in such two scenarios. The extinction of grazer for 0 < K < 0.2732 is caused by the lack of producer, while the grazer can not persist for K > 1.3459 because of a mismatch between the nutrient content in the producer and the nutrient demand of the grazer. The Rosenzweig’s paradox of enrichment [27] emerges in system (2.1) as K varies in the interval (0.2732, 0.8186): the grazer increases steadily for 0.2732 < K < 0.6543. At K = 0.6543, E2 loses its stability and a family of periodic solutions bifurcate from it, then the amplitude of the periodic solution gradually increases with K. As K increases to K = 0.8186, the periodic solution bursts into a heteroclinic orbit, and two new equilibria appear: unstable E3 and stable E4. As K further increases, the grazer density keeping on E4 declines. In other words, higher light intensity leads to lower grazer biomass, which is consistent with the results in [19]. Finally, the grazer becomes extinct at K = 1.3459 caused by a transcritical bifurcation. It can also be seen that E2 and E3 disappear at K = 0.985 through a saddle-node bifurcation. According to Fig 7 (b), we observe that the grazer survives only for 0.5680 < K < 0.9501. As K varies from 0.5680 to 0.6586, the grazer density rises at stable E2, then two new equilibria, E3 and E4, appear that coexist with E2 as 0.6586 < K < 0.7272, where the grazer density will increase if its initial value is small, but it will decrease if its initial value is large. At K = 0.7272, E2 collides with E3 and then disappears. The grazer continues to persist until K = 0.9501. STOICHIOMETRIC PRODUCER-GRAZER MODEL WITH MATURATION DELAY 37 τ p 0 50 100 150 0 0.2 0.4 0.6 0.8 1 1.2 Coexistence region (at E 2 ) Periodic solution exists τ p 0 50 100 150 200 0 0.5 1 1.5 2 2.5 Periodic solution exists Coexistence region (at E 2 ) (a) (b) τ p 0 50 100 150 200 250 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Periodic solution exists Coexistence region (at E 2 ) (c) Figure 6. The bifurcation diagrams on τ − p plane. The red curve is determined by τ = 1 µ ln cêK d(a+K) ; the blue curve stands for p(τ) = b c [1 − x ∗(τ) K ][a + x∗(τ)] + x∗(τ) with x∗(τ) = ad cêe−µτ−d; the green one is p(τ) = 1 4bcK [ (b + de µτ ê )2K2 + 2ab(b + de µτ ê ) + (ab)2 ] ; the magenta curve represents p(τ) = d(K+a) cê eµτ . The black vertical line is Hopf bifurcation line. There are periodic solutions bifurcating from E2 when parameter values are in the region between two vertical lines. Parameter values are given in (3.1) except p varies and (a) K = 0.7, (b) K = 2, (c) K = 4. Comparing with Fig 4 in [19], we find that in the absence or presence of a small delay in system (2.1), the grazer density changes following a similar route, but the ranges of light intensity supporting the grazer’s persistence are different. The persistence window of light intensity seems narrowed by the maturation delay. From the mathematical point of view, delay indeed can induce richer dynamics, for example, the saddle-node bifurcation occurring at E2 and E3. According to Fig 7, we observe that the increase of delay will reduce the survival chance of the grazer, and sustainable oscillations disappear for a very large delay. 3.2.2. The bifurcation diagrams over ê and K. Mathematical models often assume the food assimilation of the grazer is a constant. In reality, the conversion efficiency depends on a variety of factors such as cell morphology and colony architecture [7]. Some studies have pointed out that the paradox of enrichment 38 H. ZHANG, H. WANG, AND B. NIU 0 0.4 0.8 1.2 1.6 2 0 0.2 0.4 K y H o p f b if u rc at io n P er io d ic s o lu ti o n S ad d le − n o d e b if u rc at io n Saddle−node bifurcation Transcritical bifurcation Transcritical bifurcation 0 0.5 1 1.5 0 0.1 0.2 0.3 0.4 0.5 K y S ad d le − n o d e b if u rc at io n S ad d le − n o d e b if u rc at io n Transcritical bifurcation Transcritical bifurcation (a) (b) Figure 7. The bifurcation diagrams of the grazer with respect to K. All other pa- rameter values are given in (3.1) and τ is fixed at τ = 5 for (a) and τ = 100 for (b). All curves have the same interpretations as those in Fig 3. shares such an assumption, see [1, 6, 11] and references therein. How does the conversion efficiency affect the dynamics of system (2.1)? Motivated by this question, we draw a bifurcation diagram in Fig 8(a) to illustrate the existence and stability of positive equilibria on ê−K plane. The coexistence region of producer and grazer is bounded by red and magenta curves. The coexistence window of light intensity becomes wider as ê increases. However, the grazer’s extinction still occurs when the light intensity is too high or too low. 0.45 0.5 0.55 0.6 0.65 0.7 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 K ê No positive equilibrium No positive equilibrium E 4 E 2 , E 3 ,E 4 E 2 , E 3 ,E 4 Transcritical bifurcation curve Saddle−node bifurcation curve Transcritical bifurcation curve E 2 (Stable) Hopf bifur− cation curve E2 (Unstable) 0.4 0.45 0.5 0.55 0.6 0.65 0.7 −0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 y H o p f b if u rc at io n P er io d ic s o lu ti o n S ad d le − n o d e b if u rc at io n ê Transcritical bifurcation (a) (b) Figure 8. (a) The two-dimensional bifurcation diagram on ê−K plane. Take τ = 5 and other parameter values given in (3.1), except that ê varies. The red curve is determined by K(ê) = x∗(ê); the blue curve stands for K(ê) = x∗(ê)/[1 − c(p − x∗(ê))/b(a + x∗(ê))]; the magenta curve is given by K(ê) = cêp d e−µτ − a; the green one is K = K+3 (ê); the black one is Hopf bifurcation curve at E2. Here, x ∗(ê) = ad/(cêe−µτ − d) and K+3 (ê) is defined in (2.9). (b) The bifurcation diagram of the grazer with respect to ê. The curves have the same meanings as those in Fig 3. STOICHIOMETRIC PRODUCER-GRAZER MODEL WITH MATURATION DELAY 39 At medium light, increasing ê can generate or remove positive equilibria, destabilize the system. For example, a bifurcation diagram describing the grazer density over ê is depicted in Fig 8(b) with K = 0.7. Here, system (2.1) has a stable boundary equilibrium for ê < 0.4252. As ê varies from 0.4252 to 0.4861, system (2.1) has three positive equilibria: E2, E3 and E4, where two stable equilibria coexist. Bistability implies that the grazer density converges to which one stable equilibrium depending on its initial value. At ê = 0.4861, unstable E3 and stable E4 disappear simultaneously through a saddle-node bifurcation. As ê further increases, the grazer density rises steadily until ê = 0.5726, where stable E2 loses its stability to a periodic solution. The amplitude of the periodic solution increases as ê increases. In conclusion, the grazer goes extinct if the conversion efficiency is less than 0.4252 and persists for 0.4252 < ê < 0.7. However, we want to claim that the persistence of the grazer can be threatened when the limit cycle has too large amplitude because the low point will be so small that a tiny perturbation/stochasticity will drive it to extinction. Consistent with intuition, increasing the conversion efficiency of a grazer is generally beneficial for its survival, and our study further indicates subtle dependence of the grazer’s dynamics on the conversion efficiency. 3.3. Period doubling route to chaos. We have presented some theoretical and numerical results under assumption (A) in the above discussions. Actually, system (2.1) with f(x) = cx a+x can exhibit more complicate dynamics when ignoring the assumption, for example, the period doubling bifurcation route to chaos. We choose the following values of parameters to describe such phenomenon. P = 0.025, ê = 0.8, b = 1.2, d = 0.25, θ = 0.048, K = 2, µ = 0.004, a = 0.25, c = 0.81. (3.3) Fixing τ = 10, we find that when 0.0125 < q < 0.1490, there is at least one positive equilibrium E5 = (x5,y5) with x5 = 0.1677, and y5 = 1 2 [ b c (a + x5) + p− √ M ] , M = [ b c (a + x5) + p ]2 − 4b c ( p− qx5 θ ) (a + x5). Moreover, system (2.1) takes place the chaos routed by period doubling near E5 as q decreases in the interval [0.0238, 0.034], see Fig 9(a). Here we choose the component x of E5 as a Poincáre section, and draw the change of the solution on the Poincáre section over q. When q > 0.032, E5 is LAS. When 0.0288 < q < 0.032, E5 loses its stability to a stable periodic solution with period 1. When 0.0263 < q < 0.0288, the period-1 solution loses its stability to a stable periodic solution with period 2. As 0.0238 < q < 0.0263, the period-2 solution loses its stability to a stable periodic solution with period- 4, etc., and finally, E5 becomes chaotic. At q = 0.0238, the chaotic attractor suddenly disappears, which may be caused by the collision between the attractor and a periodic solution on the basin boundary of the attractor [23]. As 0.0207 < q < 0.0238, the solution converges to a stable period-1 solution. As q < 0.0207, the solution converges to a boundary equilibrium. The periodic solutions with period 1, 2 and 4 are shown in Fig 10, respectively. When the solution always satisfies x + y < p, there is no period doubling route to chaos near E5, see Fig 9(b). Here, the maxima and minima of the solutions are drawn. It can be seen that E5 is LAS for q > 0.032, then it loses its stability to a period-1 solution that is similar to Fig 10(a) or (b), we do not show it here. For system (2.1) without delay, the period doubling route to chaos is not observed in this paper. As far as we know, chaotic dynamics from LKE model without or with delay is scarce, while it has been frequently displayed in communities consisting of two or more species, see [20, 31, 29, 30] For the LKE model incorporating the maturation time of the grazer and the restriction of carbon and phosphorus elements in the producer, we see that the decrease of the minimal P:C in the producer can 40 H. ZHANG, H. WANG, AND B. NIU 0.02 0.022 0.024 0.026 0.028 0.03 0.032 0.034 0.28 0.3 0.32 0.34 0.36 q y Positive Equilibrium 0.015 0.02 0.025 0.03 0.1 0.2 0.3 0.4 q y Positive Equilibrium (a) (b) Figure 9. The bifurcation diagrams of the grazer with respect to q. (a) The solution on Poincáre section of system (2.1) with delay and the minimum function in the second equation. (b) The solution of system (2.1) with delay and min { 1, p−y(t− τ) x(t− τ) } = 1. Here the red line is E5 with the dash part being unstable and the solid part being LAS. cause chaotic oscillations of the producer-grazer population by many times binary decisions. This will make it impossible to predict the long-term population trajectories in time. 4. Biological applications In the previous sections, we have observed richer dynamics in system (2.1) through theoretical and numerical analysis. In this section, we will describe the implications of some interesting qualitative dynamical behaviors in population prediction. Time delays often destabilize an internal equilibrium and produce regular oscillations in population size for prey-predator systems, see [8, 9, 12]. However, Fig 4(a) shows that, a properly large maturation delay can make the limit cycle vanish and drive the solution to converge to an equilibrium, which implies that a suitable delay can stabilize a prey-predator system. Moreover, for a small delay, the Rosenzweig’s paradox of enrichment induced by light intensity occurs in our model, which is similar to the existing studies on LKE model in [17, 37]. While large delay impedes such phenomenon, and results in the coexistence of two different LAS positive equilibria, see Fig 7. With an intermediate light input, a transition of population size from one coexistence state to another one may occur. It seems that the periodic oscillations are harmful for the ecological balance of a predator-prey sys- tem, but this dynamical behavior corresponds to predictable evolution in population, see [11] and the references therein for more applications. Therefore, analyzing periodic behavior of system (2.1) is an important part in this paper. Somehow surprisingly, in the presence of time delay, numerical analysis on the periodic solutions near the critical line x+y = p exhibits the period-doubling route to chaos, see subsection 3.3. The chaotic behaviour is related to the irregular fluctuations and variability in nature, which can be caused by many environmental factors. In this paper, we find that the stoichiometric producer-grazer system can be chaotic via period-doubling route. Furthermore, such chaos can also be controlled as Fig 9(a) shows that chaotic oscillations will disappear when the minimum P:C ratio in producer continues to decrease. This provides an appropriate interpretation regarding the irregu- lar fluctuations in producer and grazer with stoichiometry. Although there are numerous results on STOICHIOMETRIC PRODUCER-GRAZER MODEL WITH MATURATION DELAY 41 0.1 0.14 0.18 0.22 0.26 0.26 0.28 0.3 0.32 0.34 0.36 x y 0.1 0.15 0.2 0.25 0.3 0.26 0.28 0.3 0.32 0.34 0.36 x y (a) (b) 0.1 0.15 0.2 0.25 0.3 0.26 0.28 0.3 0.32 0.34 0.36 x y 0.1 0.15 0.2 0.25 0.3 0.28 0.3 0.32 0.34 0.36 x y (c) (d) Figure 10. System (2.1) has a stable period-1 solution when (a) q = 0.0315, (b) q = 0.0295. (c) System (2.1) has a stable period-2 solution when q = 0.0275. (d) System (2.1) has a stable period-4 solution when q = 0.0259. The black line is x+y = p. the dynamics of population models with delays, the chaotic dynamics is uncommon, and the insights concerning the onset and control of chaos require further exploration. For our model with a small delay, two types of solutions coexist that converge to a limit cycle, or a LAS positive equilibrium when they have different initial values, see Fig 11. Moreover, the amplitude of the periodic solution increases over delay, and once the periodic solution collides with the critical line x + y = p, some unclear events take place such that the periodic solution vanishes and the solution finally goes to an equilibrium, see Figs 4 and 5. The biological implication behind that may be as follows. The periodic behavior is sensitive to the initial population size, the increased maturation time can lead the change of population size. Thus, with these population sizes being the new initial point, the solution finally converges to a stable constant state that is robust under tiny perturbation. In fact, our model may work as a theoretical interpretation for the switch of plankton abundances between a relatively stable state and oscillations in one year. There have been some reports reflecting such changes in view of different types of time series of plankton abundances. For instance, Fig 5 in [26] shows that abundances of some microplankton like Dinoflagellates, or some zooplankton like Decapods and Pteropods, may oscillate in several months and become relatively stationary in the remaining 42 H. ZHANG, H. WANG, AND B. NIU period. When the plankton species coexist at a constant state, sudden change of some limiting factors can lead the population size to change, for instance, the increase of temperature or light intensity from April to May causes the significant increase of the abundances of Synechococcus, Dinoflagellates and Copepods, see Fig 5 in [26]. It is worth mentioning that the coexistence of a stable constant steady state and a stable limit cycle in plankton systems has been studied in [40], where that is induced by Bautin bifurcation, while in this paper, it is the joint effect of delay and stoichiometry. 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.42 0.46 0.5 0.54 0.58 x y E 2 E 4 Figure 11. The solutions of system (2.1) with different initial values: the blue one has an initial value as (0.2715, 0.4659), and the initial value of the red one is (0.2715, 0.45). Here, K = 0.84, τ = 0.5 and other parameters are the same as those in (3.1). 5. Discussion A series of newly emerging stoichiometric population models have captured biological features of light- and nutrient-dependent species growth, but they usually neglected the time taken for various physiological processes, such as the maturation time of the grazer. Following LKE model formulated in [19] and rigorous analysis of this model in [17, 37], we formulate a DDE model and analyze the impact of the maturation delay on dynamical behavior. Bifurcation analysis not only reproduces similar behavior as those in [19, 17, 37], but also generates more exciting dynamics beyond the ODE results. For instance, delay drives a positive equilibrium to lose its stability to a stable limit cycle whose amplitude increases as delay increases. Further increasing delay allows the limit cycle to collide with the critical line x + y = p, then the solution starting from a neighborhood of an unstable positive equilibrium converges to a different positive equilibrium instead of the stable limit cycle, see Figs 4 and 5. Due to delay and the Liebig’s law of the minimum, there are period-2, period-4 solutions and the period doubling route to chaos, see Figs 9(a) and 10. Through rigorously mathematical analysis, we provide stability conditions for boundary equilibria and existence conditions for positive equilibria. Various bifurcation phenomena are presented including transcritical bifurcation, saddle-node bifurcation and Hopf bifurcation. We also obtain two types of coexistence: between two stable positive equilibria, between a stable positive equilibrium and a stable limit cycle. Taking biologically relevant parameter values, we conduct simulations to demonstrate theoretical results. Bifurcation analysis is further carried out to explain the joint effects of delay and light, delay and phosphorus on asymptotic dynamics. The increase of delay can decrease the survival STOICHIOMETRIC PRODUCER-GRAZER MODEL WITH MATURATION DELAY 43 region of the grazer on τ − K plane. When the system admits a small delay, Rosenzweig’s paradox of enrichment holds for an intermediate light intensity, while such phenomenon disappears for a large delay. Increasing phosphorus does not change sustainable oscillatory behavior, but it can elevate the survival chance of the grazer. Moreover, the bifurcation diagrams of the grazer over maturation time, conversion efficiency and the minimal P:C ratio are sketched to illustrate their impacts on asymptotic dynamics. We further explain periodic solution, chaos and coexistence of a limit cycle and an positive equilibrium form the biological point of view in section 4. These behaviors in population models are not the first studied, and some results have been applied to control the harmful algal bloom [8], manage fisheries [15], and keep the population to evolve in order [30]. We expect that our results can be used in the development of ecological stoichiometry. There are some problems worthy of further study. The first is the non-negativity and boundedness of solutions. In this paper, we only obtain the conditional non-negativity of solutions and the boundedness of x component. Due to the nonsmoothness induced by Liebig’s law of the minimum, we need phase plane fragmentation and parameter space partitioning, but the delay makes phase space to be infinite- dimensional, thus it is extremely challenging to achieve the goal. We pose this as an open mathematical question. Another perspective is a full analysis of Hopf bifurcation and all possible codimension-two bifurcations occurring at positive equilibria on the critical curve. Note that a positive equilibrium on the curve has two different governing equations on two sides of the curve, thus the Fréchet derivative of the functional differential equation at the equilibrium does not exist. Hence, one cannot obtain the linear stability of the equilibrium by the distribution of corresponding eigenvalues. As pointed out in [17], bifurcation phenomena may widely exist near the equilibrium on the critical curve for the stoichiometric models, however, the minimum functions make eigenvalue analysis, normal form method and center manifold theory [14, 10, 13] inapplicable. Therefore, novel mathematical approaches are needed to rigorously treat such special bifurcations that widely exist in many emerging non-smooth dynamical systems. Acknowledgements We are grateful to two anonymous referees for valuable comments and suggestions which greatly improved the original manuscript. Appendix Here, we study the stability of positive equilibria that are not on the critical line y = p − x. To simplify the analysis, let E∗ = (x∗,y∗) be one positive equilibrium of system (2.1). Linearizing system (2.1) at E∗ yields dx(t) dt = x∗Fx(x ∗,y∗)x(t) + x∗Fy(x ∗,y∗)y(t), dx(t) dt = −dy(t) + Qxt(y ∗,x∗,y∗)x(t− τ) + Qyt(y ∗,x∗,y∗)y(t− τ), (5.1) where one of the following cases is true: (i) Fx = c a + x∗ h ′ (x∗), Fy = − c a + x∗ , Qxt = cêay∗ (a + x∗)2 e−µτ, Qyt = cêx∗ a + x∗ e−µτ. (ii) Fx = c a + x∗ h ′ (x∗), Fy = − c a + x∗ , Qxt = − cê(p−y∗)y∗ (a + x∗)2 e−µτ, Qyt = cê(p− 2y∗) a + x∗ e−µτ, and h ′ (x∗) = b(K−a−2x∗) Kc . Therefore, the corresponding characteristic equation is λ2 + a1(τ)λ + a2(τ)λe −λτ + a3(τ) + a4(τ)e −λτ = 0, (5.2) 44 H. ZHANG, H. WANG, AND B. NIU with a1(τ) = d−x∗Fx(x∗,y∗), a2(τ) = −Qyt(y ∗,x∗,y∗), a3(τ) = −dx∗Fx(x∗,y∗), a4(τ) = x ∗[Fx(x∗,y∗)Qyt(y∗,x∗,y∗) −Fy(x∗,y∗)Qxt(y∗,x∗,y∗)]. When τ = 0, Eq. (5.2) reads λ2 + [a1(0) + a2(0)]λ + a3(0) + a4(0) = 0. (5.3) Obviously, E∗ is LAS provided that a1(0) + a2(0) > 0 and a3(0) + a4(0) > 0. System (2.1) undergoes a Hopf bifurcation at E∗ as a1(0) + a2(0) = 0 and a3(0) + a4(0) > 0. In what follows, we assume that a1(0) + a2(0) > 0 and a3(0) + a4(0) > 0 always hold and investigate the stability change of positive equilibria induced by τ. One can check that λ = 0 is not a root of (5.2) for all τ > 0. Assume that λ = iω (ω > 0) is a root of (5.2) for some τ > 0. Substituting it into (5.2) and separating the real and imaginary parts, we have cos ωτ = − a1(τ)a2(τ)ω 2 + a4(τ)(a3(τ) −ω2) a24(τ) + a 2 2(τ)ω 2 , sin ωτ = a1(τ)a4(τ)ω −a2(τ)ω(a3(τ) −ω2) a24(τ) + a 2 2(τ)ω 2 . This leads to F(ω,τ) := ω4 + [a21(τ) −a 2 2(τ) − 2a3(τ)]ω 2 + [a23(τ) −a 2 4(τ)] = 0. (5.4) Set ∆̃(τ) = [a21(τ) −a 2 2(τ) − 2a3(τ)] 2 − 4[a23(τ) −a 2 4(τ)]. It is simple to see that if (H1) a21(τ) −a 2 2(τ) − 2a3(τ) < 0, a 2 3(τ) > a 2 4(τ) and ∆̃(τ) > 0, is satisfied, then (5.4) has two differential positive roots given by ω±(τ) = 1 √ 2 [ a22(τ) + 2a3(τ) −a 2 1(τ) ± √ ∆̃(τ) ]1/2 . (5.5) Moreover, if a21(τ) −a 2 2(τ) − 2a3(τ) < 0, a 2 3(τ) > a 2 4(τ) and ∆̃(τ) = 0, is true, then ω+(τ) = ω−(τ). If a 2 3(τ) < a 2 4(τ) is true, then only ω+(τ) can exist. Otherwise, (5.4) has no positive roots. Let the set I be I = {τ > 0 : (H1) or a23(τ) < a 2 4(τ) is satisfied}. When I is nonempty and ω(τ) > 0 solves (5.4) for some τ ∈ I, we can define θ ∈ [0, 2π) by cos θ(τ) = − a1(τ)a2(τ)ω 2(τ) + a4(τ)(a3(τ) −ω2(τ)) a24(τ) + a 2 2(τ)ω 2(τ) , sin θ(τ) = a1(τ)a4(τ)ω(τ) −a2(τ)ω(τ)(a3(τ) −ω2(τ)) a24(τ) + a 2 2(τ)ω 2(τ) . Consequently, ω(τ)τ = θ(τ) + 2nπ and iω(τ∗) is a purely imaginary root of (5.2) if and only if τ∗ is the zero of functions Sn(τ) with Sn(τ) = τ − θ(τ) + 2nπ ω(τ) , τ ∈ I, n ∈ N0 = N∪{0}. It can be proved that Sn(τ), n ∈ N0 are continuous and differentiable on I, see [5] for details. Then the result established by [5] is followed. STOICHIOMETRIC PRODUCER-GRAZER MODEL WITH MATURATION DELAY 45 Lemma 5.1. (i) Assume Sn(τ∗) = 0 for some τ∗ ∈ I and n ∈ N0. Then Eq. (5.2) has at least a pair of simple purely imaginary roots λ = ±iω(τ∗). (ii) Let δ(τ∗) = sign { dReλ dτ ∣∣∣ λ=iω(τ∗) } = sign { ∂F(ω(τ∗),τ∗) ∂ω } sign { dSn(τ) dτ ∣∣∣ τ=τ∗ } . This pair of simple purely imaginary roots cross the imaginary axis from left to right if δ(τ∗) > 0 and from right to left if δ(τ∗) < 0. Remark 5.1. Note that Sn(τ) > Sn+1(τ) for all τ ∈ I, which implies that if S0(τ) = 0 has no root in I, then Sn(τ) = 0 has no root in I for all n ∈ N0. Therefore, the real parts of roots of (5.2) remain unchanged for τ ≥ 0. Remark 5.2. A direct calculation yields that at τ = τ∗ ∈ I, ∂F(ω(τ∗),τ∗) ∂ω = ±2ω(τ∗) √ ∆̃(τ∗). Therefore, if ω(τ∗) = ω+(τ∗), then δ(τ∗) = sign { dSn(τ) dτ ∣∣∣ τ=τ∗ } ; if ω(τ∗) = ω−(τ∗), then δ(τ∗) = −sign { dSn(τ) dτ ∣∣∣ τ=τ∗ } . Applying the above results, we first deal with the stability of E2. Through some calculations, we have E2 = (x ∗,h(x∗)), where h(x) = b c ( 1 − x K ) (a + x), x∗ = ad cêe−µτ −d . Moreover, a1(τ) = d− bx∗(K −a− 2x∗) K(a + x∗) , a2(τ) = −d, a3(τ) = − dbx∗ K(a + x∗) (K −a− 2x∗), a4(τ) = dbx∗ K(a + x∗) [ K −a− 2x∗ + a x∗ (K −x∗) ] . It follows that a1(0) + a2(0) = − bx∗(0)(K −a− 2x∗(0)) K(a + x∗(0)) , a3(0) + a4(0) = abd(K −x∗(0)) K(a + x∗(0)) , where x∗(0) = ad cê−d > 0. The following proposition is immediate. Proposition 5.2. At τ = 0, if ad cê−d < K < ad+acê cê−d , then E2 is LAS; if K > ad+acê cê−d , then E2 is a source-type equilibrium; if K < ad cê−d , then E2 is a saddle-type equilibrium. In particular, a Hopf bifurcation occurs at E2 when K = ad+acê cê−d . For τ > 0, we have a21(τ) −a 2 2(τ) − 2a3(τ) = [ bx∗(K −a− 2x∗) K(a + x∗) ]2 > 0, a23(τ) −a 2 4(τ) = ab2d2 (a + x∗)2 ( 1 − x∗ K )[ 4x∗2 + 3ax∗ K − (a + 2x∗) ] . This means that (5.2) has exactly one pair of purely imaginary roots if (H2) K > 4x∗2 + 3ax∗ a + 2x∗ , is satisfied, and all roots of (5.2) have negative real parts if the inequality of (H2) is reversed. 46 H. ZHANG, H. WANG, AND B. NIU Now, we turn to E3 and E4. Denote them as (x∗,h(x∗)), where h(x) = b c ( 1 − x K ) (a + x), and x∗ are two roots of êce−µτ [ p− b c ( 1 − x K ) (a + x) ] −d(a + x) = 0. We similarly obtain a1(τ) = d− bx∗(K −a− 2x∗) K(a + x∗) , a2(τ) = −d + êbe−µτ ( 1 − x∗ K ) , a3(τ) = −dbx∗(K −a− 2x∗) K(a + x∗) , a4(τ) = −bx∗ K(a + x∗) [ d(a + x∗) − êbe−µτ (K −a− 2x∗) ( 1 − x∗ K )] . Then a1(0) + a2(0) = êb ( 1 − x∗(0) K ) − cx∗(0) a + x∗(0) h ′ (x∗(0)), a3(0) + a4(0) = − cx∗(0)h ′ (x∗(0)) a + x∗(0) cêy∗(0) a + x∗(0) − cx∗(0)h ′ (x∗(0)) a + x∗(0) cê(p−y∗(0))y∗(0) (a + x∗(0))2 , = c2êx∗(0)y∗(0) (a + x∗(0))2 [ −h ′ (x∗(0)) − p−y∗(0) a + x∗(0) ] . Due to dx + cêe−µτy = cêe−µτp − ad on line l2, we see p−y∗(0) a+x∗(0) = d cê . At E3, the slope of parabola y = h(x) is larger than zero for all τ ≥ 0, thus, h ′ (x∗(0)) > 0, then a3(0) + a4(0) < 0. Therefore, E3 is unstable when τ = 0. At E4, again using the slope of parabola y = h(x) and line l2, we find h ′ (x∗(0)) < 0 and −h ′ (x∗(0)) − dcê > 0, then a1(0) + a2(0) > 0 and a3(0) + a4(0) > 0 is confirmed. Thus, E4 is always LAS when τ = 0. Considering τ > 0, we have a21(τ) −a 2 2(τ) − 2a3(τ) > [bx∗(K −a− 2x∗) K(a + x∗) ]2 > 0, a23(τ) −a 2 4(τ) = − ( bx∗ a + x∗ )2( 1 − x∗ K )[( 1 − 2x∗ + a K ) êbe−µτ + d ] {( 1 − x∗ K )[( 1 − 2x∗ + a K ) êbe−µτ + d ] − 2d ( 1 − 2x∗ + a K )} . Noticing that 2x∗ + a > dK+êbKe−µτ êbe−µτ , we have( 1 − 2x∗ + a K ) êbe−µτ + d < 0, and a3(τ) 6= a4(τ). This shows λ = 0 is not an eigenvalue and E3 remains the saddle-type of equilibrium. In addition, (5.2) has a pair of purely imaginary roots if and only if (H3) ( 1 − x∗ K )[( 1 − 2x∗ + a K ) êbe−µτ + d ] − 2d ( 1 − 2x∗ + a K ) < 0. We define I1 = {τ ≥ 0 : E2 exists, (H2) is true}, I2 = {τ ≥ 0 : E4 exists, (H3) is true}, and conclude the stability of positive equilibria E2, E3 and E4 as follows. Theorem 5.3. (i) Assume ad cê−d < K < ad+acê cê−d . STOICHIOMETRIC PRODUCER-GRAZER MODEL WITH MATURATION DELAY 47 (i-1) If I1 is empty or S0(τ) = 0 has no positive root in I1 (is nonempty), then E2 is LAS for τ ≥ 0. (i-2) For n ∈ N0, denote Jn = {τn : τn ∈ I1,Sn(τn) = 0}. Assume dSn(τ) dτ ∣∣∣ τ=τn 6= 0 and Jn1 ∩Jn2 = ∅ for n1 < n2. Collect these roots in the set J = ⋃ n∈N0 Jn = {τ0,τ1, ...,τm} with τi < τi+1, 0 ≤ i ≤ m− 1. E2 is LAS when τ ∈ ([0,τ0)∪(τm,∞))∩I1, and is unstable when τ ∈ (τ0,τm). When τ = τi ∈ J, system (2.1) takes place a Hopf bifurcation at E2. (ii) E3 is unstable whenever it exists. (iii) If I2 is empty or S0(τ) = 0 has no positive root in I2 (is nonempty), then E4 is LAS for τ ≥ 0. Moreover, the statement of (i2) is true for E4 when replacing I1 by I2. We take E2 as an example to show the stability change over τ. Still use the parameter set (3.1), by calculation, ad cê−d = 0.2648 and ad+acê cê−d = 0.7797. The proposition 5.2 indicates that for τ = 0, E2 is LAS as 0.2648 < K < 0.7797 and system (2.1) possesses periodic solutions near E2 as K > 0.7797. As τ > 0, it can be checked that condition (H2) is satisfied at E2, and then the graph of Sn on I1 is drawn in Fig 12(a). We see that S0 has two zeros: τ + 0 = 1.2367 and τ − 0 = 21.3813, and Sn, n ∈ N has no zeros, so J = J0 = {1.2367, 21.3813}. In addition, we have dS0(τ) dτ ∣∣∣ τ=τ + 0 = 0.9059 > 0, dS0(τ) dτ ∣∣∣ τ=τ − 0 = −14.1271 < 0. Accordingly, when τ ∈ [0,τ+0 ), E2 is LAS; when τ ∈ (τ + 0 ,τ − 0 ), E2 is unstable; when τ ∈ [τ − 0 , τ̂), τ̂ is the zero of K− 4x ∗2(τ)+3ax∗(τ) a+2x∗(τ) , E2 regains its stability; when τ = τ ± 0 , system (2.1) occurs Hopf bifurcations. The details are illustrated in Fig 12(b), here, we show the change of amplitudes of solutions over τ. 0 5 10 15 20 25 −50 −40 −30 −20 −10 0 10 20 τ S n S 0 S 1 τ 0 + τ 0 − 0 5 10 15 20 25 30 0.35 0.4 0.45 0.5 0.55 τ y τ 0 + τ 0 − (a) (b) Figure 12. (a) The graph of Sn and τ on I1. (b) When K = 0.7, the amplitude of the solution starting from (0.25, 0.4739). E2 is LAS for τ ∈ [0,τ+0 ); there is a stable periodic solution for τ ∈ (τ+0 ,τ − 0 ); E2 is LAS for τ > τ − 0 . Here, the blue and red curves represent the maxima and minima of amplitude of periodic solution, respectively; the magenta curve is the component y in E2 with solid part standing for stable E2 and dash part standing for unstable E2. 48 H. ZHANG, H. WANG, AND B. NIU References [1] P. A. Abrams and J. D. Roth, The effects of enrichment on three-species food chains with nonlinear functional responses, Ecology 75 (1994), 1118-1130. [2] J. F. M. Al-Omari and S. A. Gourley, Monotone travelling fronts in an age-structured reaction-diffusion model of a single species, J. Math. Biol. 45 (2002), 294-312. [3] L. Asik, M. Chen, and A. Peace, The effects of excess food nutrient content on a tritrophic food chain model in the aquatic ecosystem, J. Theor. Biol. 491 (2020), 110183. [4] M. Banerjee and Y. Takeuchi, Maturation delay for the predators can enhance stable coexistence for a class of prey-predator models, J. Theor. Biol. 412 (2017), 154-171. [5] E. Beretta and Y. Kuang, Geometric stability switch criteria in delay differential systems with delay dependent parameters, SIAM J. Math. Anal. 33 (2002), 1144-1165. [6] B. J. M. Bohannan and R. E. Lenski, The effect of resource enrichment on a chemostat community of bacteria and phage, Ecology 78 (1997), 2303-2315. [7] M. T. Brett and D. C. Müller-Navarra, The role of highly unsaturated fatty acids in aquatic foodweb processes, Freshwater Biol. 38 (1997), 483-499. [8] J. Chattopadhyay, R. R. Sarkar, and A. Abdllaoui, A delay differential equation model on harmful algal blooms in the presence of toxic substances, IMA J Math Appl Med Biol. 19 (2002), 137-161. [9] Y. Du, B. Niu, and J. Wei, Two delays induce Hopf bifurcation and double Hopf bifurcation in a diffusive Leslie- Gower predator-prey system, Chaos 29 (2019), 013101. [10] T. Faria and L. T. Magalhães, Normal forms for retarded functional-differential equations with parameters and applications to Hopf bifurcation, J. Diff. Eqns. 122 (1995), 181-200. [11] G. F. Fussmann, S. P. Ellner, K. W. Shertzer, and N. G. Hairston JR, Crossing the Hopf bifurcation in a live predator-prey system, Science 290 (2000), 1358-1360. [12] S. A. Gourley and Y. Kuang, A stage structured predator-prey model and its dependence on maturation delay and death rate, J. Math. Biol. 49 (2004), 188-200. [13] J. K. Hale and S. M. Verduyn Lunel, Introduction to functional-differential equations, Appl. Math. Sci., vol. 99, Springer-Verlag, New York, 1993. [14] B. D. Hassard, N. D. Kazarinoff, and Y. H. Wan, Theory and applications of Hopf bifurcation, London Math. Soc. Lecture Note Ser., vol. 41, Cambridge University Press, Cambridge-New York, 1981. [15] V. Jiménez López and E. Liz, Destabilization and chaos induced by harvesting: insights from one-dimensional discrete-time models, J. Math. Biol. 82 (2021), 1-28. [16] Y. Kuang, J. D. Nagy, and S. E. Eikenberry, Introduction to mathematical oncology, Chapman & Hall/CRC Math. Comput. Biol. Ser. Press, Boca Raton, FL, 2016. [17] X. Li, H. Wang, and Y. Kuang, Global analysis of a stoichiometric producer-grazer model with Holling type functional responses, J. Math. Biol. 63 (2011), 901-932. [18] Z. Liu and N. Li, Stability and bifurcation in a predator-prey model with age structure and delays, J. Nonlinear Sci. 25 (2015), 937-957. [19] I. Loladze, Y. Kuang, and J. J. Elser, Stoichiometry in producer-grazer systems: linking energy flow with element cycling, Bull. Math. Biol. 62 (2000), 1137-1162. [20] R. M. May, Biological populations with nonoverlapping generations: stable points, stable cycles and chaos, Science 186 (1974), 645-647. [21] J. A. J. Metz and O. Diekmann, The dynamics of physiologically structured populations, Lect. Notes in Bioinform., vol. 68, Springer-Verlag, Berlin, 1986. [22] S. J. Moe, R. S. Stelzer, M. R. Forman, W. S. Harpole, T. Daufresne, and T. Yoshida, Recent advances in ecological stoichiometry: insights for population and community ecology, Oikos 109 (2005), 29-39. [23] E. Ott, Chaos in dynamical systems, second ed., Cambridge University Press, Cambridge, 2002. [24] A. Peace, H. Wang, and Y. Kuang, Dynamics of a producer-grazer model incorporating the effects of excess food nutrient content on grazer’s growth, Bull. Math. Biol. 76 (2014), 2175-2197. [25] M. Rehim and M. Imran, Dynamical analysis of a delay model of phytoplankton-zooplankton interaction, Appl. Math. Model. 36 (2012), 638-647. [26] J. B. Romagnan, L. Legendre, L. Guidi, J. L. Jamet, D. Jamet, L. Mousseau, M. L. Pedrotti, M. Picheral, G. Gorsky, C. Sardet, and L. Stemmann, Comprehensive model of annual plankton succession based on the whole-plankton time series approach, PLoS One 10 (2015), e0119219. [27] M. L. Rosenzweig, Paradox of enrichment: destabilization of exploitation ecosystems in ecological time, Science 171 (1971), 385-387. STOICHIOMETRIC PRODUCER-GRAZER MODEL WITH MATURATION DELAY 49 [28] J. W.-H. So, J. Wu, and X. Zou, A reaction-diffusion model for a single species with age structure. I travelling wavefronts on unbounded domains, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 457 (2001), 1841-1853. [29] R. V. Solé, J. G. P. Gamarra, M. Ginovart, and D. López, Controlling chaos in ecology: from deterministic to individual-based models, Bull. Math. Biol. 61 (1999), 1187-1207. [30] L. Stone, Period-doubling reversals and chaos in simple ecological models, Nature 365 (1993), 617-620. [31] I. V. Telesh, H. Schubert, K. D. Joehnk, R. Heerkloss, R. Schumann, M. Feike, A. Schoor, and S. O. Skarlato, Chaos theory discloses triggers and drivers of plankton dynamics in stable environment, Sci. Rep. 9 (2019), 20351. [32] J. Urabe and R. W. Sterner, Regulation of herbivore growth by the balance of light and nutrients, Proc. Natl. Acad. Sci. USA. 93 (1996), 8465-8469. [33] H. Wang, K. Dunning, J. J. Elser, and Y. Kuang, Daphnia species invasion, competitive exclusion, and chaotic coexistence, Discrete Contin. Dyn. Syst. Ser. B 12 (2009), 481-493. [34] H. Wang, Y. Kuang, and I. Loladze, Dynamics of a mechanistically derived stoichiometric producer-grazer model, J. Biol. Dyn. 2 (2008), 286-296. [35] H. Wang, Z. Lu, and A. Raghavan, Weak dynamical threshold for the ”strict homeostasis” assumption in ecological stoichiometry, Ecol. Model. 384 (2018), 233-240. [36] H. Wang, R. W. Sterner, and J. J. Elser, On the ”strict homeostasis” assumption in ecological stoichiometry, Ecol. Model. 243 (2012), 81-88. [37] T. Xie, X. Yang, X. Li, and H. Wang, Complete global and bifurcation analysis of a stoichiometric predator-prey model, J. Dyna. Diff. Eqns. 30 (2018), 447-472. [38] X. Yang, X. Li, H. Wang, and Y. Kuang, Stability and bifurcation in a stoichiometric producer-grazer model with knife edge, SIAM J. Appl. Dyn. Syst. 15 (2016), 2051-2077. [39] S. Yuan, D. Wu, G. Lan, and H. Wang, Noise-induced transitions in a nonsmooth producer-grazer model with stoichiometric constraints, Bull. Math. Biol. 82 (2020), 1-22. [40] H. Zhang and B. Niu, Dynamics in a plankton model with toxic substances and phytoplankton harvesting, Internat. J. Bifur. Chaos Appl. Sci. Engrg. 30 (2020), 2050035. Hua Zhang, Department of Mathematics, Harbin Institute of Technology, Weihai, Shandong, 264209, P. R. China Email address: zhanghua.math@foxmail.com Hao Wang, Department of Mathematical and Statistical Sciences, University of Alberta, Edmonton, Alberta, T6G 2G1, Canada Email address: hao8@ualberta.ca Ben Niu, Corresponding author, Department of Mathematics, Harbin Institute of Technology, Weihai, Shandong, 264209, P. R. China Email address: niu@hit.edu.cn 1. Introduction 2. Basic analysis 2.1. Nonnegativity and boundedness 2.2. Stability of boundary equilibria 2.3. Existence and stability of positive equilibria 3. Numerical simulations 3.1. Biological mechanisms related to 3.2. The bifurcation diagrams over K or 3.3. Period doubling route to chaos 4. Biological applications 5. Discussion Acknowledgements Appendix References