Effect of excess food nutrient on producer-grazer model under stoichiometric and toxicological constraints

Accurately assessing the risks of contaminants requires more than an understanding of the effects of contaminants on individual organism, but requires further understanding of complex ecological interactions, elemental cycling, and the interactive effects of natural stressors, such as resource limitations, and contaminant stressors. There is increasing evidence that organisms experience interactive effects of contaminant stressors and food conditions, such as resource stoichiometry, availability and excess of nutrient. Here, we develop and analyze the first producer-grazer population model that incorporates the effects of excess nutrients, as well as nutrient limitations on grazer exposed to toxicants. We use analytical, numerical and bifurcation analysis to reduce and explore model parameterized for an aquatic system of algae and zooplankton exposed to methylmercury under varying phosphorus conditions. Under certain environmental conditions, our models predict higher toxicity than previous models that neglect the consequences excess nutrient conditions can have on grazer populations.


Introduction
Elemental imbalances between producer and grazers can influence life-history traits such as growth, survival, and reproduction.Chemical contaminants can have similar impacts on these traits [5].Recently, there have been experimental efforts to investigate any interactive effects of these stressors.Hansen et al. 2008 [8] found empirical evidence that Daphnia experience strong interactive effects of food quality and the contaminant fluoxetine on their survival, growth, and reproduction.Here they observed that Daphnia fed nutrient-rich algal food experience greater toxicity to fluoxetine.Ieromina et al. 2014 [11] observed that Daphnia fed nutrient-deficient algal food experience high toxicity to the contaminant imidacloprid.They also observed evidence that nutrient excess food can impact the reproductive output for D. magna.While Lessard et al. 2012 [13] observed no interactive effects between food nutrient content and an herbicide contaminant on Daphnia survival but observed great effects on Daphnia growth rates under nutrient-rich conditions.
The toxic contaminant mercury (Hg) can bioaccumulate in aquatic food chains as methylmercury (MeHg) posing a risk to ecosystems and humans [15].Karimi et al. 2007 [12] shows stoichiometric constraints, such as food quality, can affect the accumulation of MeHg in Daphnia.The interactive effects of nutrient availability and aqueous Hg concentration may play a significant role in the bioaccumulation of MeHg.Recently, Peace et al. 2016 [20] developed a stoichiometric producer-grazer model subject to concurrent stoichiometric constraints and contaminant stressors.They parameterized the model for algae and Daphnia exposed to MeHg under variable phosphorus conditions.This model assumes the organisms pay an energy or carbon cost when exposed to contaminants, that reduce the growth rates based on a linear dose response of growth to concentrations of the toxicant.
Following the stoichiometric assumptions of the classical LKE model [14], the stoichiometrictoxicant model developed by Peace et al. 2016 [20] allows the producer Phosphorus:Carbon(P:C) ratio to vary while keeping the grazer P:C ratio constant.This assumption is based on the fact that, although grazer stoichiometries are variable, the range of variation is small compared to the range of producer stoichiometries [22].This often results in elemental imbalances between grazers and their resources.The model assumes the producer is optimal food for the grazer if its P:C ratio is equal to or greater than the P:C of the grazer, thus incorporating the effects of low-nutrient food content on grazer dynamics.
However, it has been reported that grazer dynamics are also affected by excess food nutrient content [3,6].This phenomenon, called the stoichiometric knife edge, reflects a reduction in animal growth not only by food with low Phosphorus(P) content but also by food with excessively high P content.Peace et al. 2013 [18] extended the LKE model [14] in order to incorporate the effects of excess nutrient food content on grazer's growth in the absence of toxicants.They did this by stoichiometrically modifying the grazer's ingestion rate and conversion efficiency.The model assumes the mechanism behind this stoichiometric knife edge phenomenon is the grazer's feeding behavior.It assumes that high P content of food causes the animal to strongly decrease their ingestion rate, perhaps leading to insufficient Carbon(C) intake and thus decreased growth rate.This extended model captures the consequences of imbalances from both low P:C and high P:C food resources in the absence of toxicants.Peace et al. 2016 [20] investigate grazer growth dynamics exposed to a contaminant stresses while under nutrient limited conditions.However, this model does not incorporate the consequences of nutrient excess conditions.Here, we integrate the stoichiometric tox-mediated producer-grazer model [20] with the non-toxicant stoichiometric knife edge model [18].This result is the first model that incorporates the impact of excess nutrients on the producer-grazer dynamics when subject to a toxicant.We parameterize the model for algae and Daphnina exposed to the toxicant MeHg and investigate the impact of varying P and light conditions.

