Stoichiometric food chain model on discrete time scale

Stoichiometry-based models can yield many new insights into producer grazer systems. Many interesting results can be obtained from models continuous in time. There raises the question of how robust the model predictions are to time discretization. A discrete stoichiometric food-chain model is analyzed and compared with a corresponding continuous model. Theoretical and numerical results show that the discrete and continuous models have many properties in common but differences also exist. Stoichiometric impacts of producer nutritional quality also persist in the discrete system. Both types of models can exhibit qualitatively different behaviors with the same parameter sets. Discretization enlarges the parameter ranges for the existence of chaotic dynamics. Our results suggest that the stoichiometric mechanisms are robust to time discretization and the nutritional quality of the producer can have dramatic and counterintuitive impacts on population dynamics, which agrees with the existing analysis of pelagic systems.


Introduction
Recent advances towards the understanding of ecological interactions have been made through the development of the theory of ecological stoichiometry.Ecological stoichiometry is the study of the balance of energy (carbon) and multiple chemical elements (phosphorus, nitrogen, etc.) in ecological interactions [13].Most ecological studies have until very recently ignored the sources and consequences of this chemical heterogeneity.The stoichiometric ratios of elements, such as carbon and phosphorus, have complex dynamics as they vary within and across trophic levels.Using the fact that all organisms are structurally composed of multiple chemical elements combined in non-arbitrary proportions, it is necessary to stress the importance of incorporating the effects of food quality, as well as quantity into ecological modeling.
Recent efforts towards understanding predator-prey dynamics under nutritional constraints have been made with the development of stoichiometric population models [7,8,9,12,16].Loladze et al. 2000 [8] formulated a tractable two-dimensional Lotka-Volterra type model (LKE model) that incorporates stoichiometry into the transfer of elements between producer and grazer.The LKE model provided a foundation that many stoichiometric models built upon.It utilizes the fact that both producer and grazer are chemically heterogeneous organisms.Specifically, it tracks the amount of two essential elements, C (Carbon) and P (Phosphorus), in each trophic level.Extending the LKE model to explicitly track environmental nutrient supplies, Kuang et al. 2004 [7] mechanistically formulated a tractable model of plant-herbivore population dynamics (KHE model).An additional extension of the LKE model was developed by Loladze et al. 2004, which models two competing predators on one prey (LKEF model).All these stoichiometric models show that an increase in prey quantity accompanied by decrease in its quality can halt oscillations in predator-prey system.These effects are due to the decreased energy flow from the prey to the predator.In other words, though the predator consumes more prey, it effectively converts less of it into its own biomass [5].
It is widely recognized that the selection of time scale is very important in the study of biology and ecology [6].The stoichiometric models mentioned above are continuous in time.This raises the question about the sensitivity of the results to the discretization of time.Fan et al. [5], Sui et al. [14] and Xie et al. [17] compared the continuous stoichiometric LKE producer-grazer model [8], KHE plant-herbivore model [7] and LKEF two predators-prey model [9] to their discrete analogues, respectively.They systematically investigated the dynamic behaviors of the discrete models and conducted some comparisons with their corresponding continuous models.The results show that there are many similarities between the continuous models and their discrete analogues.However, the discrete models also exhibit new phenomenons, like chaos, not observed in the corresponding continuous models.
Peace 2015 [11] proposed a stoichiometric model in order to investigate how stoichiometric constraints affect food-web dynamics and nutrient cycling in ecosystems.Another extension of the LKE model, this is a linear food chain with three trophic levels subject to varying carbon (light) and environmental phosphorus levels.A comprehensive global analysis of the rich dynamics (i.e., chaos) of this model is explored both analytically and numerically by Chen et al. 2017 [3].Note that the food chain model proposed in [11] is defined on the continuous time scale, while experimental data are usually discretely collected, many producers such as annual plants have non-overlapping generations, and many herbivores have seasonal dynamics [17].This also raises the question about the sensitivity of the main findings of this stoichiometric food chain model to the discretization of time.
The objective of this paper is to compare the dynamics of the continuous stoichiometric tritrophic food chain model [11] with its discrete analog.In particular, we focus on these three questions: (1) whether the following most important features observed in the experiments and exhibited in the continuous stoichiometric model 'survive' the discretization, like deterministic extinction of consumers and predators when faced with low quality food; (2) whether chaotic dynamics can arise under a biologically plausible parameter set and if it does, whether the parameter ranges of the discrete model admitting chaos are the same with the continuous one; (3) whether the discrete model exhibit new mechanisms not observed by the continuous model.
In the next three sections, we construct a discrete analogue of the stoichiometric tritrophic food chain model and explore its dynamics qualitatively.In section 5, we numerically study our discrete food chain model with Holling type II functional response and systematically expound the similarities and differences between our discrete food chain model and its continuous analogue.

