Dynamics of a stoichiometric producer-grazer system with seasonal effects on light level

1.   Introduction

All living things are made of chemical elements such as hydrogen (H), oxygen (O), nitrogen (N), phosphorus (P), sulfur (S), and carbon (C), along with energy from light. The study of the balance of energy and multiple chemical elements (mainly C, P, and N) in ecological interactions to evaluate the importance of physiological constraints on organisms and, in turn, the impact of these constraints on food webs and ecosystem functioning is called ecological stoichiometry ^[1]. Ecological stoichiometry has been built to study nutrient-limited ecosystems and offer a theoretical framework to understand how food quality (generally considered through carbon/nutrient ratios) influences organisms (e.g., growth rates, reproductive success) and their roles in ecosystems (e.g., biomass production, consumer-driven nutrient recycling, nutrient and energy fluxes in food webs) ^[2].

The rate at which photosynthetic organisms produce organic compounds in aquatic ecosystems depends basically both upon light availability and upon the supplies of nutrients such as N and P ^[3]. These resources affect the total density of primary producers which then support growth and reproduction of higher trophic levels. Nonetheless, the ratio of C to P may also shape nutritional quality of the primary producer ^[4].

Loladze et al. 2000 ^[4] constructed and studied a stoichiometric Lotka-Volterra type model of the first two trophic levels of an aquatic food chain (algae-daphnia), known as the LKE model. The analysis shows that indirected competition between predator and prey for the limiting nutrient, P, can shift predator-prey interactions from a $(+, -)$ type, where the predator has a positive effect on the growth rate of its prey, and the prey has a negative effect on its predator to an exploitative competition $(-, -)$. This leads to dynamics with two equilibria, where bistability and deterministic extinction of the predator are possible. The most significant dynamics is the delivery of bistability as a result of large values of the variable-carrying capacity determined by whatever is limiting, either by C or by P. On the other hand, most organisms experience environmental seasonality, which often induces periodic fluctuations in the availability of resources (i.e., food, nutrients, or energy), and it is a well-known fact that periodic fluctuations can also cause variations in the producer-carrying capacity ^[5]. Indeed population density tends to fluctuate as a result of seasonal changes in environmental factors ^[6]. Therefore, studying the effects of seasonality in the carrying capacity will contribute to our understanding of predator-prey relationships and their dynamical evolution.

In this paper, we formulate the stoichiometric prey-predator model more mechanistically by modeling the seasonal variation of the producer's light-dependent carrying capacity. Our model is a system of two non-autonomous and non-smooth ordinary differential equations. We have done bifurcation analysis of the model with respect to average value of light-dependent carrying capacity and we have studied the effects of different amplitudes of seasonal forcing on the dynamics of the model. Through the bifurcation analysis, we observe that seasonal forcing can lead to complicated population dynamics, including periodic and quasi-periodic solutions.

The paper is organized as follows. In section2, we describe our model and provide basic analysis of positivity and boundedness. In section3, we conduct numerical simulations and present the results of bifurcation analyses. In section4, we discuss the implications of our results and suggest directions for further research. Proof of mathematical result is placed in the appendix.

2.   Model

Our seasonal producer-grazer system will be based on the following continuous time producer-grazer system developed by Loladze et al. 2000 ^[4]

where $x$ is the density of producer (in milligrams of carbon per liter, mg C/l), $y$ is the density of grazer (mg C/l), b is the intrinsic growth rate of producer (day$^{-1}$), $d$ is the specific loss rate of grazer that includes metabolic losses (respiration) and death (day$^{-1}$), $\hat{e}$ is maximal production efficiency of grazer, $K$ is the light-dependent producer constant carrying capacity in terms of C, $P$ is the total phosphorus in the system, $\theta$ is the fixed grazer P:C ratio, $q$ is the minimum producer P:C ratio and $f(x)$ is the functional response of ingestion to producer abundance which is a Holling type II functional response the following form, $f(x) = cx/(a+x)$ where $c$ is the maximal producer ingestion rate, and $a$ is the half saturation rate of grazer ingestion.

To facilitate our model formulation and its comparison with model (2.1), it is convenient to recall the main model assumptions. They are listed below:

Assumption 1. The total mass of phosphorus in the entire system is fixed, i.e., the system is closed for phosphorus with a total of $P$ (mg P/l).