Model Formulation
We start with the Stoichiometric Toxicant-mediated predator-prey model developed by Peace et al 2016 [20] based on Huang et al. 2014 [10]: where variables x and y are the prey and predator population densities (mg C/L) respectively and u and v give the toxicant body burden, or the concentration of the toxicant in the prey and predator, respectively.The function f (x) is the predator ingestion rate.By following the recommendation of the committee on toxicology of the National Research Council in 1992 and tested in Miller et al. 2000 [16], we use the power law to represent the relationship between toxicant concentrations and predator mortality rate.Predator mortality rate as a function of v, takes the form d 2 (v) = h 2 v I + m 2 where h 2 and I are positive constants for coefficient and exponent of the power function and m 2 is the natural loss rate, including both natural mortality and respiration.Parameter α 1 is the prey maximal growth rate and α 2 is the toxicant effect on prey reproduction.Parameters a 1 and a 2 are toxicant uptake rates and σ 1 and σ 2 are toxicant efflux rates for the prey and predator, respectively.Parameter β 1 where 0 < β 1 < 1 is the growth efficiency of the prey, β 2 is the toxicant effect on predator reproduction and ξ is the predator toxicant assimilation efficiency.T is the total toxicant in the system, P is the total amount of phosphorus in the system and θ is the constant predator P : C ratio.The growth rates of both organisms are influenced by the respective body burdens.Toxicological constraints appear as maximum functions in the above Model (2.1).The maximum term 0 ≤ max(0, 1− α 2 u) ≤ 1 comes from a linear dose-response for the gain rate of prey.Similarly, the term max(0, 1−β 2 v) represents a linear dose-response for the predator growth efficiency.Population growth dynamics are also influence by nutrient availability.Stoichiometric constraints appear as minimum functions in the above Model (2.1).The minimum function min K, P−θy q is used to describe the prey carrying capacity.The first input, K is the prey carrying capacity in-terms of carbon or light availability.The second input P−θy q is the carrying capacity determined by phosphorus availability.The numerator here, P − θy, is the total amount of P available for the prey.This is based on the assumption that all P is either in the prey or the predator.The model assumes prey are extremely efficient at taking up nutrients and does not allow free nutrients to be dissolved in the environment.The consumer growth rate is is described by another minimum function, min β 1 , Q θ .The first input, β 1 is used when consumer growth is limited by carbon, the second input, Q θ is used when consumer growth is limited by phosphorus.In, Peace et al. 2016 [20] they incorporate stoichiometric nutrient limitation to the toxicant mediated predator prey model developed by Huang et al. 2014 [10].The current models do not consider the effects of excess nutrient conditions.Here, we develop and analyze the first model that incorporates the impact of excess nutrients on predator-prey dynamics when subject to a toxicant.
The grazer's ingestion rate, f (x) is a monotonic increasing and differentiable function, We assume the ingestion rate is a Holling type II function of the form f (x) = cx a+x where, c is the maximal grazer's ingestion rate and a is the half saturation constant.By following Peace et al. 2013 [18] we assume that the grazer ingests P up to the rate required for its maximal growth but not more.This assumption is based on empirical data from Plath and Boersma 2001[21] which leads to a new expression for the grazer ingestion rate.The ingestion rate, f (x) is replaced by The grazer's growth rate, g(x, y) = min β 1 , Q θ u(x, y), can be written as This expression is also incorporated into the body burden equation and our system takes the following form:  We parameterize our toxicant-mediated stoichiometric knife edged producer-grazer Model (2.3) to a system of algal (producer) and Daphinia (grazer), in order to investigate the effects of co-occurring phosphorus constraints and MeHg stressors on population dynamics and MeHg bioaccumulation.All parameter values were used in Peace et al. 2016 [20] and are listed in Table 1.For the following analysis we assume I = 1 for mathematical convenience.We non-dimensionalize system (2.3) to facilitate the mathematical analysis and rescale the model with the following non-dimensional parameters: Dropping the tildes for simplicity, the dimensionless form of system (2.3) yields ) According to the parameterization in Table 1, = 0.0576.Since is sufficiently small the dynamics of the body burdens u and v are on a much faster time scale than the population dynamics of x and y.Therefore we let → 0 and apply a quasi-steady-state reduction in order to reduce the model to the slow subsystem.Applying this quasi-steady-state reduction and letting → 0 yields Now, substituting (2.6) and Q = P−θy x into (2.5a) and (2.5b) gives us the following quasi-steady-state non-dimensional system