The discrete model
Peace 2015 [11] proposed a three trophic level model incorporating stoichiometric principles that takes the following form: where x is the density of producer (in milligrams of carbon per liter, mgC/L), y is the density of consumer (mgC/L), z is the density of predator (mgC/L), b is the intrinsic growth rate of producer (day −1 ).Parameters d y and d z are the specific loss rate of consumer and predator respectively that includes metabolic losses (respiration) and death (day −1 ).Parameters êy and êz are constant production efficiencies (yield constant) of the consumer and predator, respectively.K is the producer's constant carrying capacity that depends on some external factors such as light intensity.Functions f (x) and g(y) are the consumer's and predator's ingestion rate, which may be a Holling type II functional responses.
To better understand our discrete model formulation and its comparison with model (2.1), it is necessary to recall here some main assumptions of the continuous model.They are listed below: A1: The total mass of phosphorus in the entire system is fixed, i.e., the system is closed for phosphorus with a total of P T (mgP/L).A2: Producer's P (phosphorus) : C (carbon) varies, but never falls below a minimum q (mgP/mgC); consumer and predator maintain a constant P : C, θ y and θ z (mgP/mgC), respectively.
A3: All phosphorus in the system is divided into three pools: P in the consumer, P in the predator and P in the producer.
Note that assumption A2 allows the producer stoichiometric P:C ratio to vary but holds the ratios higher trophic levels fixed.This strict homeostatic assumption is based on the fact that, although consumer and predator stoichiometries are variable, the range of variation is small compared to the range of producer stoichiometries.Assumption A3 does not allow for P to be dissolved in the environment, but assumes all available P for the producer is taken up immediately.
Our discrete food chain stoichiometric model is based on the above continuous time model (2.1).There are several ways of deriving discrete time versions of dynamical systems corresponding to a continuous time formulation.Here we follow the method used in [5].

Analysis
First, we prove some preliminary results for (2.4), such as boundedness and positive invariance.It is easy to see that all components of the solution of (2.4) stay positive when they exist.By carrying out similar arguments like those in Theorem 3.1 and 3.2 in [5], we reach the following theorems.In Loladze [8], it is assumed that the function f (x) in system (2.4) is a bounded smooth function satisfying Also we assume the function g(y) has the same properties that g(0) = 0, g (y) > 0, g (0) < ∞ and g (y) < 0 for y > 0.
In this paper, for convenience we assume that f (x) = xp(x) and g(y) = ys(y).It has been proven in [8] that lim Theorem 3.1.For system (2.4), we have for all n > 0, where v and w are any numbers satisfying êy f (e b−p(U)v ) < d y and êz g(Ve êy −s(V)w−d y ) < d z .

Mathematical Biosciences and Engineering
Volume 16, Issue 1, 101-118 Proof.It is easy to see that the both components of the solution of equation (2.4) stay positive when exist.Notice that the function F(x) = x exp b(1 − x/K) achieves its maximum at x = K/b, we have Hence, for all nonnegative integers n, we have If êy f (U) ≤ d y , then it is clear that for all n > 0, we have y(n) ≤ y(0).We thus assume below that êy f (U) > d y .Let v be large enough so that êy f (e b−p(U)v ) < d y .
We claim that for all nonnegative integers n, we have