Assumption 2. Phosphorus to carbon ratio $(P:C)$ in the producer varies, but it never falls below a minimum $q\ (mg\ P/ mg\ C);$ the grazer maintains a constant $P:C, \theta \ (mg\ P/ mg\ C).$

Assumption 3. All phosphorus in the system is divided into two pools: phosphorus in the grazer and phosphorus in the producer.

From these assumptions, new constraints in terms of P and C can be introduced. According to Liebig's law of the minimum, which states that growth is controlled by the strongest limiting environmental factor ^[1], two minimum functions have been obtained. The first minimum function is used to describe the producer's carrying capacity, $\min\{K, (P-\theta y)/q\}$. The first input, $K$, is the carrying capacity determined by light availability. The second input, $(P- \theta y)/q$ is the carrying capacity determined by P availability. It is also noticeable that the carrying capacity of the producer depends on grazer density. The second minimum function is used to describe the grazer growth rate, $\min\{1, (P-\theta y)/x \theta\}$. The first input, $1$, is used when grazer growth is limited by carbon. The second input, $(P-\theta y)/x \theta$, where $(P-\theta y)/x$ is the variable producer P:C ratio, is used when grazer growth is limited by P.

In ecosystems, the light-dependent carrying capacity of primary producers is season-dependent. Therefore, one additional assumption is necessary.

Assumption 4. $K$ is a periodically varying function of time and is considered as

where $K_0$ is the average value of $K$ and $\varepsilon$ is the amplitude of the seasonal forcing; $\varepsilon K_0$ is the magnitude of the perturbation and $\omega$ is the angular frequency of the fluctuations caused by seasonality.

Based on model system (2.1) and assumptions (A1)-(A4), the Seasonal LKE model (called the SLKE) takes the following form:

with period $2\pi/ \omega =$half-year.

In the SLKE model, the carrying capacity of the producer also depends on grazer density. In the absence of the grazer, the carrying capacity of the producer depends only on seasonal light and phosphorus availability, which it is denoted as

Since $K$ cannot be negative, $\varepsilon$ lies in the interval $[0, 1]$. $\varepsilon = 0$ corresponds to absence of seasonality, and $\varepsilon = 1$ means the maximum value of the parameter is twice its average value, i.e., the maximum value of the seasonal light availability, $K_{max}$, corresponds to $2K_0$. Now, by considering the maximum seasonal light availability and the grazer's absence, we can redefine the producer carrying capacity as

Here, we state a lemma and a theorem on the boundedness and positive invariance of the solutions of Model (2.2) respect to $K_{max}$.

Lemma 2.1. Solutions with initial conditions in the open rectangle $\{(x, y):0 < x < k = \min\{K_{max}, P/q\}, 0 < y < P/ \theta \}$ remain there for all forward times.

Proof. The proof is provided in Appendix A.

Theorem 2.1. Solutions with initial conditions in the open trapezoid (or triangle if $K_{max} \geq P/q$)

remain there for all forward times.

Proof. The proof is very similar to that in Loladze et al. 2000 ^[4] (Appendix C), so the details are omitted here.

3.   Numerical Simulations

In this section, we provide numerical simulations for Model (2.2) to support and complement our analytical findings and compare with those of the LKE system. In particular, we will discuss the effect of seasonal variation on the modeling of the dynamical interaction between producer and grazer. We choose the Monod type function $f(x) = cx/(a+x)$ as the functional response of the grazer. We will vary the two parameters in the forcing term, $K_0$ and $\varepsilon$, for fixed values of the other parameter values which are biologically realistic values. These parameter values listed in Table 1 are adapted from Andersen 1997 ^[7] and Urabe and Sterner 1996 ^[8] except $K_0$, $\varepsilon$, and $\omega$.

The parameter $K_0$ varies from 0 to 3 mg C/l, which enables us to compare the solutions of the SLKE model with the solutions of the LKE model. The parameterization of $\omega$ determines the length of the seasonal period or the number of peaks for the producer carrying capacity a year. Cebrián and Valiela 1999 ^[9] investigated seasonal trends in phytoplankton biomass from multiple northern temperate coastal ecosystems via the compilation of multiple published data sets and observed that bimodal cycles, displaying two annual peaks are most commonly observed. Therefore, we assume the angular frequency of the fluctuations is $4\pi/365$, which results in a seasonal light-dependent producer carrying capacity with two peaks a year. The SLKE model system is analyzed with the following initial conditions: $x = 0.5$ (mg C/l), $y = 0.25$ (mg C/l).

