Dynamics of a stoichiometric producer-grazer model with maturation delay

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 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.


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 zooplanktonphytoplankton 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.
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, 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 Note that lim x→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 (y j ) 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 follows: dx(t) dt = bx(t) 1 − x(t) min{K, (P − θy(t))/q} − f (x(t))y(t), dy j (t) dt =ê min 1, P − θy(t) θx(t) f (x(t))y(t) − µy j (t) −ê min 1, 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 and the mature grazer density with w s (s) = y(s). This leads to w s (t) = e −µ(t−s) y(s). Therefore, 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, (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 stoichiometry 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 intensity 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 stoichiometric 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.
(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 Rosenzweig'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.

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 (2.1) 2.1. Nonnegativity and boundedness. In the biological perspective, the initial conditions are given as Denote k = min{K, P/q}, the basic property of solutions with initial values (2.2) is stated in the following theorem.
, y(t)) be any solution of system (2.1) subject to initial conditions (2.2). Then, Proof. Solving x(t) from the first equation of (2.1) gives Thus, x(t) ≥ 0 for all t > 0, and x(t) > 0 if x(0) > 0. Moreover, we see that From the standard comparison argument, we have lim sup t→∞ x(t) ≤ k.
We can repeat the process for any [nτ, (n + 1)τ ], n ≥ 2. 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: E 0 = (0, 0) and E 1 = (k, 0), and their local stability is obtained by the distribution of eigenvalues corresponding to the linearized system. (2.4) this suggests that two eigenvalues have different signs, so E 0 is always unstable. At E 1 , if k < p, then the characteristic equation takes the following form Note that the stability of E 1 depends on the real parts of the zeros of λ + d −êf (k)e −µτ e −λτ . When 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 Similarly, we obtain E 1 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 As a consequence, the positive equilibria are determined by (2.7) When the parabola y = h(x) and the line l 1 intersect in the region {(x, y) ∈ R 2 : 0 < x < K, 0 < y, x + y < p}, there exists a positive equilibrium, denoted by E 2 . As the parabola y = h(x) enters the region {(x, y) ∈ R 2 : 0 < x < K, 0 < y, x + y > p}, when it is tangent with the line l 2 , a newly positive equilibrium exists; when it intersects with the line l 2 at two different points, system (2.1) has two newly positive equilibria, denoted by E 3 and E 4 . 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: (ii) If the line l 1 , the parabola y = h(x) and the critical line , where x * is given in (2.7). (iii) If the parabola y = h(x) is tangent to the line l 2 at (x,ỹ), then we have (iv) Denoting the intersection point of the line l 2 and x−coordinate as K 4 , we have K 4 (τ ) = cêp d e −µτ − a. One can check that if τ satisfies assumption (A), then 0 < K 1 (τ ) < p < K 4 (τ ) for each τ . If cp < ab or cp > ab,êe −µτ < ad+pd cp−ab is further true, then 0 . Based on [37, Theorems 3.1-3.7], we state the existence of positive equilibria in the following theorem.
Then system (2.1) has no positive equilibria.
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 I 1 , I 2 and S n (τ ) are defined in Appendix.

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 Hopf bifurcation occurs if S n (τ ), n ∈ N 0 has zeros in I 1 Hopf bifurcation occurs if S n (τ ), n ∈ N 0 has zeros in I 2 Hopf bifurcation occurs if S n (τ ), n ∈ N 0 has zeros in I 2 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 The pair (τ, K) satisfying K 1 (τ ) < K < K 2 (τ ) is in regions D3, D4, D5 and D6, thus E 2 exists in these regions. We further know that in regions D3 and D4, E 2 is unstable induced by Hopf bifurcation, and it remains LAS in regions D5 and D6.
When K + 3 (τ ) < K < K 2 (τ ), (τ, K) locates in regions D3 and D5, and E 3 exists in these regions. When K + 3 (τ ) < K < K 4 (τ ), (τ, K) belongs to regions D2, D3 and D5, E 4 exists in these regions. All bifurcation curves described in Remark 2.2 are also exhibited in Fig 1(a). Observe that three curves: K = K 1 (τ ), K = K 2 (τ ) and K = K 4 (τ ) intersect at P 3 , at which two positive equilibria with one merged by E 2 and E 3 and the other E 4 collide with E 1 from its two sides, then all positive equilibria disappear and only one boundary equilibrium stays. In fact, the position of P 3 on τ − K plane is (τ , K 1 (τ )), whereτ = 1 µ ln cêp d(a+p) . It is also seen that the curve K = K + 3 (τ ) is tangent with the curve K = K 2 (τ ) at P 2 , thus E 2 , E 3 and E 4 collide into one positive equilibrium. Meanwhile, it can be calculated that there is a zero eigenvalue for the corresponding characteristic equation. Therefore, we assert P 2 is a cusp bifurcation point. The Hopf bifurcation curve near E 2 intersects K-axis at (0, 0.7797) and K = K 2 (τ ) at P 1 , respectively. At point P 4 , one can see that E 3 and E 4 collide into one equilibrium, E 2 goes to a boundary equilibrium and disappears. As a result, system (2.1) has exactly one positive equilibrium.
Biologically, K = K 1 (τ ) (the red curve) and K = K 4 (τ ) (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 K 1 (τ ) is an increasing function and K 4 (τ ) 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  Figure 1. The bifurcation diagrams on τ − K plane. The red curve is determined by K = K 1 (τ ); the blue curve stands for K = K 2 (τ ); the magenta curve is given by K = K 4 (τ ); the green one is K = K + 3 (τ ); the black one is Hopf bifurcation curve at E 2 . quantity of the producer which drives the grazer to go extinct as well. Moreover, the existence of P 3 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 = K 2 (τ ) at two distinct points, denoted by P 0 and P 1 . 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 E 2 .
In view of Fig 1(a) and (b), system (2.1) undergoes a saddle-node bifurcation when E 2 and E 3 (or E 3 and E 4 ) collide into one positive equilibrium. This phenomenon is presented in  , the variation of the grazer density is complex. The grazer density slowly reduces until τ = 1.278 where a supercritical Hopf bifurcation occurs at E 2 such that E 2 loses its stability, and the grazer density changes periodically for 1.278 < τ < 21.308. At τ = 21.308, a periodic solution disappears and E 2 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,  K = 0.585185 is an important threshold value, at which E 2 is unstable and system (2.1) without delay has at least one limit cycle around E 2 , and E 3 and E 4 merge into a saddle-node equilibrium E 3,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, E 3 and E 4 leave the critical line x + y = p and separate, E 3 is unstable and E 4 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 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   It can be observed that system (2.1) admits the unique positive equilibrium E 2 as τ and p vary simultaneously. Moreover, when parameter values are in the region between two black vertical lines, E 2 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,  3.2. The bifurcation diagrams over K orê. In this part, we study the relationship between dynamical 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, E 2 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 E 3 and stable E 4 . As K further increases, the grazer density keeping on E 4 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 E 2 and E 3 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 E 2 , then two new equilibria, E 3 and E 4 , appear that coexist with E 2 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, E 2 collides with E 3 and then disappears. The grazer continues to persist until K = 0.9501. 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 E 2 and E 3 . 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 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. 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: E 2 , E 3 and E 4 , 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 E 3 and stable E 4 disappear simultaneously through a saddle-node bifurcation. Asê further increases, the grazer density rises steadily untilê = 0.5726, where stable E 2 loses its stability to a periodic solution. The amplitude of the periodic solution increases aŝ e 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. Fixing τ = 10, we find that when 0.0125 < q < 0.1490, there is at least one positive equilibrium E 5 = (x 5 , y 5 ) with x 5 = 0.1677, and Moreover, system (2.1) takes place the chaos routed by period doubling near E 5 as q decreases in the interval [0.0238, 0.034], see Fig 9(a). Here we choose the component x of E 5 as a Poincáre section, and draw the change of the solution on the Poincáre section over q. When q > 0.032, E 5 is LAS. When 0.0288 < q < 0.032, E 5 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, E 5 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 E 5 , see Fig 9(b). Here, the maxima and minima of the solutions are drawn. It can be seen that E 5 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 Here the red line is E 5 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.

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 system, 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 irregular fluctuations in producer and grazer with stoichiometry. Although there are numerous results on 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 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.

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 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 infinitedimensional, 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.
It can be proved that S n (τ ), n ∈ N 0 are continuous and differentiable on I, see [5] for details. Then the result established by [5] is followed.
Remark 5.1. Note that S n (τ ) > S n+1 (τ ) for all τ ∈ I, which implies that if S 0 (τ ) = 0 has no root in I, then S n (τ ) = 0 has no root in I for all n ∈ N 0 . Therefore, the real parts of roots of (5.2) remain unchanged for τ ≥ 0.
Proposition 5.2. At τ = 0, if ad cê−d < K < ad+acê cê−d , then E 2 is LAS; if K > ad+acê cê−d , then E 2 is a source-type equilibrium; if K < ad cê−d , then E 2 is a saddle-type equilibrium. In particular, a Hopf bifurcation occurs at E 2 when K = ad+acê cê−d . For τ > 0, we have 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.
(i-1) If I 1 is empty or S 0 (τ ) = 0 has no positive root in I 1 (is nonempty), then E 2 is LAS for τ ≥ 0.
(ii) E 3 is unstable whenever it exists.
(iii) If I 2 is empty or S 0 (τ ) = 0 has no positive root in I 2 (is nonempty), then E 4 is LAS for τ ≥ 0. Moreover, the statement of (i2) is true for E 4 when replacing I 1 by I 2 .