Model Analysis
Note that if T > σ 2 1 then the toxicant levels are too high for the prey to reproduce and grow and (x, y) → (0, 0).For the following analysis we assume that T < σ 2  1 and let k = min K, P q .
Theorem 3.1.Solution of the system (2.7) with the initial conditions in the set will remain there for all forward time.
Proof.First we show the solutions remain in the rectangle R defined by 0, k × 0, P θ , then we show solutions are also bounded by the inequality qx + θy < P. Let S (t) = (x(t), y(t)) be a solution of System (2.7) with S (0) ∈ R. Assume there exists a time t 1 > 0 such that S (t 1 ) touches or crosses a boundary of R for the first time.The following cases prove solutions remain in R by contradiction.Case 1 left boundary: This implies that x(t 1 ) ≥ x(0)e αt 1 > 0, where α is a constant.This contradicts x(t 1 ) = 0 and proves that S (t 1 ) does not reach this boundary.Case 2 right boundary: x(t 1 ) = k = min K, P q Then for every t ∈ [0, t 1 ], Then x(t) < k by a standard comparison argument.This contradicts x(t 1 ) = k and proves that the trajectory does not cross this boundary.Case 3 bottom boundary : y(t 1 ) = 0 Then for every t ∈ [0, t 1 ],

≡ αy
This implies that y(t 1 ) ≥ y(0)e αt 1 > 0, where α is a constant.This contradicts y(t 1 ) = 0 and proves that S (t 1 ) does not reach this boundary.Case 4 top boundary: y(t 1 ) = P θ Then for every t ∈ [0, t 1 ], Then y(t) < P θ by a standard comparison argument.This contradicts y(t 1 ) = P θ and proves that the trajectory does not cross this boundary.The above cases prove the trajectories are bounded in R. Now assume S (0) ∈ Ω and there exists a time t 1 > 0 such that S (t 1 ) touches or crosses a boundary of Ω for the first time.The final case proves solutions remain in Ω by contradiction.Case 5: qx(t 1 ) + θy(t 1 ) = P Then qx(t) + θy(t) < P for every t ∈ [0, t 1 ) and qx (t) + θy (t) ≥ 0 x(t)y(t) x(t)y(t) Then x(t)y(t) = 0 a contradiction.
To investigate the equilibria we first rewrite system (2.7) in the following form where, The Jacobian takes the following forms

Boundary Equilibria
We investigate the equilibria E 0 = (0, 0) The local stability of E 0 = (0, 0) is determined by the Jacobian in the following form, Jacobian, J(E 0 ) has the eigenvalues with opposite signs, thus E 0 is saddle point.In the absence of grazer, the carrying capacity of the producer depends only on the light and phosphorus availability, which we denote as ) is determined by the Jacobian in the following form, The stability of E 1 (k, 0) depends on the sign of G(k, 0).E 1 is locally asymptotically stable if G(k, 0) < 0 and E 1 is saddle point if G(k, 0) > 0.

Interior Euilibria
In this subsection, we analyze the existence and stability of the interior equilibrium numerically with phase plane analysis, Figure 1.The solution of the system are bounded by the trapezoidal region shown below.The phase plane is divided into three biologically significant regions by the two lines This represent the case where P is neither limiting nor in excess.Region II is defined by β 1 > Q θ ; here, growth is limited by the deficiency of P. Region III is defined by , where P is in excess and reduces grazer growth rate.Note that these same regions are described in Peace et al. 2013 [18] in the absence of toxicants.

Numerical Analysis
Peace et al. 2016 [20] investigate the effect of toxicant under stoichiometric constraints on the population dynamics for different light level.In this section we are comparing our result with [20].We use the non-dimensionalized system of both models for simulations.The Black color will represent the results of Peace et al. 2016 [20] and the red color will represent our proposed model.Some instances the models agree and the red line is pllotted on top of the black line.These simulations use the parameter listed in Table 1 with K=1.5mg C/L for varying Phosphorus level.All simulations used the Holling type II function f (x) = cx a+x for the ingestion rate.

