Mechanistically derived Toxicant-mediated predator-prey model under Stoichiometric constraints

1.   Introduction

While the development of models in Ecotoxicology has advanced understandings of how contaminants cycle through food webs and impact organisms ^[1,2,3,4], many existing models do not consider the important role that environmental nutrients can have. The interactive roles that essential elements, such as nitrogen and phosphorus, can play with environmental toxicants, such as mercury and cadmium, on population dynamics and toxicities can have significant implications ^[5]. Recent research on aquatic food webs has observed that the impact of toxicants on the population dynamics of Daphnia grazers depends on food quality ^[6,7,8]. These examples make it clear that the strength and type of any interactive effects of contaminants and nutrient stressors depends on the elements in the system. Considering a producer-grazer system of algae and Daphnia exposed to methyl mercury (MeHg) ^[8] showed that bioaccumulation dynamics depends on phosphorus (P) limitation. The accumulation of MeHg in Daphnia depends on the nutritional value (phosphorus to carbon ratio, P:C) of its algal diet.

Researchers ^[9,10] developed and analyzed models of this same system parameterized with data from ^[8]. These models capture the somatic growth dilution phenomenon, where organisms experience a greater than proportional gain in biomass relative to toxicant concentrations when consuming high vs low quality food, as observed in empirical data ^[8]. These modeling efforts are some the first dynamical models formulated using the theory of Ecological Stoichiometry ^[11] to investigate both stoichiometric and toxicological constraints on population dynamics. While tracking the transfer of both P and MeHg across trophic levels, these models make a simplifying assumption that environmental aqueous MeHg loads are constant and not influenced by the population dynamics of algae or Daphnia. While this simplifying assumption reduced model complexity it may lead to consequences on predicting the important role nutrient cycling can have on toxicant bioaccumulation.

Here we relax the assumption that environmental toxicant loads are constant by explicitly modeling the amount of aqueous toxicant as it depends on the organisms' uptake and efflux rates. We develop and analyze a stoichiometric-toxicant mediated producer-grazer model that explicitly tracks environmental toxicant loads. We explore parameter ranges analytically and numerically to investigate when this simplifying assumption has major implications on population dynamics and toxicity predictions.

2.   Model formulation

We start with the stoichiometric toxicant-mediated knife edged predator–prey model developed by ^[10]:

where $Q = (P-\theta y)/x$ is the prey variable phosphorus to carbon ratio (P:C). Here $x(t)$ and $y(t)$ are the biomass of the prey and predator, respectively, measured in terms of C. $u(t)$ and $v(t)$ give the toxicant body burden, or the concentration of the toxicant in the prey and predator, respectively. $\sigma_1$ and $\sigma_2$ are toxicant efflux rates for the prey and predator, respectively. $T$ is the total toxicant in the system. $\xi$ is the predator toxicant assimilation efficiency. $P$ is the total amount of phosphorus in the system. $q$ is the minimum P:C ratio of the prey. Parameter $\theta$ is the constant predator P:C ratio, $\alpha_1$ is the prey maximal growth rate, $\alpha_2$ is the effect of the toxicant on the prey growth rate, and $a_1$ and $a_2$ are toxicant uptake rates. Parameter $\beta_1$, which is between $0$ and $1$, is the growth efficiency of the prey and $\beta_2$ is the toxicant effect on predator reproduction. The grazer's ingestion rate, $f(x)$ is a monotonic increasing and differentiable function, $f'(x) > 0$ and $f''(x) < 0$ with $\lim\limits_{x\rightarrow\infty}f(x) = c$. Model (2.1) is based on previous modeling works by ^[4,9] and ^[12].

By following the recommendation of the committee on toxicology of the National Research Council in 1992 and tested in ^[9,10,13] uses the power law to represent the relationship between toxicant concentrations and predator mortality rate. Predator mortality rate as a function of $v$, takes the following form

where $h_2$ and $I$ are positive constants for the coefficient and exponent of the power function and $m_2$ is the natural loss rate, including both natural mortality and respiration.

Population growth dynamics are influence by nutrient availability. Stoichiometric constraints appear as minimum functions in the above model (2.1). The minimum function $\min \left\{K, (P-\theta y)/q \right\}$ is used to describe the prey carrying capacity. The first input, $K$, is the prey carrying capacity in terms of C and represents the light intensity. The second input, $(P-\theta y)/q$, is the carrying capacity determined by phosphorus availability where $P-\theta y$ is the prey's phosphorus content. This is based on the assumption that all available nutrients are 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 described by another minimum function $\min\{\beta_1, Q/\theta\}$. The first input, $\beta_1$, is used when consumer growth is limited by carbon. The second input, $Q/\theta$ is used when consumer growth is limited by phosphorus.