This is obviously true for
Assume the above claim is false.Assume first the case that (I): and y(n 1 ) > V.In this case, we have (using the property of p (x) < 0) a contradiction, indicating the claim is true for this case.
(II): Assume now the case of y(0) > v.In this case, we have which implies that y(2) < y(1).
In other words, as long as y(n) > v, then y(n + 2) < y(n + 1).Hence there are two possibilities, either (i) for some y(n * ) ≤ v for some n * > 0; or (ii) y(n) > v for all n > 0. In case of (ii), we have y(n) strictly decreasing for n > 1 and the claim is obviously true.In case of (i), from the proof of case (I), we see that for y(n) < V for n > n * and hence the claim is also true.We can also use the same method to discuss the boundedness of z(n).Let w be large enough so that êz g(Ve êy −s(V)w−d y ) < d z .
We claim that for all nonnegative integers n, we have A straightforward consequence of the above theorem is that is positively invariant.Indeed, we have Theorem 3.2.For (2.4), ∆ is globally attractive with respect to initial values x(0), y(0) and z(0) such that x(0) > 0, P T /θ y > y(0) > 0 and P T /θ z > z(0) > 0.
Proof.It is easy to see from the model equations and the proof of the above theorem that for large values of n, x(n) ∈ (0, (K/b)e b(1−1/b) ).Notice that if y(n) > v for large values of n > 0, then y(n) must have a limit y * ≥ v. Hence for large values of n, we have This contradicts y * > v > 0. Then for large values of n, y(n) ∈ (0, v).On the other side, notice that if z(n) > w for all n > 0, then z(n) must have a limit z * ≥ w.Hence for large values of n, we have This contradicts z * > w > 0, proving the theorem.

Equilibria
To facilitate the equilibria analysis in the following discussion, we rewrite equation (2.4) as where To find equilibria of equation (2.4) we solve The variational matrix of equation (2.4), after some straightforward algebraic calculations, reads where H z = ∂H ∂x = 0. (4.9)

Boundary equilibria
The possible equilibria or steady states of equation (2.4) solve the system of algebraic equations It is easy to observe that the equilibria of equation (2.4) are exactly the same as those of equation (2.1).
The only boundary equilibria are E 0 (0, 0, 0), E 1 (k, 0, 0) and E 2 ( x, ȳ, 0).At the origin E 0 (0, 0, 0), the Jacobian matrix takes the form It is clear that the two characteristic roots e −d y and e −d z have magnitude less than 1 while the other has magnitude e b larger than 1.Consequently, the origin E 0 is always unstable.
Proof.If we let λ 1 , λ 2 and λ 3 be the characteristic roots of J(E 1 ), then the condition 0 < b < 2 and êy min(1, implies λ i > 1 for some i. Remark 1. From [8], we know that G(x, 0, 0) = 0 has a unique solution x 1 ∈ (0, P T /θ y ) and a unique solution ] then E 1 is stable.But for equation (2.4), the situation is quite different since E 1 can be possibly stable or possibly unstable when k x), then one can easily derive sufficient conditions for the stability of E 0 and E 1 in equation (2.4).
The boundary equilibrium E 2 ( x, ȳ, 0) indicates the extinction of predator.Thus E 2 can be viewed as the internal equilibria of the degenerating simpler two-dimensional system with only one consumer on one producer.E 2 ( x, ȳ, 0) can be studied in the x-y phase plane.Since both equations (2.1) and (2.4) have the same equilibria, from [5] we claim that equation (2.4) possibly has multiple E 2 -type boundary equilibria.In the following, we assume that E 2 is such an equilibrium and discuss its local stability.
The Jacobian at the predator-free equilibria E 2 ( x, ȳ, 0) reads We have the following theorem on the local asymptotical stability of the boundary equilibrium E 2 by referenced [5]., then in region i (i.e., ȳ < P T /θ y − x), the following are true.
• If the producer's nullcline is decreasing at E 2 (i.e., F x > 0), and then E 2 is locally asymptotically stable.
• If the producer's nullcline is decreasing at E 2 (i.e., F x > 0), and In region ii (i.e., ȳ < P T /θ y − x), the following are true.• If the slope of the producer's nullcline at E 2 is greater than the consumer's (i.e., −G x /G y < −F x /F y ), then E 2 is unstable.
• If the slope of the consumer's nullcline at E 2 is greater than the producer's (i.e., −G x /G y > −F x /F y ), and then E 2 is locally asymptotically stable.
• If the slope of the consumer's nullcline at E 2 is greater than the producer's (i.e., −G x /G y > −F x /F y ), and then E 2 is unstable.