Numerical Simulations
Numerical simulations of the reduced model, System (2.7) are presented for the producer in Figure 2 and the grazer in Figure 3 for varying values of P. The red simulations presented in Figures. 2 and 3 correspond to the phase planes presented in Figure 1.
Figure 1.Phase planes for the reduced system (2.7) using parameters found in Table 1 and varying values for P: (a) low total phosphorus P=0.035 mg P/L, (b) P=0.04 mg P/L, (c) P=0.07 mg P/L, (d) excess phosphorus P=0.1 mg P/L.Dashed curves are producer nullclines and solid curves are grazer nullclines.Open circles denote unstable equilibria and filled circles denote stable equilibria.The dotted curve in (b) shows a stable periodic orbit.
For low P values we got the similar results as Peace et al. 2016 [20], see Figures.2(a) and 3(a).As P increases, we start to see the differences between these two models.Figures 2(b) and 3(b) show the bistability for our proposed model, with a stable equilibria and a stable limit cycle.Under these same conditions, the populations oscillate for the Peace et al. 2016 [20] model.As we continue to increase P, the Peace et al. 2016 [20] model predicts similar oscillatory dynamics for moderate and high P values, however in the proposed model System (2.7), the oscillations collapse for higher P values, Figures 2(c,d) and 3(c,d).Very large values of P causes the extinction of grazer population, Figure 3(d).
Figure 2. Numerical simulations for producer densities using parameter found in Table 1 and varying values for P: (a) low total phosphorus P=0.035 mg P/L, (b) P=0.04 mg P/L, (c) P=0.07 mg P/L, (d) excess phosphorus P=0.1 mg P/L.Black curves are solutions of System (2.1) Peace et al. 2016 [20] and red curves are solutions of our reduced model System (2.7).

Bifurcation Analysis
Here we perform a bifurcation analysis of the full model, System (2.3) using K mg C/L as the bifurcation parameter for both low and high P mg P/L, see Figure 4.
Qualitatively, the bifurcation diagrams for low P (Figures 4ab) and high P (Figures 4cd) are similar for both models.For low light levels, the grazer is at a stable equilibrium of a low density.As K increases this interior equilibrium loses stability and stable limit cycles emerge at a Hopf bifurcation.The limit cycles start off with small amplitudes, but as K continues to increase the amplitudes get bigger.Eventually for higher levels of K these limit cycles collapse and another interior equilibrium point gains stability.This occurs at a saddle node bifurcation.As K increases further to very high values the grazer density begins to decrease.Here, the grazer is limited by food quality, or P-limitation.Figure 3. Numerical simulations for grazer densities using parameter found in Table 1 and varying values for P: (a) low total phosphorus P=0.035 mg P/L, (b) P=0.04 mg P/L, (c) P=0.07 mg P/L, (d) excess phosphorus P=0.1 mg P/L.Black curves are solutions of System (2.1) Peace et al. 2016 [20] red curves are solutions of our reduced model System (2.7).
Under low P conditions we see differences between the two models for low K conditions (Figures 4ab).First, under very low K, there is region where the previous model predicts a stable grazer population at a low density, whereas our model predicts grazer extinction (Figure 4a).Under low light conditions Q, the algae P:C ratio, is large.In both models, there is a low quantity of food available for the grazers.In the previous model, the grazer is limited by food quantity.However, our model incorporates consequences of high P content, and the grazer is limited by food quality, This is the influence of the stoichiometric knife edge.It is important to note the differences in the predicted grazer body burdens at low K levels (Figure 4b).Here, there is range of low K values, where the previous models predicts lower body burdens than our model.We also observe a shift in the location of the Hopf bifurcation and the amplitude of the cycles.The previous model exhibits limit cycles at a lower value of light than our model, and the cycles have larger amplitude.For high light levels, the two models predict qualitatively similar dynamics.Indeed, here the low P and high K conditions yield low Q values.Here the grazer is limited by P in both models.
The differences between the predictions of the model are more significant under high P levels, where effects of the stoichiometric knife edge are more likely to be occur (Figures 4cd).Here, the entire bifurcation diagram is shift significantly to the right for our model, note the bifurcation diagram takes K from 0 to 5 mg C/L, in these figures.Light levels above 5 mg C/L are unrealistic in natural and laboratory settings.The Hopf bifurcation occurs at much higher light levels.Consequently, the dynamics of the two models are qualitatively different even for high light levels.