The growth rates of both organisms are also influenced by the respective body burdens. Toxicological constraints appear as maximum functions in the above model (2.1). The term $\max(0, 1-\alpha_2u)$, which is a fraction between $0$ and $1$, represents a linear dose response for the prey's gain rate. The term $\max(0, 1-\beta_2v)$ represents linear dose response for the predator's reproduction efficiency.

Researchers of ^[10] considered the total toxicant $(T)$ as constant over the time, they didn't distinguish between the toxicant in the environment and total toxicant in the system, both are considered as $T$. In our model we consider the environmental toxicant as $T_e$ and the total toxicant in the system as $T$ where $T = ux+vy+T_e$.

We extend system (2.1) by allowing the environmental toxicant in the system to change over time rather than being constant, as it will give more insight about the population dynamics and risk assessment. For our system we still assume that the total toxicant is constant, in other words the system is closed for any input or output of toxicant. Toxicant is taken up by the prey and predator at rates of $a_1 T_e$ and $a_2 T_e$ respectively. Therefor aqueous toxicant is removed from the environment due at rates $a_1T_ex$ and $a_2T_ey$. The prey and predator populations release toxicant back into the environment at rates $\sigma_1ux$ and $\sigma_2vy$ respectively. Additionally, we assume that predator death, $d_2(v)$, results in toxicant released back into the environment at rate $d_2(v)vy$. We can express the change of toxicant in the environment over time by a differential equation as follows

The last term in equation (2.2) accounts for the toxicant released back into the environment following the predator's assimilation of ingested prey. The predator's toxicant assimilation coefficient is $\xi$, the predator gains $f(x)u$ of toxicant from predation but only assimilates $\xi f(x)u$ amount of toxicant and fails to assimilate $\left(1-\xi\right)f(x)u$ amount of the toxicant, so $\left(1-\xi\right)f(x)uy$ amount of the toxicant returns back to the environment. Since, we assume the system is closed the total amount of toxicant is constant, $\frac{d}{dt}\left(T\right) = 0$. This is verified below:

Combining the system (2.1) and (2.2) we get

Again, we assume the total toxicant of the system is constant,

and therefore,

which allows us to substitute $T_e = T-ux-vy$ into system (2.4) and reduce the above system to a four dimensional system for convenience. Model (2.5) gives the final form of our mechanistically derived toxicant-mediated stiochiometric predator-prey model and Table 1 lists the parameter values.

3.   Model analysis

To aid in model analysis we first non-dimensionalize and employ a quasi-steady-state assumption to reduce the model to a two-dimensional system. We then investigate boundary equilibria analytically and use phase-plane analysis to investigate interior equilibria.

3.1.   Model reduction

Here, we non-dimensionalize the system (2.5) and assumed $I = 1$ for mathematical convenience to facilitate the mathematical analysis and re-scale the model with the following non-dimensional parameters:

Dropping the tildes for simplicity, the dimensionless form of system (2.5) yields

According to the parameterization in Table 1, $\epsilon = 0.0576$. Since $\epsilon$ 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$. Here we make the assumption that the exchange of the toxicant between the water and the organisms is faster than the other processes including assimilation, dilution, and organismal growth. This assumption is similarly made by previous ecotoxicology models ^[3,4,19]. Considering the fast and slow systems in the model and a quasi-steady-state assumption we let $\epsilon\to0$ in order to reduce the model to the fast-slow subsystem.

Applying this quasi-steady-state reduction and letting $\epsilon\to0$ yields

Now, substituting (3.3) and $Q = \frac{P-\theta y}{x}$ into (3.2a), (3.2b) gives us the following quasi-steady-state non-dimensional system

Considering the fast-slow subsystems and using the quasi-steady-state assumptions essentially decouples the dynamics of the trophic transfer of the toxicant from population dynamics and has been used in previous work to write an expression for the bioconcentration factor ^[19].

3.2.   Equilibria analysis

To investigate the equilibria we first rewrite system (3.4) in the following form

where,

The Jacobian takes the following forms

3.2.1.   Boundary equilibria

We investigate the equilibria

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

So, $x^* = {\bf{k}}$. The local stability of $E_1 = ({\bf{k}}, 0)$ is determined by the Jacobian in the following form,