Internal equilibria
If system (2.4) has an internal equilibrium E * (x * , y * , z * ), then it satisfies To determine the type of species interactions occurring at this equilibrium, it is convenient to examine the Jacobian matrix at E * taking the following form, as the ecosystem matrix of (2.4) at E * .Due to the complexity of the system (2.4), we can not study the stability of positive equilibrium routinely by Routh-Hurwitz criteria, however we can conduct a qualitative analysis of the ecosystem matrix.In this matrix, the i j-th term measures the effect of the j-th species on the i-th species's per capita growth rate.
If food quality is good the variable P:C ratio of the producer is larger than the fixed P:C ratio of the consumer, i.e., (P T − θ y y * − θ z z * )/x * > θ y , and the ecosystem matrix at E * takes the following form , which means that the interaction between the consumer and the producer is of conventional (+, −) type.If food quality is bad the variable P:C ratio of the producer is less than the fixed P:C ratio of the consumer, i.e., (P T − θ y y * − θ z z * )/x * < θ y , and the ecosystem matrix takes the form , which implies that the interaction between the consumer and the producer is changed from the traditional (+, −) to (−, −).In other words, both the producer and the consumer compete with each other and the competition lies in two species not only internally but also interspecifically for the limited available P. Here, increasing the quantity of a producer population of bad quality will not promote the consumer's growth.An unusual and important property of (2.4) is that the sign of ∂G/∂x changes from positive to negative as soon as the producer's P : C is less than that of the consumer's.This negative derivative means that, higher producer density reduces the growth rate of consumer.If the quality of the producer is bad, then any further increase in the producer's density deteriorates the quality of producer, which offsets the benefit provided by the higher producer density to the consumer that has already been limited by P.This effect is in sharp contrast to many conventional population dynamics models.