The producer and grazer solutions of system (2.2) are plotted against time (days) for various values of light-dependent carrying capacity $K_0:0.25, 0.75, 1, 2$ mg C/l and different forcing amplitudes $\varepsilon:0, 0.1, 0.3, 0.7$. The $\varepsilon$ values chosen represent the cases: (a) for no forcing ($\varepsilon = 0$), which corresponds exactly to the LKE model; (b) for weak ($\varepsilon = 0.1$) seasonal forcing; (c) for intermediate ($\varepsilon = 0.3$) seasonal forcing; and (d) for high ($\varepsilon = 0.7$) seasonal forcing.

In Figure 1, $K_0$ is varied with no forcing ($\varepsilon = 0$). Under low light conditions, the producer is limited by carbon (see Figure 1a). Here the low light levels correspond with a low producer carrying capacity and producer densities are low. In these conditions the producer P:C ratio is high, however producer density is low enough that grazers are maintained at fairly low levels. When light and thus $K_0$ is increased, the predator-prey interaction becomes unstable (see Figure 1b), due to Rosenzweig's paradox of enrichment ^[10]. Further increase in $K_0$ leads back to a stable system with high density in both the producer and grazer, but both of them reach much higher densities than in Figure 1a (see Figure 1c). As $K_0$ is increased still further, grazer density steadily decrease due to low quality food. Finally, under high light condition (see Figure 1d) the producer density is high, but its low quality drives the grazer to extinction.

In Figure 2-4, we give four plots for each of three different values of $\varepsilon:0.1, 0.3$, and $0.7$, respectively.

Low light conditions ($K_0 = 0.25$): The grazer density oscillates periodically above the producer density with weak seasonal forcing (see Figure 2a). When the seasonal forcing is increased, these periodic oscillations continue and the grazer density gradually decreases until it undergoes prolonged near-extinction state and then recovers from it (see Figure 3-4a).

Medium light conditions ($K_0 = 0.75$): At weak seasonal forcing quasi-periodic oscillations (QPOs) can be observed in population densities (see Figure 2b). When the seasonal forcing is increased, the population densities continue to oscillate quasi-periodically (see Figure 3b). With further increased in the seasonal forcing, QPOs transform into periodic oscillations (see Figure 4b).

High light conditions ($K_0 = 1$): The system displays periodic oscillations in population densities for varying seasonal forcing (see Figure 2-4c).

Very high light conditions ($K_0 = 2$): For weak seasonal forcing, the grazer is driven to extinction (see Figure 2d). However, when the seasonal forcing is increased, the grazer density undergoes prolonged near-extinction state and then recovers from it (see Figure 3-4d). It is shown that seasonal forcing plays a very important role in helping the grazer to survive.

3.1.   Bifurcation Analysis

Here, we present bifurcation diagrams, which highlight the extreme grazer population densities observed during oscillatory solutions (see Figure 5), as well Poincaré bifurcation diagrams, which characterizes the oscillatory dynamics as periodic and quasi-periodic (see Figure 7). We consider $K_0$ as the bifurcation parameter that varies from 0 to 3. We will also vary the other parameter in the forcing term, $\varepsilon,$ for fixed values of the other parameter values; following Table 1. In Figure 5, the local minimums and maximums of the grazer population densities are shown in black and red, respectively.

Figure 5a shows bifurcation diagram for $\varepsilon = 0$, i.e. no seasonality, which corresponds exactly to the LKE model (2.1). Under low light conditions, the grazer extincts due to starvation. When light and thus $K_0$ is increased, the grazer density steadily rises and a stable periodic orbit can be observed following a Hopf bifurcation. As $K_0$ continues to increase the periodic orbit collapses and an interior equilibrium gains stability following a saddle-node bifurcation. The grazer population begins to decline due to stoichiometric constraints and eventually low quality food drives the grazer to extinction.

With low level of seasonality ($\varepsilon = 0.1$) the SLKE model (2.2) exhibit similar dynamics, however the equilibrium are now small amplitude periodic orbits (see Figure 5b). As seasonality increases, the region of large amplitude oscillations widens and the amplitudes become more extreme (see Figures 5b-d). Under high level of seasonality the grazer densities approach very low values, where they are in danger of stochastic extinction. The highlight threshold that cause grazer extinction increases as seasonality increases. Despite the large amplitude oscillations, under high light conditions seasonality can rescue grazers from extinction.