Conclusion
Elemental mismatches between trophic levels often leads to consequences on population growth.Food resources are rarely optimally suited for grazers in terms of their nutrient content.Daphnia grazers are often faced with P limitation, when their diet consists of algae with low P:C [22].On the other hand, when their diet consists of algae with high P:C they have to deal with this P excess and can become limited by C [21,18,19].These types of nutritional constraints influence life-history traits, such as growth rates and conversion efficiencies.Toxicants can have similar impacts on these life-history traits.We developed and analyzed the first model that incorporates the impact of both low and excess nutrients on the producer-grazer dynamics when subject to a toxicant.Previous models that consider concurrent nutrient and toxicological constraints consider only nutrient limited conditions [20].Here, we expanded the model by Peace et al. 2016 [20] to incorporate the effects of nutrient excess conditions on population dynamics, a phenomenon called the stoichiometric knife edge.Peace et al. 2013 [18] modeled this phenomenon without the presence of toxicants.The proposed model (2.3) yields additional insights for risk assessment compared to previous work, especially under excess food nutrient conditions.
To facilitate the analysis of the model we assumed population metabolism occurs on a faster time scale than population growth dynamics.Considering the fast and slow systems in the model and a quasi-steady-state assumption we were able to reduce system (2.3) down to the two dimensional system (2.7).We showed the boundedness of the reduced system (2.7) analytically.We investigated the existence and local stability of boundary equilibria analytically.We numerically observed the existence and stability of interior equilibria and limit cycles.
We compared the proposed model (2.7) results with those of system (2.1) developed by Peace et al. 2016 [20].For low values of P the population dynamics are similar (Figure2a,Figure3a), however  [20] and red corresponds to our full model, System (2.3).Solid curves are stable equilibria and the maximum and minimum of stable limit cycles.Dashed curves are unstable branches.
they begin to deviate for moderate and high P loads (Figure2b-d,Figure3b-d).The Peace et al. 2016 [20] model predicts periodic cycles for moderate and high P loads, however our models captures the effects of excess P content on the populations and predicts extinction of the grazer population for high P value.The most important differences between the two model predictions can be seen in the bifurcation diagrams (Figure 4).In particular see Figure 4(b,d), under low light levels K near 0.5 mg C/L our model predicts a high body burden for the grazer.Under these conditions previous models may be under predicting the effects of toxicants which can have important implications on risk assessment protocols developed using these models.
In order to incorporate the stoichiometric knife edge phenomenon we following the assumption of Peace et al. 2013 [18] and considered a reduction in grazer ingestion rate under high P food conditions.However, the mechanisms underlying this phenomenon is likely more complicated than a simple re-duction in foraging behaviors [7].For example nutrient and carbon absorption rates of ingested food, gut passage time, and other metabolic processes such as respiration rates may play a role.It is possible that these processes are also influences in the presence of a toxicant.Additional experiments on metabolic rates under multiple constraints can help us modify the models presented here to consider these extensions.
While these models are parameterized for a system of Daphnia and algal under varying P and MeHg conditions, the models have broader applicability.Reduction in growth rates when exposed to high nutrient conditions has been observed for a variety of species [21,4,17].At the same time many system experience exposure to a multitude of anthropogenic and natural toxicants [25].Reparameterization of these models can yield additional insight into other ecological producer-grazer systems subject to concurrent stoichiometric and toxicological constraints.

Figure 4 .
Figure 4. Bifurcation diagrams of System (2.3) for the grazer biomass (y mg C/L) and body burden (v µg MeHg/mg C) using K as bifurcation parameter between for (a,b) low P = 0.03 mg P/L and (c,d) high P = 0.07 mg P/L.Black corresponds to model (2.1), developed by Peace et al. 2016[20] and red corresponds to our full model, System (2.3).Solid curves are stable equilibria and the maximum and minimum of stable limit cycles.Dashed curves are unstable branches.