Numerical simulations
In this section, we carry out systematic numerical simulations to verify and deepen our analytical findings and to compare them with those of the corresponding continuous system.In particular, we will discuss the impact of the discretization of the model (2.1).
In the following, we assume that p(x) = c 1 /(a 1 + x) and s(y) = c 2 /(a 2 + y).We will numerically study the following food chain model with Holling type II functional response In order to compare the dynamical behaviors of the discrete model (5.1) and continuous model (5.2), we carry out numerical simulations with the same parameters used in [3], which are listed in Table 1.These values are mainly selected from [2] and [15] and have been applied by [8] for guidance in setting parameters at biologically realistic values.
The parameter K, indirectly represents the light intensity or energy input into the system.Bifurcation diagrams provide a clear and visual way to understand how the behavior of a dynamical system depends on a parameter.The bifurcation diagrams in Figure 1 illustrate how the dynamics of the populations in equations (5.1) and (5.2) vary with K. Numerical simulations of the time series for both models along the gradient of light levels K from 0 to 10 mgCL −1 are shown in Figure 2 to assist in intuitive analysis.
When the light intensity is low (e.g., 0 < K < 0.15 in Figure 1), the consumers and predators of the discrete and continuous models can not persist due to starvation.When K increases further (e.g., 0.15 < K < 0.5 in Figure 1), the consumer and predator of the two models persist at an stable equilibrium one after the other.Figure 2 (a 1 ) and (b 1 ) shows such a typical case with K = 0.4.As K continues to increase, the stable equilibrium loses its stability to a limit cycle.This limit cycle starts with zero amplitude and grows as K is further increased within a reasonable interval (e.g., 0.5 < K < 0.7 in Figure 1).Both systems (5.1) and (5.2) exhibit the Rosenzweig's paradox enrichment.Figure 2 (a 2 ) and (b 2 ) show periodic oscillating dynamics when K = 0.6.As K further increases, both dynamics of discrete and continuous models seem 'chaotic'.As a typical example, for K = 1, the solution curves in Figure 2 (a 3 ) and (b 3 ) show all populations oscillate irregularly and admit chaos.After that, Figure 1.The bifurcation curves versus K for the discrete model (a i , i = 1, 2, 3) and continuous model (b i , i = 1, 2, 3).The parameters are defined in Table 1 and the initial conditions are x = 0.5 mgCL −1 , y = 0.5 mgCL −1 and z = 0.5 mgCL −1 .populations in both models all come back to regularly periodic status (e.g., Figure 2 (a 4 ) and (b 4 ) with K = 4).As K increases from intermediate levels to high levels (e.g., 2.4 < K < 6.5 in Figure 1), the mean values of the oscillatory density of the producers of two models have a tendency to rise while that of the predators catch an opposite trend due to the bad quality of the producers.The amplitudes of the limit cycle solutions of both models gradually increase from zero, and then decrease to zero again.At a threshold value (e.g., K = 6.5 in Figure 1), the limit cycles of the two models collapse and interior equilibrium E * gains stability.All three species coexist at this equilibrium for a large enough range of K (e.g., 6.5 < K < 9.8 in Figure 1).Figure 2 (a 5 ) and (b 5 ) show case studies with K = 8.Extremely high light intensity (e.g., 9.8 < K < 10 in Figure 1) leads to a catastrophic crash of both predator and consumer due to low food quality of the producer.Figure 2 (a 6 ) and (b 6 ) are case studies with K = 10 where both predator and consumer perish.Due to the extinction of higher trophic levels, the primary level of the food chain (i.e., producer) in both discrete and continuous models tends towards its carrying capacity.This type of scenario is often observed as large algal blooms.
The numerical simulations in Figure 1 and Figure 2 reveal that both the continuous model and its discrete analogue exhibit similar dynamics, but differences exist.There are some parameter sets that yield different dynamics for model (5.1) and (5.2).When K = 1.13, populations of the continuous model tend to a positive stable equilibrium, while those of the discrete model show 'chaotic' behaviors (Figure 2 (a 7 ) and (b 7 )).This phenomenon also contributes to the conclusion drawn by Fan et al 2005 [5] that discretization can destabilize the internal equilibria.
There are some subtle but more important differences between the dynamics of equations (5.1) and (5.2).Both models exhibit chaos but at different parameter intervals (Figure 1).The simulations of the spectra of maximum Lyapunov exponent (MLE) of models (5.1) and (5.2) against K are shown in Figure 3.The maximum Lyapunov exponent of model (5.1) is positive when K varies within the inter- val K ∈ (0.9, 2.4), while that of model (5.2) is positive for K ∈ (0.7, 1.1).The internal equilibrium of equation (5.1) undergoes chaotic behavior later than the corresponding one of (5.2).Furthermore, the parameter range of model (5.1) admitting chaos is bigger than that of model (5.2).That is to say discretization can enlarge the parameter ranges for the existence of chaotic dynamics of the stoichiometric tritrophic models.