Finally, to complement our use of bifurcation diagrams, we present a number of Poincaré diagrams (Figure 6) and Poincaré based bifurcation diagrams (Figure 7) as a means to qualitatively characterize the behavioral differences amongst solutions for a variety of $\varepsilon$ and $K_0$ values. This will allow for a more precise discussion of the rich behavior one can see in Figures 5b-d, which is not present in Figure 5a. Additionally, we use $\varepsilon = 0.5$ to investigate how the dynamics transition from low to high seasonality (see Figure 7c).

The Poincaré diagrams are constructed by using Runge-Kutta-Fehlberg numerical integration with a given set of initial conditions $x = 0.5$ (mg C/l), $y = 0.25$ (mg C/l), and parameters $K_0$ and $\varepsilon$. After waiting for transient effects to die out (after 50 years in this case), $x$ and $y$ pairs are plotted once per period of the forcing function (every half year) in the resulting Poincaré section until the thousands of points necessary to make a clear picture are achieved. We use the Poincaré map to classify system behavior. Distinct points indicate periodicity (pictured in the Poincaré bifurcation diagrams below as single $y$-values for a given $K_0$, or in a typical Poincaré diagram as a single point in which every one of the thousands of points plotted coincides exactly), continuous curves indicate quasi-periodicity (see Figure 6), and filled in areas potentially chaos (not found in any solutions explored) according to Lynch 2007 ^[11].

We do not present a Poincaré diagram for a periodic case, as the Poincaré bifurcation diagrams capture this behavior, and there is little to be gained from a picture of a single point. All periodic cases found in the analysis of this model have the same period as the forcing function or that period divided by a natural number (most likely the former from inspecting many of the cases in Figure 4), as all periodic cases have only one distinct point. Multiple points (not found here) would instead indicate period doubling, tripling, or...etc. from the period of the forcing function. We found several cases that exhibited quasi-periodicity and depict two examples in Figure 6.

The bifurcation diagrams of Poincaré sections presented in Figure 7 take a wide range of Poincaré sections, plotting only the grazer densities $y$ for a given parameter set along the $x$-axis, showing for what combinations of $K_0$ and $\varepsilon$ periodicity, quasi-periodicity, or extinction occur. That is each point from a Poincaré diagram is mapped from $(x, y)$ to $(K_0, y)$ for the value of $K_0$ that corresponds to that diagram, and all data having to do with the producer is disregarded. This is performed for many different individual Poincaré diagrams where each one shares the same value of $\varepsilon$. The $x$-axis is sampled at $0.01$ increments of $K_0$. Thus, there will be $300$ Poincaré diagrams represented in a Poincaré bifurcation diagram which varies $K_0$ from $0$ to $3$.

Any sample with only one dot at a given value of $K_0$ represents a periodic solution (the positioning of that dot along the $y$-axis has no physical significance other than that the grazer population reaches that density at some point during a given half year), while any sample with multiple points along the $y$-axis represents a quasi-periodic solution. Multiple isolated spikes occur in the $\varepsilon = 0.5, 0.7$ cases that are only single points. These too represent very small regions of quasi-periodicity that may be smaller than the $K_0$ step size allows for multiple points along the $y$-axis to show as in the $\varepsilon = 0.7, K_0 = 1.65$ case. However, taking a finer step size does reveal that hidden in these spikes are indeed small ranges of $K_0$ for which quasi-periodicity does occur.

For low level of seasonality ($\varepsilon = 0.1$) the systems exhibit a wide range of quasi-periodicity: $K_0\in[0.57, 0.88]$ (see Figure 7a). As seasonality increases this zone of quasi-periodicity decreases to $K_0\in[0.6, 0.77]$ (see Figure 7b); however, smaller isolated regions of quasi-periodicity emerge for higher light conditions ($\varepsilon = 0.7, K_0 = 1.65$ for example). The amplitudes ($K_0$ values) of these isolated occurrences of quasi-periodicity increase as seasonality increases (see Figure 7c-d).

For low $\varepsilon$ value, the grazer is able to persist only for $K_0 \in (0.15, 1.92)$. As $\varepsilon$ is increased, it manages to persist only for $K_0 \in (0.16, 1.99)$. As $\varepsilon$ continues to increase, it is able to persist only for $K_0 \in (0.17, 2.17)$. As $\varepsilon$ increases still further, the population is capable of persisting for $K_0 \in (0.22, 2.52)$. One can observe that increasing the value of forcing amplitude ($\varepsilon$) increases the parameter space for which the predatory population persists.