The stability of $E_1({\bf{k}}, 0)$ depends on the sign of $G({\bf{k}}, 0)$. $E_1$ is locally asymptotically stable if $G({\bf{k}}, 0) < 0$ and $E_1$ is saddle point if $G({\bf{k}}, 0) > 0$.

3.2.2.   Interior equilibria

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 $\beta_1 = \frac{Q}{\theta}$ and $f(x) = \frac{c\theta}{Q}$. Region Ⅰ is defined by $\beta_1 < \frac{Q}{\theta}$ and $f(x) < \frac{c\theta}{Q}$. This represents the case where P is neither limiting nor in excess. Region Ⅱ is defined by $\beta_1 > \frac{Q}{\theta}$; here, growth is limited by the deficiency of P. Region Ⅲ is defined by $\beta_1 < \frac{Q}{\theta}$ and $f(x) > \frac{c\theta}{Q}$, where P is in excess and reduces grazer growth rate. Note that these same regions are described in ^[10] without tracking environmental toxicants and in ^[12] in the absence of toxicants.

It is important to note the the model is nonsmooth and there are minimum operators in the elements of the Jacobian matrices, as they are derivatives of non-differentiable functions. This leads to the observed phase plane fragmentation and partitioning of the parameter space (Figure 1). These dynamics are common in nonsmooth stoichiometric models and have been thoroughly analyzed in previous models without toxicants. Robust analyses and global bifurcation has been conducted on the classical LKE model ^[20] with varying Holling type functional responses ^[21,22,23] as well as the modified LKE model incorporating the stoichiometric knife-edge ^[24] with Holling type Ⅱ functional response ^[25]. The above use of the Jacobian shows only local stability of boundary equilibria, which was verified via numerical simulations. More rigorous analyses confirming the Jacobian is well defined and global bifurcation may shed more light into these dynamics and is left for future work.

4.   Numerical analysis

In this section we perform numerical simulation for our reduced model (3.4) and compare the results with a previous model by ^[10], which neglects to track environmental toxicant loads. Parameters used are given in Table 1 with varying phosphorous conditions. In the following subsection we present bifurcation diagrams of the reduced model (3.4) considering $P$ and $T$ as bifurcation parameters.

4.1.   Numerical simulations

Numerical simulations of the reduced model, System (3.4) are presented for the producer in Figure 2 and the grazer in Figure 3. We use the parameters given in Table 1 with $P$ varied form $P = 0.01-0.15$ mg P/L, $K = 1.5$ mg C/L and $T = 0.01\mu$g. The solid curves presented in Figures 2 and 3 are solutions that correspond to the phase planes presented in Figure 1. The dashed curves are the solutions from reduced model of ^[10].

We observed that for low $P$ value ($P = 0.03$) and moderate $P$ values ($P = 0.05, 0.1$) our model provides different results than ^[10] (see Figure 2a, b.c and Figure 3a, b.c) but for high $P$ value ($P = 0.15$) both models give us the same results for the producer population dynamics (see Figure 2d and Figure 3d). Also, our model predicts oscillations when $P = 0.05$ but we don't see any oscillation for ^[10] (see Figure 2c and Figure 3c). In addition to examining equilibria predictions, it is worth noting the interesting transient behaviors exhibited under some parameter sets. Certain indices can be used to measure the extent and duration of transient growth in cases where the solutions approach asymptotically stable equilibria ^[26]. In particular we observe solutions with high reactivity, which is a measurement of the maximum possible growth rate captured during transient dynamics ^[26], see for example the solid curves in Figures 2c, 3c and the dashes curve in Figures 2b, 3b.

4.2.   Bifurcation analysis

In this subsection we perform bifurcation analysis to our System (3.4) by using total phosphorus ($P$) and total toxicant ($T$) as a bifurcation parameter separately. Red curves corresponds to our proposed model and Black curves corresponds to the previous model ^[10]. Solid curves are stable equilibria and the maximum and minimum of stable limit cycles. Dashed curves are unstable branches.