Conclusion and discussion
Continuous stoichiometric models can be used to study the balance of energy and multiple chemical elements between populations with overlapping generations.However, experimental data are usually collected on discrete time intervals, many plants in wild and agricultural settings have non-overlapping generations, and many herbivores have clear annual or seasonal dynamics [14].Such scenarios call for discrete models.In this paper, a stoichiometric discrete model with three trophic levels is developed and analyzed.The model is a MacArthur-Rosenweig type food chain model under stoichiometric constraints and has rich dynamics.The possible attractors of model (2.4) include boundary equilibrium, internal equilibrium, limit cycles, or even a strange attractor (chaos) (see Figure 2).The theoretical and numerical results of our discrete model (2.4) show that they retain most of the important dynamical features exhibited by its continuous analogs (2.1). Figure 1 and Figure 3 suggest that the discrete model (2.4) admits the paradox of enrichment and also chaotic dynamics.These are typical dynamics  1 and the initial conditions are x = 0.5 mgCL −1 , y = 0.5 mgCL −1 and z = 0.5 mgCL −1 .
of resource explicit population models, like stoichiometric predator-prey systems [3].The common and notable point is that very low producer quality can drive the consumer and predator to deterministic extinction.Some important differences between the continuous and discrete models exist.The models predict qualitatively different dynamic behaviors for particular parameter sets (Figure 2 (a 7 ) and (b 7 )).These results highlight the fact that discretization can destabilize the internal equilibria.When light levels are low or high (e.g., 0 < K < 0.7 or 2.4 < K < 10 in Figure 1), the discrete model (2.4) and the continuous model (2.1) exhibit similar dynamical behaviors.When light levels are at intermediate levels (e.g., 0.7 < K < 2.4 in Figure 1) different dynamics are predicted.Maximum Lyapunov exponents of models (2.1) and (2.4) in this region (Figure 3) show that the biologically plausible parameter range for chaotic dynamics of the discrete model (2.4) is larger than that of the continuous model (2.1).This supports Andersen's [2] and Fan's [5] argument on the possibility of chaos in discrete stoichiometric systems.
In fact, the discretization can significantly impact the dynamics of system (2.1).We further compare the dynamics of the continuous model (2.1) and discrete model (2.4) by choosing the bifurcation parameter b, the producer's intrinsic growth rate.From the bifurcation diagram in Figure 4, one can easily observe that the dynamics of model (2.4) are very different from those of the corresponding continuous model for b larger than 2. Figure 4 exhibits another interesting limiting dynamics of model (2.4). Figure 4 (a 1 ) shows the familiar periodic doubling route to chaos for the producer's density.However, the densities of consumer and predator crash and quickly reach zero shortly after the producer enters the chaos zone (Figure 4 (a i ), i = 2, 3).This suggests that chaotic producer densities provide an important factor for the extinction of consumer and predator populations.This interesting and important biological insight is totally lost in the continuous model (2.1) (see Figure 4 (b i ), i = 1, 2, 3).
In summary, our analytical and numerical analysis suggest the following implications: • For the most part, the continuous and the discrete stoichiometric models exhibit the same phenomena.Our discrete model retains most of the important dynamical features exhibited by its continuous analogs.Our analysis shows that qualitatively, the stoichiometric effects of food quality of producer are robust to discretization of time for particular parameter sets.• Our discrete stoichiometric tritrophic model exhibits chaotic dynamics for different parameter ranges than the corresponding continuous model.Discretization enlarges the parameter ranges for the existence of chaotic dynamics.• The discrete model captures some interesting results that the continuous model can not.In the discrete model, when the producer's intrinsic growth rate (b) is large, the consumer and predator can not cope with the sharp fluctuations of the producer's density and tend to perish (Figure 4).However, the continuous model predicts simplier dynamics and fails to capture this phenomenon.While stoichiometric effects of low food quality are presented in both continuous and discrete models, there are quantitative and qualitative effects that arise from discretization.The discretization of time illuminates the mechanism of difference between the continuous model and its discrete analog.The time scaling in ecology is very important and has crucial implications.For practical problem, it is critical that the time scale should be appropriately selected in ecology and the corresponding model can be more effective.This study can possibly serve as such an example of important insights.

Figure 4 .
Figure 4.The bifurcation curves versus b, (a i ) and (b i ) for the discrete model and continuous model, respectively, where i = 1, 2, 3.The parameters are defined in Table1and the initial conditions are x = 0.5 mgCL −1 , y = 0.5 mgCL −1 and z = 0.5 mgCL −1 .

c
2018 the Author(s), licensee AIMS Press.This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)

Table 1 .
Parameters of model (2.4) with default values K Producer carrying capacity limited by light 0 − 10 mgCL −1 Theorem 4.1.For equation (2.4), the origin E 0 is always unstable.The boundary equilibrium E 1 is locally asymptotically stable if Figure 3. Spectrum of the maximum Lyapunov exponent (MLE) versus K (a) for the discrete model and (b) for continuous model.The parameters are defined in Table