4.   Conclusion and Future Directions

In the present study, a predator-prey model with periodically varying light-dependent carrying capacity is proposed and the effects of amplitude of seasonal forcing ($\varepsilon$) on the dynamics of the aquatic food chain are studied. Conditions for positivity and boundedness of solutions are obtained in Lemma 2.1 and Theorem 2.1, respectively. The simulations shown in figures 3 and 4 show seasonal patterns in algae and Daphnia population densities. These types of patterns are observed in natural systems. Indeed, a commonly observed phenomenon in seasonal cycles of plankton communities is the clear-water phase ^[12], which correspond to periods of time with low algal densities. Additionally, peaks in zooplankton communities have been observed to occur seasonally, with peaks in the Spring and the fall for eutrophic lakes ^[13].Bifurcation analysis of the model with average value of light-dependent carrying capacity and the amplitude of seasonal forcing is performed. Numerical experiments reveal that in the absence of seasonal variation (LKE model), the dynamics of system are either stable equilibrium or limit cycles. However, the modified predator-prey system with periodically varying light-dependent carrying capacity exhibits rich dynamics, including periodic and quasi-period solutions. Also, it can be observed that increasing the value of the forcing amplitude increases the probability of persistence for the grazers.

While our model captures the dynamics of incorporating light-dependent seasonality into producer carrying capacity, seasonal effects can influence many other aspects of population dynamics. In addition to the effects of light, the effects of temperature are a considerable opportunity for research in the near future, because temperature also varies periodically during the year. So, by following the procedure used by Scheffer et al. 1997 ^[12], the effects of temperature can also be included in the Seasonal LKE model (SLKE).

It is important to note that the SLKE model does not explicitly track phosphorus in the prey nor in the media that supports the prey. Therefore, it can be extended by mechanistically deriving and accounting for phosphorus in both the prey and the media following the procedure used by Wang et al. 2008 ^[14]. In addition to these two directions, we also note the opportunity to extend the SLKE model to include multiple consumer species and to examine subsequent impacts of coexistence and exclusion, as in the analysis of Loladze et al. 2000 ^[15].

Acknowledgments

We would like to thank the referees and the editor for their careful reading of the original manuscript and many valuable comments and suggestions that greatly improved the presentation of this work. AP was partially supported by NSF grant DMS-1615697.

A.   Proof of Lemma 2.1

Proof. We will use the technique of proof by contradiction, i.e., assume there exists a time $t_1>0$ such that a trajectory with initial conditions in the open rectangle $(0, k) \times (0, \frac{P}{\theta})$ touches or crosses the boundary of the closed rectangle $[0, k] \times [0, \frac{P}{\theta}]$ for the first time. We will have four cases to consider because the boundary consists of four sides.

Case 1 left border: $x(t_1) = 0$. Let $\overline{f} = f'(0) = \lim_{x\to 0} \frac{f(x)}{x}$ and $\overline{y} = \max _{t \in [0, t_1]} y(t) < \frac{P}{\theta}$. Then for every $t \in[0, t_1]$

where $\alpha$ is a constant. Thus, $x(t) \geq x(0)e^{\alpha t}$, which implies that $x(t_1) \geq x(0)e^{\alpha t_1}>0$. Therefore, no trajectory can touch the left-hand-side border of the rectangle.

Case 2 right border: $x(t_1) = k$. Then for every $t \in[0, t_1]$

By the standard comparison argument we find $x(t) < k$, thus, no trajectory touches the right-hand-side border.

Case 3 bottom border: $y(t_1) = 0$. Then for every $t \in[0, t_1]$

This implies that $y(t_1) \geq y(0)e^{-d t_1}>0$. Therefore, no trajectory can touch the bottom border.

Case 4 top border: $y(t_1) = \frac{P}{\theta}$. $0 \leq y < P/ \theta$ for $0 \leq t < t_1$. Let $\overline{f} = f'(0) = \lim_{x\to 0}$.

Then for every $t \in[0, t_1]$

The standard comparison argument yields a contradiction, $y(t) < \frac{P}{\theta}$, thus, no trajectory touches the top border.

Conflict of interest

All authors declare no conflicts of interest in this paper.