We observed there are similarities to the bifurcation structure between the two models, however important differences also exist. For low values of P, the solution of both models exhibit a stable equilibrium with a high producer ($x$) density and low grazer ($y$) density, Figure 4. Here, the grazer density is low due to low food quality. As P increases, this interior equilibrium looses its stability and stable limit cycles emerge. The limit cycles emerge at a lower P value in our System (3.4) than they do for the previous model. Additionally, the amplitudes of the limit cycles in our System (3.4) are much larger. As P continues to increase, the limit cycles loose stability and another interior equilibrium gains stability at a Hopf bifurcation. This bifurcation occurs at a higher value of P in our model then the previous model. For very large values of P, eventually the grazer density begins to decline in both models, however this decline occurs sooner and more rapidly decline in the previous model. For extremely large P values the grazer is driven to extinction. The bifurcation structure for P are similar to previous works that incorporate the stoichiometric knife edge phenomenon into producer-grazer models ^[10,12,24], and we refer the reader to those works for additional details on the bifurcation.

Figure 5 shows the bifurcation analysis using $T$ as the bifurcation parameter. While the bifurcation dynamics are qualitatively similar to those of ^[10] the bifurcation points on the proposed model are shifted to the right. When the contamination increases from left to right, according to the previous model predator goes extinct quickly after the contamination level reached around $T = 0.01\mu$g T/L, whereas the grazer continues to persist as toxicant reach high values ($T = 0.07 \mu$g T/L) for the our model.

5.   Conclusion

Dynamical mathematical models have made significant contributions to ecotoxicology and predicting risk due to bioaccumulation of natural and anthropogenic toxicants. Toxicants and essential elements interactively influence population dynamics and the trophic transfer of each other. Predicting accurate body burdens is a major step is toxicological risk assessment protocols. The mechanistically derived two-trophic food web model we developed predicts different population dynamics than previous models, which may have significant impacts on accurate predictions.

While biologically realistic parameter values are often dynamic, modeling approaches often need to consider many parameters as constants in order to reduce model complexity and remain trackable. Initial Ecological Stoichiometric population models assumed that the total nutrient levels in the environment were constants ^[12,20]. Over the last decades, models dropped this simplifying assumption and incorporated dynamic environmental nutrient levels ^[1,24]. These provided much better qualitative information about the predator–prey dynamics under certain nutrient conditions. Here, we take the same approach to track environmental toxicants explicitly.

We relaxed the assumption of environmental toxicant being constant, and we were still able to perform similar mathematical analyses of our model. Furthermore, our model extension did not cause unnecessary model complexity or add any new parameters. We developed a four-dimensional knife–edged stoichiometric predator–prey model (2.5) which explicitly tracks the toxicant load of the populations and as well as in the environment. We used the slow-fast subsystem and quasi-steady state reduction to reduce the system (2.5) to a two-dimensional system (3.4) for the mathematical analysis. We investigated stability of the boundary equilibria analytically and interior equilibria with phase plane analysis (see Figure 1). Numerical simulation of the population dynamics by varying nutrient conditions are presented and compared to the previous model ^[10] (see Figures 2, 3).

Bifurcation diagrams of the populations are investigated and compared to the previous model ^[10], wihch deepen understanding and may improve prediction of risk assessments (see Figures 4, 5). Here, the presented analyses were performed on the reduced model (3.4) rather than the full model (2.5). The results of the numerical simulations and bifurcation analyses incorporate the quasi steady-state assumption and reduction to the fast and slow subsystems. We note that the full model may produce modified predictions. One major consequence that previous models neglect to capture are the wide amplitude limit cycles observed in our model. This is especially true for producer densities under Phosphorus enrichment, see Figure 4a, which exhibit values very close to zero where they are endanger of stochastic extinction. Although this modeling approach is parameterized for an algal-Daphnia system under the toxicant effect of Methylmercury, this modelling framework can be used in other predator–prey systems and as well to predict the risk of bioaccumulation of other toxicants. Also, this works can be extended to more complicated food web models with multiple species and higher trophic levels.

Although the model assumes all the parameters are constant in time it is well known that ecological parameters often exhibit seasonality ^[27,28]. In this study we did not consider the seasonal variation of the environments but focused efforts on relaxing the assumption that the environmental toxicant is not constant and varies due to population dynamics. In this work we explicitly track environmental toxicant loads. We have done recent work on the influence of seasonal producer carrying capacities on a similar system that assumes environmental toxicant loads are constant ^[29]. Indeed seasonality can play a major role on model predictions. Future models can consider seasonal effects and dynamics environmental toxicant loads. Seasonality has been shown to have significant impacts on transient dynamics ^[30], where the reactivity of systems can be magnified when parameters vary over time.

Acknowledgments

This work was supported in part by NSF under Grant No 1615697. Authors thank Dr.Linda Allen for useful comments on developing and reducing the model.

Conflict of interest

All authors declare no conflicts of interest in this paper.