Full article: Contrasting stoichiometric dynamics in terrestrial and aquatic grazer

Formulae display: $MathJax Logo$ ?Mathematical formulae have been encoded as MathML and are displayed in this HTML version using MathJax in order to improve their display. Uncheck the box to turn MathJax off. This feature requires Javascript. Click on a formula to zoom.

Abstract

The turnover rate of producer biomass in aquatic ecosystems is generally faster than in terrestrial. That is, aquatic producer biomass grows, is consumed, and is replaced considerably faster than terrestrial. The WKL model describes the flow of phosphorus and carbon through a grazer–producer system, hence varying the model parameters allows for analysis of different ecosystems of this type. Here we explore the impacts of the intrinsic growth rate of the producer and the maximal ingestion rate of the grazer on these dynamics because these parameters determine turnover rate. Simulations show that for low intrinsic growth rate and maximal ingestion rate, the grazer goes extinct; for higher values of these parameters, coexistence occurs in oscillations. Sensitivity analysis reveals the relative importance of all parameters on asymptotic dynamics. Lastly, the impacts of changing these two parameters in the LKE model appears to be quantitatively similar to the impacts in the WKL model.

Keywords:

2010 Mathematics Subject Classifications:

1. Introduction

Earth has many diverse habitats. One of the more fundamental dichotomies in the biosphere is between aquatic- and terrestrial-based ecosystems. These two groups vary greatly in average scale, both spatially and temporally. For example, consider the difference between an acacia tree, which grows at a rate of 44.2 cm per year [Citation13] and feeds giraffes with an average individual ingestion rate of 16.6–19.0 kg per day [Citation14], and a microscopic phytoplankton species which supports zooplankton grazers, both of which cannot be individually discerned by the naked human eye. These two systems also vary in time scale, particularly in the rate of turnover of the producer species at the base of these food chains. While a cyanobacteria bloom on the surface of a lake can appear in a matter of days or weeks, trees take years to reach maturity. But there are still fundamental similarities between these two seemingly opposite systems which allow us to compare them at the elemental level.

The framework of ecological stoichiometry was developed to study such interactions. Ecological stoichiometry applies the law of conservation of mass to ecological interactions and studies the balance of the elements that make up life [Citation16]. This approach quantifies relationships between organisms made up of several measurable elements. These elements can be neither created nor destroyed in ordinary chemical reactions, imposing a balance on the amount in a closed biological system throughout its interactions and processes. Consideration of this balance in mathematical modelling of ecological systems allows for study of the flow of nutrients and energy in a system.

There are many elements that are required for growth, reproduction, and survival of living organisms. The three main elements in biological structural molecules are carbon, nitrogen, and phosphorus, despite their scarcity in the Earth's crust relative to other elements. Carbon provides energy as well as structure, while nitrogen and phosphorus are crucial constituents of proteins and nucleic acids [Citation16]. However, due to different requirements for biomolecules for structural, metabolic, and reproductive components, the stoichiometric ratios of organisms vary. In general, herbivores are assumed to have higher, more rigid nutrient requirements than the producer they consume [Citation16]. This nutrient imbalance means that grazers can be limited either by the quantity or quality of their food. This adds a degree of complexity to modelling grazing.

A common application of ecological stoichiometry is in the study of grazer–producer systems. Initially, these were modelled using the Lotka–Volterra predator–prey equations [Citation4]: $\begin{aligned} \frac{d x}{d t} & = b x - a x y \\ \frac{d y}{d t} & = c x y - d y \end{aligned}$ where x and y correspond to the prey and predator respectively, b is the net growth rate of the prey in the absence of predators, d is the net death rate of the predators in the absence of prey, and c/a is the conversion efficiency from prey to predator biomass (a>c).

The Lotka–Volterra predator–prey model has been widely used to model predator–prey and grazer–producer interactions, explaining examples such as the changes in fishing during WWI that sparked Volterra's interest in the topic, as well as cycles in the lynx and snowshoe hair pelts traded in the 1840s by the Hudson Bay company [Citation4]. However, there are cases where considering all prey (e.g. producers) to be identical at the elemental level fails to capture realistic dynamics.

For example, an experiment involving Daphnia and a green alga showed that at very high light intensity, the algal population boomed due to increased photosynthetic rate, but the grazer abundance remained low [Citation5]. While the light-limited growth of the algae at low light levels and the resulting low grazer abundance could be explained by the Lotka–Volterra equations, the model cannot explain a case where high algal abundance does not result in a high grazer abundance. This is because in any non-stoichiometric form of the Lotka–Volterra equations, increased algal growth can only be beneficial to the grazer. But, as experimentally demonstrated, this is not always true – when there is too much growth of the algae, their more flexible nutrient requirements allow them to become phosphorus-poor, and therefore they can limit grazer growth due to being poor quality food relative to the requirement of the grazer.

One model developed using the framework of ecological stoichiometry to deal with this counterexample is the LKE model [Citation9], which is a predator–prey model. In this case, the ‘prey’ is a primary producer, such as phytoplankton, and the ‘predator’ is a zooplankton grazer, such as Daphnia. This model tracks only two elements, carbon (C) and phosphorus (P), where all others are assumed to be sufficiently abundant so as to be nonlimiting – that is, there is enough in the environment for the requirements of the organisms considered. Carbon is often included in ecological stoichiometry models, since it can be used to represent energy or biomass. The producer population is quantified by the density of carbon in the producer, x, and the grazer population by the density of carbon in the grazer, y. In this case, phosphorus was chosen as a focal nutrient since it is often a limiting nutrient in freshwater systems, and it is used in construction of several biological molecules for structure and energy metabolism.

The LKE model [Citation9] is $\begin{aligned} \frac{d x}{d t} & = b x (1 - \frac{x}{m i n (K, (P - θ y) / q)}) - f (x) y \\ \frac{d y}{d t} & = \hat{e} m i n (1, \frac{(P - θ y) / x}{θ}) f (x) y - d y \end{aligned}$ where x and y are the producer and grazer carbon densities, b is the intrinsic growth rate of the producer, d is the specific loss rate of the grazer (including respiration and death), K is the constant light dependent carrying capacity, $\hat{e}$ is the maximal production efficiency, $f (x)$ is the grazer's ingestion rate, q is the minimum P:C in the producer, θ is the fixed P:C of the grazer, and P is the total phosphorus in the system.

In the absence of the grazer, the producer exhibits logistic growth limited either by energy or by the availability of phosphorus. On the other hand, the grazer carbon levels undergo exponential decay in the absence of the producer. The growth of the grazer is limited either by food quantity or food quality. That is, by either the amount of producer carbon available or by the amount of producer phosphorus available relative to their needs. This in particular allows for the model to exhibit such dynamics in the very high light case as those presented by Elser and Kuang (2002) [Citation5], where phosphorus limitation of the producer growth causes the algae to become poor quality food for the grazer, and thus to limit the grazer abundance.

This model's applications are limited by one of its main assumptions, which states that all phosphorus in the system is divided into two pools: phosphorus in the grazer and phosphorus in the producer. This requires immediate recycling of phosphorus and immediate utilization by the producer, and does not allow for any free phosphorus in the medium. The relaxation of this assumption yielded the WKL model [Citation18], which is presented in Section 2 and is the basis for this work.

Despite the difficulty in analysing the nonsmooth LKE model, some analysis has been completed. If f and g are assumed to be Holling type I functional responses, then the system has no limit cycles and the internal equilibrium is globally asymptotically stable [Citation8]. With Holling type II functional responses, bifurcation analysis of the parameter K revealed the potential for bistability and several bifurcations [Citation8]. A global analysis of the LKE model with Holling type II functional responses was also completed, revealing four types of bistability as well as many possible bifurcations [Citation21]. These analyses illuminate the rich and complicated dynamics this relatively simple stoichiometric model [Citation9] can exhibit.

There are several other possible applications of ecological stoichiometry, many of which involve extensions of the LKE model. These are primarily focussed on looking at the effect of food quality on population dynamics, since explicit consideration of more nutrients than carbon allows other nutrient limitations to impact the model dynamics in realistic ways. For example, there are models incorporating the stoichiometric knife edge – a theory that there is an ideal nutrient richness, due to evidence that grazers are affected by both insufficient and excess food nutrient content [Citation11, Citation12, Citation22]. There is also a model which explicitly tracks phosphorus loading of the aquatic environment [Citation1]. This is globally relevant because of nutrient loading by humans due to agricultural fertilization and industrial emissions.

There is also evidence that elemental mismatches between trophic levels can influence population growth and foraging behaviours. Ecological stoichiometry has been used to capture this influence through consideration of varied energetic costs of foraging dependent upon food nutrient content. Such models have been used to show that grazers can benefit from compensatory feeding behaviours when consuming non-optimal food [Citation10].

However, despite its global applications, stoichiometry comes with many associated challenges. As with all models, one must balance the realisticness of the model with the ease of analysis. Stoichiometric systems are often nonsmooth, and require consideration of multiple cases due to applications of Liebig's Law of the Minimum [Citation22].

Some analysis has already been completed for the WKL model [Citation18]. The analysis is focussed around K, the resource carrying capacity determined by light. This particular parameter is controllable in a laboratory setting. However, there are other parameters that remain to be investigated. These parameters help uniquely define the conditions both within and surrounding the biological interaction we consider.

This paper intends to compare the dynamics in a terrestrial ecosystem versus those in an aquatic ecosystem by specifically focussing on r and c. The maximal growth rate of producers, r, tends to be higher in aquatic systems than in terrestrial [Citation15]. Similarly, the maximum ingestion rate tends to be higher in aquatic grazer–producer systems than terrestrial. There is evidence that producer biomass can be consumed at a rate four times higher by aquatic grazers than terrestrial [Citation15]. Comparison of the life cycle of an acacia tree to that of a phytoplankton clearly exemplifies this phenomenon. Investigating these parameters can allow us to rigorously compare such terrestrial and aquatic ecosystems, despite the extensive biological differences.

Another related contrast in parameter values lies in the tradeoff between r and c. Often, prey species must ‘choose’ between investing energy into their growth, increasing r, or their defences against predation, decreasing c [Citation6, Citation17]. For example, some Caribbean coral reef sponges have been found to exhibit increased growth when they are not under the threat of predation, even if the predation is prevented by cages [Citation7]. Hence, we expect to naturally see producers with low r and low c, and producers with high r and high c. This is similar to the above contrast between terrestrial and aquatic ecosystems. Thus, the investigative efforts within this paper can also allow for comparison between organisms' survival/reproductive strategies.

2. Model formulation

The WKL model [Citation18] tracks the amounts of carbon and phosphorus in the producer and the grazer. The model uses the following assumptions:

The total mass of phosphorus in the entire system is fixed, i.e. the system is closed to phosphorus but open to carbon.
The phosphorus to carbon ratio (P:C) in the producer varies, but never falls below a fixed minimum q (mg P/mg C); the grazer maintains a constant P:C ratio, denoted by θ (mg P/mg C), which is known as homeostasis [Citation16].

Let x be the density of carbon content in the producer and y be the density of carbon content in the grazer, both measured in (mg C)/l. For phosphorus contents, we use p for the density of phosphorus content in the producer and P for the density of free phosphorus in the media, both measured in (mg P)/l. Hence, p/x is the P:C ratio of the producer. Note that due to assumption 2, we do not explicitly track the phosphorus content in the grazers – the instantaneous phosphorus content in the grazer is simply $θ y$ .

The resulting equations are $\begin{aligned} \frac{d x}{d t} & = \underset{p r o d u c e r g r o w t h l i m i t e d b y n u t r i e n t a n d l i g h t}{\underset{⏟}{r x (1 - \frac{x}{m i n {K, p / q}})}} - \underset{u p t a k e b y g r a z e r s}{\underset{⏟}{f (x) y}} \\ \frac{d y}{d t} & = \underset{g r a z e r g r o w t h l i m i t e d b y f o o d q u a l i t y a n d q u a n t i t y}{\underset{⏟}{\hat{e} m i n \{1, \frac{p / x}{θ}\} f (x) y}} - \underset{g r a z e r d e a t h a n d r e s p i r a t i o n l o s s}{\underset{⏟}{\hat{d} y}} \\ \frac{d p}{d t} & = \underset{P u p t a k e b y p r o d u c e r}{\underset{⏟}{g (P) x}} - \underset{P l o s s d u e t o g r a z i n g}{\underset{⏟}{\frac{p}{x} f (x) y}} - \underset{P l o s s d u e t o p r o d u c e r r e c y c l i n g}{\underset{⏟}{d p}} \\ \frac{d P}{d t} & = \underset{P u p t a k e b y p r o d u c e r}{\underset{⏟}{- g (P) x}} + \underset{P r e c y c l i n g f r o m p r o d u c e r}{\underset{⏟}{d p}} \\ + \underset{P r e c y c l i n g f r o m d e a d g r a z e r}{\underset{⏟}{θ \hat{d} y}} + \underset{P r e c y c l i n g f r o m g r a z e r f a e c e s}{\underset{⏟}{(\frac{p}{x} - \hat{e} m i n \{θ, \frac{p}{x}\}) f (x) y}} \end{aligned}$ From assumption 1, the total phosphorus in the system is constant, i.e. $\begin{aligned} \frac{d T}{d t} & = 0 \end{aligned}$ for $T = p + P + θ y$ . Thus, we can write $P = T - p - θ y$ , and we can reduce the system to three equations [Citation18]: (1) $\begin{aligned} \frac{d x}{d t} & = r x (1 - \frac{x}{m i n {K, p / q}}) - f (x) y \end{aligned}$ (1) (2) $\begin{aligned} \frac{d y}{d t} & = \hat{e} m i n \{1, \frac{p / x}{θ}\} f (x) y - \hat{d} y \end{aligned}$ (2) (3) $\begin{aligned} \frac{d p}{d t} & = g (T - p - θ y) x - \frac{p}{x} f (x) y - d p \end{aligned}$ (3) The parameters are r, the intrinsic growth rate of the resource (day $^{- 1}$ ); K, the resource carrying capacity determined by light ((mg C)/l); q, the minimal possible P:C of the producer ((mg P)/(mg C)); $\hat{e}$ , the maximal conversion rate of the grazer; θ, the homeostatic constant P:C of the grazer; $\hat{d}$ , the loss rate of the grazer (day $^{- 1}$ ); and d, the phosphorus loss rate of the producer (day $^{- 1}$ ). Due to the second law of thermodynamics, $\hat{e} < 1$ , and in reality, $θ > q$ . The model also uses two functions: f(x), which is the ingestion rate of the grazers, and g(P), which is the per capita P uptake by the producer: $\begin{aligned} f (x) & = \frac{c x}{a + x} \\ g (P) & = \frac{\hat{c} P}{\hat{a} + P} \end{aligned}$ Here both f(x) and g(P) are assumed to take the form of Holling type II functional response. In general, f and g are assumed to be bounded smooth functions that are zero at zero, strictly increasing, and concave down. For the Holling type II functional responses chosen here, let c be the maximal ingestion rate of the grazer (day ⁻¹ ); $\hat{c}$ , the maximal phosphorus uptake rate of the producer ((mg P)/(mg C)/ day); a, the half-saturation constant of the grazer ((mg C)/l); and $\hat{a}$ , the phosphorus half-saturation constant of the producer ((mg P)/l).

For analyses, parameter values from Wang et al. (2008) [Citation18] are used, as shown in Table . The baseline values used for r and c are 0.93 day $^{- 1}$ and 0.75 day $^{- 1}$ respectively. Therefore, we consider values ranging from 0.1 to 2.0 for these two parameters, since this provides us up to more than double the realistic parameters used previously. With these analyses, all other parameters are assumed to be equal between terrestrial and aquatic ecosystems, which clearly somewhat limits the applicability of the results. The system was originally parametrized for a freshwater system. Hence, we consider the parameter values for r and c less than the specified baseline values to represent terrestrial ecosystems, and those greater than or equal to the baselines to represent aquatic.

Table 1. The parameter (P) values (V) used for simulations [Citation18].

Display Table

3. Mathematical analysis

3.1. Invariant set

Wang et al. (2008) [Citation18] presented the following theorem for forward invariance.

Theorem 3.1

Solutions with initial conditions in the set $Ω = {(x, y, p) : 0 < x < min {K, T / q}, 0 < y, 0 < p, p + θ y < T}$ remain there for all forward times.

This is proven by way of contradiction. The full proof was provided by Wang et al. (2008) [Citation18] but in brief, one considers a solution X(t) with initial condition in Ω, then assumes there is a time $t_{1}$ such that X(t) touches or crosses the boundary of the closure of Ω for the first time. Then one considers cases for different segments of boundary and reaches a contradiction in each case.

This set is biologically meaningful. Densities of elements cannot be negative, so we require positivity, which we have here. Also, we know that growth of the producer is limited either by light (K) or by phosphorus availability relative to their needs (T/q). This is Liebig's Law of the Minimum – the resource which is least abundant relative to an organism's needs becomes limiting [Citation16]. Hence, the bounds on x in Ω make sense. The bounds on phosphorus levels also make sense: the phosphorus contained in the producers and grazers ( $p + θ y$ ) cannot exceed what is available in the system (T).

For the purposes of this paper, this set should be kept in mind when considering stability. Equilibria outside of this set cannot be globally attracting, since no solution starting in this set will ever leave it. Also, for numerical simulations, once a solution enters this set, we know the general location of the solution for all forward times.

3.2. Equilibria

Wang et al. (2008) [Citation18] found the equilibria for the model when f and g are Holling type I functions, then analysed the stability of the boundary steady states. For this simplified case, there were two boundary equilibria: the extinction equilibrium $E_{0}$ = (0, 0, 0), and the grazer extinction equilibrium $E_{1}$ where the form depends on if light or nutrients are limiting for the producer. $E_{1}$ is given by $\begin{aligned} (K, 0, \frac{T K}{K + d / α}) K < \frac{T}{q} - \frac{d}{α} \\ (\frac{T}{q} - \frac{d}{α}, 0, q (\frac{T}{q} - \frac{d}{α})) K > \frac{T}{q} - \frac{d}{α} \end{aligned}$ where f(x) = $β x$ and g(P) = $α P$ [Citation18].

As stated in Section 2, for the model, f and g are always assumed to be bounded smooth functions that are zero at zero, strictly increasing, and concave down. While the Holling type I functional responses used in [Citation18] do satisfy this requirement, they are not realistic. Holling type I requires the assumption that there are no physical limits to the amount of food the grazer can consume. Clearly metabolic restrictions make this unrealistic. Thus, while Holling type I makes mathematical analysis more manageable, it limits applicability of the results.

For the numerical analyses, we assume that f and g are Holling type II functions. Hence, we find equilibria that will match our numerical analyses. Then, the equilibria $(\bar{x}, \bar{y}, \bar{p})$ satisfy (4) $\begin{aligned} 0 & = \bar{x} (r (1 - \frac{\bar{x}}{m i n {K, \bar{p} / q}}) - \frac{c \bar{y}}{a + \bar{x}}) \equiv \bar{x} F (\bar{x}, \bar{y}, \bar{p}) \end{aligned}$ (4) (5) $\begin{aligned} 0 & = \bar{y} (\hat{e} m i n \{1, \frac{\bar{p} / \bar{x}}{θ}\} \frac{c \bar{x}}{a + \bar{x}} - \hat{d}) \equiv \bar{y} G (\bar{x}, \bar{y}, \bar{p}) \end{aligned}$ (5) (6) $\begin{aligned} 0 & = \frac{\hat{c} (T - \bar{p} - θ \bar{y}) \bar{x}}{\hat{a} + (T - \bar{p} - θ \bar{y})} - \frac{c \bar{p} \bar{y}}{a + \bar{x}} - d \bar{p} \equiv H (\bar{x}, \bar{y}, \bar{p}) \end{aligned}$ (6) To find the equilibria, we split into cases based on the minimum functions included in $F (\bar{x}, \bar{y}, \bar{p})$ and $G (\bar{x}, \bar{y}, \bar{p})$ , given by

$\bar{p} < K q$ and $\bar{p} < θ \bar{x}$ : producer is nutrient limited and grazer is limited by food quality.
$\bar{p} < K q$ and $\bar{p} > θ \bar{x}$ : producer is nutrient limited and grazer is limited by food quantity.
$\bar{p} > K q$ and $\bar{p} < θ \bar{x}$ : producer is light limited and grazer is limited by food quality.
$\bar{p} > K q$ and $\bar{p} > θ \bar{x}$ : producer is light limited and grazer is limited by food quantity.

All four cases have between one and three boundary equilibria, the extinction equilibrium $E_{0} = (0, 0, 0)$ , and the grazer extinction equilibrium/equilibria $E_{1}$ . As with the Holling type I case, the form of $E_{1}$ depends on what resource is limiting for the producer, i.e. it is the same for case 1 as case 2, and the same for case 3 as case 4. When the producer is nutrient limited, E $_{1}$ is $(\frac{d q (\hat{a} + T) - \hat{c} T}{q (d q - \hat{c})}, 0, \frac{d q (\hat{a} + T) - \hat{c} T}{d q - \hat{c}})$ Note that for the above, $\bar{x} = (\bar{p} / q)$ . Also, for the given baseline parameter values, this equilibria is non-negative and thus biologically feasible ( $E_{1} = (7.4980, 0, 0.0300)$ ). When the producer is light limited, $E_{1}$ actually has two possible values of $\bar{p}$ , both of which are positive: $(K, 0, \frac{(\hat{a} d + T d + \hat{c} K) \pm \sqrt{(\hat{a} d + T d + \hat{c} K)^{2} - 4 d \hat{c} T K}}{2 d})$ For cases 1 and 3, we also have an additional mathematically possible boundary equilibrium: $(0, - \frac{a d}{c}, \frac{\hat{d} a θ}{\hat{e} c})$ However, this is not biologically feasible. All parameters are assumed to be positive and thus for this equilibrium, $\bar{y} < 0$ which is not realistic since we cannot observe negative densities of carbon. Also, this equilibrium would not satisfy the quality limitation condition if we multiply the terms within the minimum by $f (x)$ , since $f (0) = 0 < ((c \bar{p}) / (a θ)) = \hat{d} / \hat{c}$ .

For the parameter values used in the numerical simulations, $K < (T / q) = 7.50$ for all K in 0.25 –2.00. Also, for K<7.5, $K q < \bar{p}$ for either form of the boundary equilibrium. Thus for the range of K we consider and the values of T and q used, we will never have $K q > \bar{p}$ and thus these equilibria will always fall in the producer light limited region of the phase space. Hence we are restricted to cases 3 and 4 for the parameters given in Table 1. For the condition that distinguishes case 3 from case 4, we note that $\bar{x} θ = K θ$ . Also, $\bar{p}$ does not depend on r or c, thus we can fix r and c at their baseline values and check the condition only varying K. Therefore, we only need to check the sign of $\bar{p} - θ \bar{x}$ for our various values of K, and for both of the grazer extinction equilibria (varying $\bar{p}$ ). For the version of $\bar{p}$ that uses the plus sign, we observe that $\bar{p} - θ \bar{x} > 0$ for all $K \in (0, 2]$ . Hence, this boundary equilibrium is always in case 4. On the other hand, the equilibrium applying the negative sign has $\bar{p} - θ \bar{x} > 0$ until a value of K between 0.74 and 0.75, when it becomes positive. Hence, this boundary equilibrium is in case 4 until approximately K = 0.75, at which point it switches to case 3.

Observe that in all cases, the forms of the biologically feasible boundary equilibria have no explicit dependence on either r or c. Given our assumption that all other parameters are the same between terrestrial and aquatic ecosystems, this means that the value of the boundary equilibria will not depend on whether the ecosystem is terrestrial or aquatic. However, the asymptotic state of the system will still depend on r and c, since they will determine which equilibrium is stable.

Lastly, there may exist coexistence equilibria, which satisfy $\begin{aligned} 0 & = r (1 - \frac{\bar{x}}{m i n {K, \bar{p} / q}}) - \frac{c \bar{y}}{a + \bar{x}}, \\ 0 & = \hat{e} m i n \{1, \frac{\bar{p} / \bar{x}}{θ}\} \frac{c \bar{x}}{a + \bar{x}} - \hat{d}, \\ 0 & = \frac{\hat{c} (T - \bar{p} - θ \bar{y}) \bar{x}}{\hat{a} + (T - \bar{p} - θ \bar{y})} - \frac{c \bar{p} \bar{y}}{a + \bar{x}} - d \bar{p}, \end{aligned}$ i.e. $F (\bar{x}, \bar{y}, \bar{p}) = 0$ , $G (\bar{x}, \bar{y}, \bar{p}) = 0$ , and $H (\bar{x}, \bar{y}, \bar{p}) = 0$ . Clearly there is likely to be an explicit dependence on r and c for coexistence equilibria, although analytically determining the explicit dependence is incredibly time consuming and may be impossible due to the extensive nonlinearity.

The results of this section are summarized in the following theorem.

Theorem 3.2

Equations (Equation1(1) $\begin{aligned} \frac{d x}{d t} & = r x (1 - \frac{x}{m i n {K, p / q}}) - f (x) y \end{aligned}$ (1) )–(Equation3(3) $\begin{aligned} \frac{d p}{d t} & = g (T - p - θ y) x - \frac{p}{x} f (x) y - d p \end{aligned}$ (3) ) with Holling type II functional responses possess the trivial extinction equilibrium $E_{0} = (0, 0, 0)$ , up to two boundary grazer extinction equilibria $E_{1}$ , and may have coexistence equilibria, where the boundary equilibria $E_{1}$ have two possibilities:

If $\bar{p} < K q,$ then there is one boundary equilibrium $E_{1} = (\frac{d q (\hat{a} + T) - \hat{c} T}{q (d q - \hat{c})}, 0, \frac{d q (\hat{a} + T) - \hat{c} T}{d q - \hat{c}})$
If $\bar{p} > K q,$ then there are two boundary equilibria $E_{1} = (K, 0, \frac{(\hat{a} d + T d + \hat{c} K) \pm \sqrt{(\hat{a} d + T d + \hat{c} K)^{2} - 4 d \hat{c} T K}}{2 d})$

3.3. Stability

For the extinction equilibrium, Wang et al. (2008) [Citation18] proved a stability theorem for the extinction steady-state. Here we improve upon the theorem, while following a very similar proof:

Theorem 3.3

The extinction steady-state $E_{0} = (0, 0, 0)$ in Equations (Equation1(1) $\begin{aligned} \frac{d x}{d t} & = r x (1 - \frac{x}{m i n {K, p / q}}) - f (x) y \end{aligned}$ (1) )–(Equation3(3) $\begin{aligned} \frac{d p}{d t} & = g (T - p - θ y) x - \frac{p}{x} f (x) y - d p \end{aligned}$ (3) ) is globally asymptotically stable if $d > \tilde{m} g (T),$ where $\tilde{m} = m i n {x (0) / p (0), [1 + d / r] / q}$ .

Proof.

Let $u = x / p$ , then applying quotient rule as well as Equations (Equation1(1) $\begin{aligned} \frac{d x}{d t} & = r x (1 - \frac{x}{m i n {K, p / q}}) - f (x) y \end{aligned}$ (1) )–(Equation3(3) $\begin{aligned} \frac{d p}{d t} & = g (T - p - θ y) x - \frac{p}{x} f (x) y - d p \end{aligned}$ (3) ) $\begin{aligned} \frac{d u}{d t} & = \frac{d}{d t} \frac{x}{p} \\ = \frac{(d x / d t) p - x (d p / d t)}{p^{2}} \\ = \frac{d x / d t}{p} - x \frac{d p / d t}{p^{2}} \\ = \frac{1}{p} (r x (1 - \frac{x}{m i n {K, p / q}}) - f (x) y) - \frac{x}{p^{2}} (g (T - p - θ y) x - \frac{p}{x} f (x) y - d p) \\ = \frac{r x}{p} (1 - \frac{x}{m i n {K, p / q}}) - \frac{f (x) y}{p} - \frac{x^{2}}{p^{2}} g (T - p - θ y) + \frac{f (x) y}{p} + \frac{d x}{p} \\ = r u (1 - \frac{x}{m i n {K, p / q}}) - \frac{f (x) y}{p} - u^{2} g (T - p - θ y) + \frac{f (x) y}{p} + d u \\ = r u (1 - \frac{x}{m i n {K, p / q}}) - u^{2} g (T - p - θ y) + d u \end{aligned}$ Note that $m i n {K, p / q} \leq p / q ⟺ \frac{1}{p / q} \leq \frac{1}{m i n {K, p / q}} ⟺ - \frac{1}{m i n {K, p / q}} \leq - \frac{1}{p / q}$ Also, since $g (0) = 0$ and $g^{'} (P) > 0$ $- u^{2} g (T - p - θ y) \leq 0$ Thus $\frac{d u}{d t} \leq r u (1 - \frac{x}{p / q}) + d u = r u (1 - q u) + d u$ Then $\frac{d u}{d t} \leq r u (1 + d / r - q u)$ Therefore, $u \leq m i n {x (0) / p (0), [1 + d / r] / q} \equiv \tilde{m}$ . From Equation (Equation3(3) $\begin{aligned} \frac{d p}{d t} & = g (T - p - θ y) x - \frac{p}{x} f (x) y - d p \end{aligned}$ (3) ) $\begin{aligned} \frac{d p}{d t} & = g (T) x - \frac{p}{x} f (x) y - d p \\ \leq g (T) x - d p \\ \leq g (T) \tilde{m} p - d p \\ = (g (T) \tilde{m} - d) p \end{aligned}$ since $u = x / p$ , and $u \leq \tilde{m}$ . Since we assume $d > \tilde{m} g (T)$ , and this implies $g (T) \tilde{m} - d < 0$ , then $p \to 0$ as $t \to \infty$ .

Now, consider Equation (Equation1(1) $\begin{aligned} \frac{d x}{d t} & = r x (1 - \frac{x}{m i n {K, p / q}}) - f (x) y \end{aligned}$ (1) ) $\begin{aligned} \frac{d x}{d t} & = r x (1 - \frac{x}{m i n {K, p / q}}) - f (x) y \\ \leq r x (1 - \frac{x}{p / q}) = r x (1 - \frac{q x}{p}) \end{aligned}$ Therefore, ${l i m s u p}_{t \to \infty} x (t) \leq p / q$ . Since $p \to 0$ as $t \to \infty$ , then this implies $x \to 0$ as $t \to \infty$ . Then as $t \to \infty$ in Equation (Equation2(2) $\begin{aligned} \frac{d y}{d t} & = \hat{e} m i n \{1, \frac{p / x}{θ}\} f (x) y - \hat{d} y \end{aligned}$ (2) ), the first term goes to 0 and thus $y \to 0$ .

Thus, the extinction steady-state is globally asymptotically stable if $d > \tilde{m} g (T)$ .

The original theorem from Wang et al. (2008) [Citation18] had $d > m g (T)$ with $m = m i n {x (0) / p (0), [1 + (d + f^{'} (0) T / θ) / r] / q}$ . Since we assume f is strictly increasing and all parameters are assumed to be positive, then $f^{'} (0) T / θ > 0$ , and thus $1 + (d + f^{'} (0) T / θ) / r] / q > [1 + d / r] / q .$ Therefore, $m \geq \tilde{m}$ , and so there is a potentially larger range of values of d for which the extinction equilibrium is stable if we use $\tilde{m}$ instead of m.

This theorem was proven for the general forms of f and g, and thus applies for Holling type II. We observe the possible dependence on r and c. However, for the parameter regime considered here, this condition requires $m = x (0) / p (0) < 0.3167$ , which is highly unrealistic as it requires the initial cell quota of the producer to exceed 3.1579. However, this condition is not necessary, and we may still observe stability of the extinction steady state.

It remains to investigate the stability of the other boundary equilibria. The Jacobian matrix is $A = [\begin{matrix} F + x F_{x} & x F_{y} & x F_{p} \\ y G_{x} & G + y G_{y} & y G_{p} \\ H_{x} & H_{y} & H_{p} \end{matrix}]$ By Routh–Hurwitz criterion, all eigenvalues of A have strictly negative real parts if the following conditions hold: trA < 0; detA < 0; and detA - (trA)( $\sum_{k = 1}^{3} A_{k k}$ ) > 0, where $A_{k k}$ is the determinant of the matrix produced by removing row k and column k from matrix A.

We compute the Jacobian by computing the necessary terms: $\begin{aligned} \frac{\partial F}{\partial x} & = - \frac{r}{m i n {K, p / q}} + \frac{c y}{(a + x)^{2}} & F + x F_{x} & = \{\begin{cases} r - \frac{2 r x}{K} - \frac{a c y}{(a + x)^{2}} & K < p / q \\ r - \frac{2 r q x}{p} - \frac{a c y}{(a + x)^{2}} & K > p / q \end{cases} \\ \frac{\partial F}{\partial y} & = - \frac{c}{a + x} & x F_{y} & = - \frac{c x}{a + x} \\ \frac{\partial F}{\partial p} & = \{\begin{cases} 0 & K < p / q \\ \frac{r x q}{p^{2}} & K > p / q \end{cases} & x F_{p} & = \{\begin{cases} 0 & K < p / q \\ \frac{r q x^{2}}{p^{2}} & K > p / q \end{cases} \\ \frac{\partial G}{\partial x} & = \{\begin{cases} \frac{a c \hat{e}}{(a + x)^{2}} & 1 < \frac{p / x}{θ} \\ - \frac{\hat{e} p c}{θ (a + x)^{2}} & 1 > \frac{p / x}{θ} \end{cases} & y G_{x} & = \{\begin{cases} \frac{a c \hat{e} y}{(a + x)^{2}} & 1 < \frac{p / x}{θ} \\ - \frac{c \hat{e} p y}{θ (a + x)^{2}} & 1 > \frac{p / x}{θ} \end{cases} \\ \frac{\partial G}{\partial y} & = 0 & G + y G_{y} & = \{\begin{cases} \frac{c \hat{e} x}{a + x} - \hat{d} & 1 < \frac{p / x}{θ} \\ \frac{c \hat{e} p}{θ (a + x)} - \hat{d} & 1 > \frac{p / x}{θ} \end{cases} \\ \frac{\partial G}{\partial p} & = \{\begin{cases} 0 & 1 < \frac{p / x}{θ} \\ \frac{c \hat{e}}{θ (a + x)} & 1 > \frac{p / x}{θ} \end{cases} & y G_{p} & = \{\begin{cases} 0 & 1 < \frac{p / x}{θ} \\ \frac{c \hat{e} y}{θ (a + x)} & 1 > \frac{p / x}{θ} \end{cases} \\ \frac{\partial H}{\partial x} & = \frac{\hat{c} (T - p - θ y)}{\hat{a} + (T - p - θ y)} + \frac{c p y}{(a + x)^{2}} \\ \frac{\partial H}{\partial y} & = - \frac{\hat{a} \hat{c} θ x}{(\hat{a} + T - p - θ y)^{2}} - \frac{c p}{a + x} \\ \frac{\partial H}{\partial p} & = - \frac{\hat{a} \hat{c} x}{(\hat{a} + T - p - θ y)^{2}} - \frac{c y}{a + x} - d \end{aligned}$ Therefore, for the four cases described in Section 3.2, we have different Jacobian matrices for E₁.

CASE 1: $\bar{p} < K q$ and $\bar{p} < θ \bar{x}$ : $E_{1} = (\bar{p} / q, 0, \bar{p})$ $A_{1} = [\begin{matrix} - r & - \frac{c \bar{p}}{a q + \bar{p}} & \frac{r}{q} \\ 0 & \frac{c \hat{e} \bar{p} q}{θ (a q + \bar{p})} - \hat{d} & 0 \\ \frac{\hat{c} (T - \bar{p})}{\hat{a} + T - \bar{p}} & - \frac{\hat{a} \hat{c} θ \bar{p}}{q (\hat{a} + T - \bar{p})^{2}} - \frac{c \bar{p} q}{a q + \bar{p}} & - \frac{\hat{a} \hat{c} \bar{p}}{q (\hat{a} + T - \bar{p})^{2}} - d \end{matrix}]$ CASE 2: $\bar{p} < K q$ and $\bar{p} > θ \bar{x}$ : E₁ = $(\bar{p} / q, 0, \bar{p})$ $A_{2} = [\begin{matrix} - r & - \frac{c \bar{p}}{a q + \bar{p}} & \frac{r}{q} \\ 0 & \frac{c \hat{e} \bar{p}}{a q + \bar{p}} - \hat{d} & 0 \\ \frac{\hat{c} (T - \bar{p})}{\hat{a} + T - \bar{p}} & - \frac{\hat{a} \hat{c} θ \bar{p}}{q (\hat{a} + T - \bar{p})^{2}} - \frac{c \bar{p} q}{a q + \bar{p}} & - \frac{\hat{a} \hat{c} \bar{p}}{q (\hat{a} + T - \bar{p})^{2}} - d \end{matrix}]$ CASE 3: $\bar{p} > K q$ and $\bar{p} < θ \bar{x}$ : E₁ = $(K, 0, \bar{p})$ $A_{3} = [\begin{matrix} - r & - \frac{c K}{a + K} & 0 \\ 0 & \frac{c \hat{e} \bar{p}}{θ (a + K)} - \hat{d} & 0 \\ \frac{\hat{c} (T - \bar{p})}{\hat{a} + T - \bar{p}} & - \frac{\hat{a} \hat{c} θ K}{(\hat{a} + T - \bar{p})^{2}} - \frac{c \bar{p}}{a + K} & - \frac{\hat{a} \hat{c} K}{(\hat{a} + T - \bar{p})^{2}} - d \end{matrix}]$ CASE 4: $\bar{p} > K q$ and $\bar{p} > θ \bar{x}$ : E ₁ = $(K, 0, \bar{p})$ $A_{4} = [\begin{matrix} - r & - \frac{c K}{a + K} & 0 \\ 0 & \frac{c \hat{e} K}{a + K} - \hat{d} & 0 \\ \frac{\hat{c} (T - \bar{p})}{\hat{a} + T - \bar{p}} & - \frac{\hat{a} \hat{c} θ K}{(\hat{a} + T - \bar{p})^{2}} - \frac{c \bar{p}}{a + K} & - \frac{\hat{a} \hat{c} K}{(\hat{a} + T - \bar{p})^{2}} - d \end{matrix}]$ To determine the stability of $E_{1}$ for case i (i ∈{1, 2}) we need to find the trace and determinant of $A_{i}$ . However, this is not particularly illuminating mathematically – determining conditions such that $A_{1}$ and $A_{2}$ satisfy the Routh–Hurwitz criterion seems quite complicated (see Appendix).

However, for cases 3 and 4, we can decompose the matrix to be a $2 \times 2$ sub-matrix and one negative eigenvalue $- \frac{\hat{a} \hat{c} K}{(\hat{a} + T - \bar{p})^{2}} - d < 0$ . The $2 \times 2$ sub-matrix is upper triangle whose eigenvalues are $- r < 0$ and $A_{i}^{22}$ , i = 3, 4. The sign of the eigenvalue $A_{i}^{22}$ determines the stability of $E_{1}$ .

Since all parameters are positive, then $E_{1}$ is a stable node for $\hat{d} > \frac{c \hat{e} \bar{p}}{θ (a + K)}$ in case 3; and for $\hat{d} > \frac{c \hat{e} K}{a + K}$ in case 4.

In the opposite situation, $E_{1}$ is a saddle with a one-dimensional unstable manifold and a two-dimensional stable manifold for $\hat{d} < \frac{c \hat{e} \bar{p}}{θ (a + K)}$ in case 3; and for $\hat{d} < \frac{c \hat{e} K}{a + K}$ in case 4.

We can find the value of c such that $\hat{d}$ is equal to the right hand side in our conditions for a bifurcation to occur. The values are in Table . We observe that the sufficient condition for stability of our equilibria involves a low maximal ingestion rate, which suggests that terrestrial systems would be more likely to trend toward the grazer extinction equilibria based purely on r and c.

Table 2. The value of c for the appearance of a zero eigenvalue in cases 3 and 4. Note that there are two possible values of $\bar{p}$ for case 3, thus there is a row for each of these equilibria.

Display Table

The stability results for cases 3 and 4 are summarized in the following theorem.

Theorem 3.4

The following stability results hold for the grazer extinction equilibrium $E_{1}$ .

If $\bar{p} > K q$ and $\bar{p} < θ \bar{x}$ , the boundary equilibrium $E_{1}$ is a stable node for $\hat{d} > (c \hat{e} \bar{p}) / (θ (a + K)),$ is a saddle for $\hat{d} < (c \hat{e} \bar{p}) / (θ (a + K)),$ and a bifurcation occurs when $\hat{d} = (c \hat{e} \bar{p}) / (θ (a + K)) .$
If $\bar{p} > K q$ and $\bar{p} > θ \bar{x}$ , the boundary equilibrium $E_{1}$ is a stable node for $\hat{d} > (c \hat{e} K) / (a + K),$ is a saddle for $\hat{d} < (c \hat{e} K) / (a + K),$ and a bifurcation occurs when $\hat{d} = (c \hat{e} K) / (a + K) .$

4. Numerical dynamics and their implications

4.1. Numerical simulations

Using Matlab, the system was simulated using the parameter values in Table . Since the model is nonsmooth, ode45 was not reliable for certain combinations of parameters. A singularity produced negative densities, thus a stiff solver ode23s was used instead. The initial condition was always held at $(x_{0}, y_{0}, p_{0})$ = (0.3, 0.3, 0.01), which is not in the invariant set from Section 3.1 but is at least biologically feasible and orbits starting here can still enter the set in time. This initial condition was selected since it was the one used in the paper where the model was presented.

First the simulations were run for t ranging from 0 to 50 days, varying one of the two focal parameters at a time. For each K in {0.25, 0.75, 1.00, 2.00}, the system was numerically simulated for r in {0.1, 0.2,…, 2.0}, with c held constant at the baseline value 0.75. Then the process was repeated with r held constant at 0.93 and c in {0.1, 0.2,…, 2.0}. The lower values of r and c were used to represent the slower average turnover rate of terrestrial systems, and the higher parts of the ranges were used for the faster average turnover rate of aquatic systems.

For the intrinsic growth rate of the producer (r), we see that as r increases with K held at the following level:

K =0.25: grazers benefit; producers harmed until they plateau; coexistence at a steady state.
K =0.75, 1.00: oscillations appear and then replaced by coexistence at a steady state.
K =2.00: grazers benefit; producers harmed until a point where grazers go extinct.

For the maximal ingestion rate of the grazer (c), we see that as c increases with K held at the following level:

K =0.25, 0.75, 1.00: grazers benefit and producers harmed until oscillations appear.
K =2.00: grazers do better until the lines cross and the system trends toward oscillations.

Figures and show the shifts in dynamics described above.

Figure 1. Sample dynamics for various r values, r = 0.1, 1.0, 2.0, for $K = 0.25, 0.75, 1.00, 2.00$ and baseline c (as well as all other parameters). For low to intermediate light, we observe that increasing r in general increases the population densities for both producers and grazers. In the high and very high light cases, increasing r to 2.0 is detrimental to the grazer, likely due to nutrient limitation of grazer growth.

Figure 2. Sample dynamics for various c values, c = 0.1, 1.0, 2.0, for $K = 0.25, 0.75, 1.00, 2.00$ and baseline r (as well as all other parameters). In general, increasing c is detrimental to the producer population, and beneficial for the grazer in all cases except low light.

Then the simulations were run for t in 0 to 200, varying both of the focal parameters. The time limit was extended since there appeared to be some dynamics that had not completely ‘settled’ by t = 50. The values of the focal parameters were selected based on the previous simulations to align with where dynamics shifts occurred. Thus the simulations were run for all possible combinations of K in {0.25, 0.75, 1.00, 2.00}, r in {0.1, 0.5, 1.0, 1.5, 2.0}, and c in {0.1, 0.5, 1.0, 1.5, 2.0}, again using ode23s in Matlab. The behaviour at t = 200 was then classified according to the possible dynamics seen in prior papers [Citation9, Citation18]. The dynamics were classified as either grazer extinction; coexistence at a nonzero steady state; coexistence with oscillations; or coexistence with oscillations with decreasing amplitude, leading to coexistence at a steady state. Examples of the above dynamics are shown in Figure . These observed dynamics were then used to develop classifications for two-parameter bifurcation diagrams and expectations for what dynamics may be produced.

Figure 3. Sample dynamics: (a) grazer extinction (terrestrial r and c); (b) coexistence oscillations to coexistence steady state (mixed r and c); (c) coexistence oscillations (aquatic r and c); (d) coexistence at a steady state (mixed r and c).

4.2. Sensitivity analysis

Local sensitivity analysis was completed for the various light intensities. For K = 0.25, 1.00, and 2.00, the baseline parameters produce a solution which tends to an equilibrium. After visually assessing the dynamics for each of the required parameter combinations using a plot in Matlab, the variables selected for sensitivity analysis for these K values were x, y, and p at 200 days. The normalized forward sensitivity indices for $x (200)$ , $y (200)$ and $p (200)$ for all parameters were computed using the formula: $Υ_{ρ}^{u} \equiv \frac{\partial u}{\partial ρ} \times \frac{ρ}{u}$ where u is the variable and ρ is the parameter [Citation2]. To estimate the partial derivative, a central difference approximation was used: $\frac{\partial u}{\partial ρ} = \frac{u (p a r + h) - u (p a r - h)}{2 h} + O (h^{2})$ where u is the variable, ρ is the parameter, and par is the baseline values of the parameter, as given in Table . For the parameters that vary between simulations, r = 0.93 and c = 0.75 were used for the baseline values. To approximate the values of $x (200)$ , $y (200)$ , and $p (200)$ for par + h and par−h, ode23s was used, with all other parameters set at baseline values except the one for which the index was calculated. h was taken to be 1 percent of the baseline value. The results for $K = 0.25$ are shown in Table . This was also repeated for K = 1.00 and K = 2.00, with the results shown in Table and Table respectively.

Table 3. The sensitivity of variables to the parameters, sorted from largest to smallest absolute value, for low light (K = 0.25; ode23s).

Display Table

Table 4. The sensitivity of variables to the parameters, sorted from largest to smallest absolute value, for high light (K = 1.00; ode23s).

Display Table

Table 5. The sensitivity of variables to the parameters, sorted from largest to smallest absolute value, for very high light (K = 2.00; ode23s).

Display Table

For K = 0.75 (intermediate light), the baseline parameters produce a limit cycle. Accordingly, local sensitivity indices were computed instead for the amplitude and period of the oscillations in the variables, using the same formulas to estimate the partial derivative and the normalized forward sensitivity indices. Each of the necessary simulations was run using ode23s, and the timespan $[0, 200]$ was determined to be adequate to capture the settled behaviour. The period of the oscillations was estimated to be between 25 and 30 days, and in order to make sure that enough troughs and crests existed in the tail for each variable, the tail used was $[\geq 140, 200]$ , i.e. roughly the last 60 days, depending on the intervals established by ode23s. Amplitude was determined by finding the maximum and minimum values in the tail for the variable. The period was determined by first finding places where the sign of the change in the variable value from one element of the array to the next changed, then finding the time difference between the last and third last such point. The resulting sensitivity indices are shown in Table and Table for the amplitude and period respectively.

Table 6. The sensitivity of the amplitude of oscillations to the parameters, sorted from largest to smallest absolute value, for intermediate light (K = 0.75; ode23s).

Display Table

Table 7. The sensitivity of the period of oscillations to the parameters, sorted from largest to smallest absolute value, for intermediate light (K = 0.75; ode23s).

Display Table

For K = 0.25 (Table ), we see that the intrinsic growth rate of the producer, r, ranks fifth for $x (200)$ and $y (200)$ , and seventh for $p (200)$ . For K = 0.75 (Tables and ), it ranks sixth for the amplitude of oscillations x and seventh for the amplitude of y and p; it ranks thirds, fourth, and second respectively for the periods of oscillations of x, y, and p. For K = 1.00 (Table ), it ranks fifth for $x (200)$ , and seventh for both $y (200)$ and $p (200)$ . For K = 2.00 (Table ), it ranks eighth for $x (200)$ and $y (200)$ , and eleventh for $p (200)$ .

Comparatively, the grazer ingestion rate c ranks second, fourth, and fourth for K = 0.25 x, y, and p respectively. It ranks first for all amplitudes of oscillation for K = 0.75, and fifth, fifth, and sixth for the periods. For K = 1.00, c ranks third, third, and second. Finally, for K = 2.00, c ranks third, fourth, and fourth for x, y, and p respectively. The above rankings are summarized in Table .

Table 8. The ranking of r and c among the 12 parameters for the sensitivity indices. For K = 0.75, the rankings are amplitude first, then period.

Download CSV Display Table

The sign of the indices corresponds to the direction of the relationship between the parameter and the target value. We see that increasing r consistently increases the steady-state producer carbon density; for low light, it increases the steady-state grazer carbon density and decreases the steady-state producer phosphorus density, and vice versa for (very) high light. This suggests that increasing the intrinsic growth rate of the producer is consistently good for the producer, while it is only good for the grazer at low light and harmful for higher light conditions. This is likely due to the resulting decrease in nutrient quality of the producer with higher growth rates but limited phosphorus resources. At intermediate light, increasing r decreases the amplitude of oscillations in all three variables, while it increases the period for carbon but decreases the period for producer phosphorus. This suggests that overall r has a dampening effect on oscillations. In this case, increasing K consistently increases both amplitude and period of oscillations.

Overall, we notice that the system is more sensitive to the grazer ingestion rate than to the intrinsic growth rate of the producer. Therefore, the grazer's impact on the turnover rate is more influential in the contrast between terrestrial and aquatic ecosystems than the producer's. In general, the sensitivity rank of r decreases as light intensity increases, and the system is overall more sensitive to c for intermediate to high light levels. Given this is local sensitivity analysis, it is entirely dependent on the baseline parameters. Note that for the parameters used, the systems with K = 0.25 and K = 1.00 approach the coexistence equilibrium; K = 0.75 produces coexistence oscillations; and K = 2.00 approaches the grazer extinction boundary equilibrium. Since we explicitly know the form of the grazer extinction equilibrium, the sensitivity results for K = 2.00 are not unexpected. However, the results for K = 0.25 and K = 1.00 give us insight into the coexistence equilibrium that was not solved for explicitly.

4.3. Bifurcation analysis

Bifurcation analysis was performed using MatCont [Citation3]. For all one parameter bifurcation diagrams, a solid blue curve represents a stable equilibrium point; a magenta dashed curve is an unstable equilibrium point; and a cyan dotted curve represents a stable limit cycle.

4.3.1. One parameter

In Figure , we see that the boundary equilibria are unstable for all $r \in (0, 2]$ . The coexistence equilibrium is always stable for $K = 0.25, r \in (0, 2]$ , and c held at its baseline value. There is a Neutral Saddle Equilibrium that occurs at r = 0.0574999999999425, as well as a branch point, a Hopf point, and another branch point which occur at numbers that are e−13, e−11, and e−12 (i.e. essentially 0 and too small to continue in the two-parameter diagrams).

Figure 4. One parameter bifurcation diagrams for r with low light ( $K = 0.25$ ); c is held at its baseline value. x is on the left and y is on the right. There are no bifurcations. The coexistence equilibrium is stable throughout.

In Figure , for K = 0.25, we see that there is a transcritical bifurcation, which occurs at c = 0.594559616369951. At this bifurcation, the grazer extinction equilibrium becomes unstable and the coexistence equilibrium becomes stable. Note that before this bifurcation point, the coexistence equilibrium is not biologically feasible. In Section 3.3, we found that the determinant of the Jacobian evaluated at the grazer extinction equilibrium changed signs at c = 0.5946, since K = 0.25 falls under case 4. This validates our result.

Figure 5. One parameter bifurcation diagrams for c with low light ( $K = 0.25$ ); r is held at its baseline value. x is on the left and y is on the right. There is a transcritical bifurcation around c = 0.59. For $c < 0.59$ , we observe a stable grazer extinction equilibrium; for c>0.59, the coexistence equilibrium is stable.

As shown in Figure , the complete extinction equilibrium and the grazer extinction equilibrium are both unstable for the full range of r values for K = 0.75. However, there are several bifurcations. There is a saddle-node bifurcation at r = 1.151380, where an unstable saddle coexistence equilibrium collides with a stable node coexistence equilibrium. There is a Hopf bifurcation at r = 1.23910, where a stable limit cycle disappears, and a branch of the coexistence equilibrium becomes stable. Then there is another saddle-node bifurcation at r = 1.23951, at which this stable branch of the coexistence equilibrium collides with an unstable branch. Note that this interval of stability is too small to see in the diagram. It is possible that between r = 1.23910 and r = 1.23951 there is tristability. There is also a neutral saddle equilibrium point, and a branch point at a very low parameter value.

Figure 6. One parameter bifurcation diagrams for r with intermediate light (K= 0.75); c is held at its baseline value. x is on the left and y is on the right. There is a saddle-node bifurcation around $r = 1.15$ , a Hopf bifurcation around 1.23910, and another saddle-node bifurcation around r = 1.23951. For $r < 1.23910$ , there is a stable limit cycle; for $r > 1.15$ , there is a stable coexistence equilibrium. There may be bistability between a coexistence equilibrium and a limit cycle for $r \in (1.15, 1.23910)$ .

From Figure , for K = 0.75, the grazer equilibrium is stable until the transcritical bifurcation at c = 0.397464101826043, where the coexistence equilibrium becomes biologically feasible and stable. There are two saddle-node bifurcations: the first is at c = 0.62783 and the second is at c = 0.650870. Note that only the second is clearly visible in . Also, a stable limit cycle appears with an increasing amplitude at a Hopf bifurcation around c = 0.62800. As with the bifurcation diagram for r, there is clearly an interval of bistability between the Hopf bifurcation and the second saddle-node bifurcation, and there may be an interval of tristability around c ∈ [0.62783, 0.62800]. In Section???, we found that the determinant of the Jacobian evaluated at the grazer extinction equilibrium changed signs at c = 0.3975, since K = 0.75 is in Case 3. This agrees with the transcritical bifurcation point observed here.

Figure 7. One parameter bifurcation diagrams for c with intermediate light (K= 0.75); r is held at its baseline value. x is on the left and y is on the right. There is a transcritical bifurcation around $c = 0.40$ , saddle-node bifurcations at c = 0.62783, 0.65, and a Hopf bifurcation at c = 0.62800. For $c < 0.40$ , grazer extinction is stable; for $0.40 < c < 0.65$ , coexistence; and for $c > 0.62800$ , oscillations.

For all values of r for K = 2.00, the grazer extinction equilibrium is stable, as shown in Figure .

Figure 8. One parameter bifurcation diagrams for r with very high light (K = 2.00); c is held at its baseline value. x is on the left and y is on the right. There are no bifurcation points. The grazer extinction equilibrium is stable throughout.

For K = 2.00 and c, the grazer extinction equilibrium is stable until the coexistence equilibrium becomes biologically feasible and stable at $c = 0.892787137713355$ (Figure ), which matches the value found in Section 3.3. There is also a saddle-node bifurcation at $c = 1.302372$ , and a neutral saddle equilibrium at $c = 0.698479$ . Based on the dynamics, there should be another saddle-node bifurcation and a Hopf bifurcation between the unstable coexistence equilibrium and the stable limit cycle (as in ), but they could not be found using MatCont. Period doubling happens at the left end of the oscillations.

Figure 9. One parameter bifurcation diagrams for c with very high light (K = 2.00); r is held at its baseline value. x is on the left and y is on the right. There is a transcritical bifurcation around c = 0.89 and a saddle-node bifurcation at c = 1.30. The grazer extinction equilibrium is stable for c<0.89; then the coexistence equilibrium is stable for $0.89 < c < 1.30$ ; then there is a stable limit cycle for c>1.30.

4.3.2. Two parameter

For Figure (a), the vertical line is the transcritical bifurcation. The stable behaviour in the regions is: (1) grazer extinction equilibrium; and (2) coexistence equilibrium.

Figure 10. Two-parameter bifurcation diagrams, varying K. Green solid curves correspond to transcritical bifurcations, and magenta dashed curves to saddle-node bifurcations. For regional stable behaviour, (1) corresponds to the grazer extinction equilibrium, (2) to a coexistence equilibrium, (3) to coexistence oscillations, and (4) to coexistence oscillations.

For Figure (b), the vertical line is the transcritical bifurcation and the diagonal line corresponds to the saddle-node bifurcation along the coexistence equilibria curve. There is a cusp point where the two curves intersect. Note that from the one parameter analysis, we know there should also be a Hopf branch, but it could not be continued in two parameters using MatCont. There should also be another saddle-node curve. There is probably a region of bistability missing from this diagram, but the saddle-node bifurcation likely provides a decent approximation of the transition between a stable coexistence equilibrium and stable coexistence oscillations. The regional stable behaviour is: (1) grazer extinction equilibrium; (2) coexistence equilibrium; and (3) coexistence oscillations.

For Figure (c), we see the similar curves and regions to Figure (b). This is due to the smaller difference between $K = 0.75$ and K = 1.00 relative to the other increments in K. However, the shift in the saddle-node curve is sufficient to justify the different behaviours observed at baseline for K = 0.75 and K = 1.00.

For Figure (d), the vertical line is the transcritical bifurcation. The magenta curve is from the saddle-node bifurcation along the coexistence equilibrium curve. The point where the curves intersect is a cusp point. As in the one parameter bifurcation analysis, there is likely a missing Hopf curve and a missing saddle-node curve. However, these curves could not be located properly to be extended in the two parameter diagrams. The regions are: (1) grazer extinction equilibrium; (2) coexistence equilibrium; (3) coexistence oscillations; and (4) coexistence oscillations. Note that the lower part of the saddle-node branch does not appear as a limit point when one parameter bifurcation diagrams are created for r, using $c = 1.00$ or c = 2.00.

We assumed that all parameters other than r and c are the same between terrestrial and aquatic ecosystems. Since the baseline values for r and c (0.93 and 0.75 respectively), were for an aquatic system, then we consider roughly the lower left hand corner of the diagrams to represent terrestrial ecosystems and the opposite to represent aquatic. Thus, in low light conditions, we expect either grazer extinction or coexistence at an equilibrium for terrestrial ecosystems, and coexistence at equilibrium for aquatic (Figure (a)). For intermediate to high levels, we would expect to see a variety of possible dynamics for terrestrial systems including grazer extinction, coexistence at equilibrium, and coexistence oscillations; for aquatic, these results suggest coexistence would occur (Figure (b)–(c)). Lastly, for very high light levels, achievable only in a laboratory setting, we expect the terrestrial grazer to die out completely, and the aquatic system to exhibit some form of coexistence (Figure (d)).

We observe that from the grazer equation with Holling type II: $\frac{d y}{d t} < y * (\hat{e} c - \hat{d})$ Since $y \geq 0$ biologically, then this means that for $c < \hat{d} / \hat{e} = 0.2973$ , the grazer's population is decreasing, regardless of the size of the producer population or the light intensity. So for low c, we can only ever see extinction of the grazer, as found in above diagrams. Given that terrestrial populations persist, this would seem to suggest that the value of c must be above this threshold, even in terrestrial ecosystems. This may also suggest that the grazer's loss rate should also be lower in terrestrial populations than in aquatic.

4.4. WKL vs. LKE model

Wang et al. (2008) [Citation18] found that solutions of the WKL model are almost identical to those of the LKE model for small or large K, while they slightly differ for intermediate K. However, when K is near the homoclinic bifurcation point (K = 0.95), they are completely different.

To investigate if there are any similar differences between the two models as c and r vary, simulations were completed using ode23s for the WKL model and the LKE model for all possible combinations of K in {0.25, 0.75, 1.00, 2.00, 0.95}, r in {0.1, 1.0, 2.0}, and c in {0.1, 1.0, 2.0}.

The two differ slightly quantitatively for

K = 0.25: r = 1, c = 2; r = 2, c = 2.
K = 0.75: r = 1, c = 1; r = 1, c = 2; r = 2, c = 1; r = 2, c = 2.
K = 1.00: r = 1, c = 1; r = 1, c = 2; r = 2, c = 1; r = 2, c = 2.
K = 2.00: r = 1, c = 1; r = 1, c = 2; r = 2, c = 2.
K = 0.95: r = 1, c = 1; r = 1, c = 2; r = 2, c = 1; r = 2, c = 2.

We notice that there is never discernible differences for r = 0.1 or c = 0.1. Thus for very low intrinsic growth rate and grazer ingestion rate, there is no noticeable difference between the two models, regardless of the value of K. For none of the simulated combinations were there completely different dynamics as there were for the baseline values of r and c, and K = 0.95. This indicates that the bifurcation point in K shifts for different values of r and c.

The differences between the WKL and LKE models relate to the relaxation of the assumption that there is no free phosphorus in the media. Given the above, we can conclude that this assumption matters more for aquatic ecosystems with intermediate to high turnover rates than for terrestrial ecosystems. This is likely because at very low values of c, the phosphorus-related quality of the producer is less likely to be the controlling factor in the grazer's population since the inability to keep up with their death rate is more important.

5. Discussion

Multiple mathematical models have been developed to study the flow of nutrients and energy through a grazer–producer system. In particular, the WKL model allows for free phosphorus in the media. The impact of changing the light dependent carrying capacity of the producer on the dynamics of this system have been studied in the past. However, other parameters are also of interest.

In particular, the intrinsic growth rate of the producer (r) and the maximal ingestion rate of the grazer (c) vary between aquatic and terrestrial ecosystems, particularly for terrestrial systems with large producers and herbivores. In general, both rates are lower in terrestrial ecosystems than aquatic, resulting in a slower turnover rate for producer biomass in terrestrial-based than aquatic-based ecosystems. All other parameters are assumed to be the same, regardless of whether the system is terrestrial or aquatic.

For very low r and c, extinction of the grazer is observed numerically; for very high r and c, oscillatory coexistence or coexistence at a steady state is observed, depending on the value of K. This seems to suggest that aquatic ecosystems are more prone to exhibiting coexistence than terrestrial ecosystems.

Overall, local sensitivity analysis implies that r and c are not the most important parameters in determining the asymptotic behaviour of the system. Other parameters have more influence over the results for this particular parameter regime, but these may differ in general between aquatic and terrestrial ecosystems. Generally c has more influence than r, and changing K has more of an influence on the sensitivity of the system to r than it does on c. We also observe that the grazer loss rate is more influential for intermediate to high light levels, and that the system becomes more sensitive to the light intensity dependent carrying capacity as it increases. Note that this is likely because at (very) high light levels, the system for these parameter values reaches the grazer extinction equilibrium, where the prey population is at their light intensity dependent carrying capacity.

Part of the reason that the analysis indicates that terrestrial populations should not persist could be due to other parameters that should differ. As mentioned in the bifurcation analysis, the grazer's loss rate is likely to have a strong impact, which is reaffirmed by the fact that it was one of the parameters the system was the most sensitive to at extreme values of K and the fact that it is part of the stability condition of the grazer extinction equilibrium. Light intensity likely also differs considerably between terrestrial and aquatic ecosystems, and we see that the grazer extinction region shrinks with increasing K until a point and then grows again. Intermediate light levels in terrestrial ecosystems may explain the persistence of grazer populations observed naturally, while a lower grazer loss rate may explain terrestrial grazer persistence in low light conditions (e.g. in the shade of a dense rainforest canopy).

Future work could include examination of larger parameter ranges. During bifurcation analysis, bifurcations were observed at higher values of the parameters in some cases, but were not included due to the a priori parameter restrictions. It bears mentioning that the baseline values used for r and c, which lie in the middle of the investigated parameter ranges, are for an aquatic ecosystem. Further investigation of data to determine other parameter regimes to test would help to more definitively contrast these ecosystems. This could help with determining if there are other parameters that differ largely that may contribute to the results observed naturally. Given the system's sensitivity to grazer loss rate, further explicit consideration of this parameter in addition to those examined here may help explain the unrealistic results implying that land herbivores cannot persist. Global stability and sensitivity analyses focussing on r and c also have yet to be completed.

Most stoichiometric models, including the ones mentioned here – LKE and WKL [Citation9, Citation18] – assume strict homeostasis for heterotrophs, that is, the grazer in the system must maintain a specific, fixed nutrient ratio within its cells, regardless of the nutrient availability in its food [Citation16]. This is in contrast to the larger variation in chemical content in the organisms it consumes. However, this assumption is not completely realistic. It has been shown to be reasonable when the stoichiometric variability of heterotrophs is sufficiently narrow, independent of variation in their food source, and to be not valid for herbivores with small mortality rates [Citation20]. Therefore, in the case that explicit consideration of grazer loss rate is taken into account, a model may need to be used that relaxes or removes this assumption [Citation19].

There are several unsolved open mathematical problems for the nonsmooth WKL model [Citation18]. As mentioned previously, global stability analysis still has yet to be completed. Even within local stability analysis, there are cases where no conclusions could be reached (see Appendix). The Holling type II functional responses used here produce equilibria which depend on many parameters, as well as similarly complicated Jacobian matrices. These complexities make local stability analysis challenging. A necessary condition for stability of the extinction equilibrium may be intriguing to find, given the one presented here is only sufficient. In addition, all bifurcation analysis was completed using the bifurcation software Matcont. As such, rigorous bifurcation analysis for r and c still needs to be performed mathematically. Investigation of the higher codimension bifurcations observed in the two-parameter bifurcation diagrams could also be illuminating.

Acknowledgements

The first author would like to thank the invaluable advice of Wang Research Group, as well as their help in using Matcont. The second author would like to thank an NSERC Discovery Grant for partial support on this research. We also would like to thank both anonymous reviewers for helpful suggestions.

Disclosure statement

No potential conflict of interest was reported by the authors.

Correction Statement

This article has been republished with minor changes. These changes do not impact the academic content of the article.

Additional information

Funding

This work was supported by the Natural Sciences and Engineering Research Council of Canada [grant number RGPIN-2015-04581].

References

L. Asik and A. Peace, Dynamics of a producer-grazer model incorporating the effects of phosphorus loading on grazer's growth, Bull. Math. Biol. 81(5) (2019), pp. 1352–1368. doi: 10.1007/s11538-018-00567-9
PubMed Web of Science ®Google Scholar
N. Chitnis, J.M. Hyman, and J.M. Cushing, Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model, Bull. Math. Biol. 70(5) (2008), pp. 1272–1296. doi: 10.1007/s11538-008-9299-0
PubMed Web of Science ®Google Scholar
A. Dhooge, W. Govaerts, Y.A. Kuznetsov, H.G.E. Meijer, and B Sautois, New features of the software MatCont for bifurcation analysis of dynamical systems, Math. Computer Model. Dyn. Syst. 14(2) (2008), pp. 147–175. doi: 10.1080/13873950701742754
Web of Science ®Google Scholar
L. Edelstein-Keshet, Mathematical Models in Biology, Soc. Indust. and Appl. Math., Philadelphia, 2005.
Google Scholar
J.J. Elser and Y. Kuang, Ecological stoichiometry, Encyclopedia Theor. Ecol. 84 (2002), pp. 718–722. Doi:10.1016/B978-008045405-4.00309-8
Google Scholar
D.A. Herms and W.J. Mattson, The dilemma of plants: to grow or defend, Quarterly Rev. Biol. 67(3) (1992), pp. 283–335. doi: 10.1086/417659
Web of Science ®Google Scholar
W. Leong and J.R. Pawlik, Evidence of a resource trade-off between growth and chemical defenses among Caribbean coral reef sponges, Marine Ecol. Progress Ser. 406 (2010), pp. 71–78. doi: 10.3354/meps08541
Web of Science ®Google Scholar
X. Li, H. Wang, and Y. Kuang, Global analysis of a stoichiometric producer-grazer model with holling type functional responses, J. Math. Biol. 63(5) (2011), pp. 901–932. doi: 10.1007/s00285-010-0392-2
PubMed Web of Science ®Google Scholar
I. Loladze, Y. Kuang, and J.J. Elser, Stoichiometry in producer-grazer systems: linking energy flow with element cycling, Bull. Math. Biol. 62(6) (2000), pp. 1137–1162. doi: 10.1006/bulm.2000.0201
PubMed Web of Science ®Google Scholar
A. Peace and H. Wang, Compensatory foraging in stoichiometric producer-grazer models, Bull. Math. Biol. 81(12) (2019), pp. 4932–4950. doi: 10.1007/s11538-019-00665-2
PubMed Web of Science ®Google Scholar
A. Peace, H. Wang, and Y. Kuang, Dynamics of a producer-grazer model incorporating the effects of excess food nutrient content on grazer's growth, Bull. Math. Biol. 76(9) (2014), pp. 2175–2197. doi: 10.1007/s11538-014-0006-z
PubMed Web of Science ®Google Scholar
A. Peace, Y. Zhao, I. Loladze, J.J. Elser, and Y. Kuang, A stoichiometric producer-grazer model incorporating the effects of excess food-nutrient content on consumer dynamics, Math. Biosciences 244(2) (2013), pp. 107–115. doi: 10.1016/j.mbs.2013.04.011
PubMed Web of Science ®Google Scholar
R.A.P. Pellew, The impacts of elephant, giraffe and fire upon the Acacia tortilis woodlands of the Serengeti, African J. Ecol. 21(1) (1983), pp. 41–74. doi: 10.1111/j.1365-2028.1983.tb00311.x
Web of Science ®Google Scholar
R.A. Pellew, Food consumption and energy budgets of the giraffe. J. Appl. Ecol. 21(1) (1984), pp. 141–159. doi: 10.2307/2403043
Web of Science ®Google Scholar
J.B. Shurin, D.S. Gruner, and H. Hillebrand, All wet or dried up? Real differences between aquatic and terrestrial food webs, Proc. Royal Soc. B: Biol. Sci. 273(1582) (2005), pp. 1–9. doi: 10.1098/rspb.2005.3377
Web of Science ®Google Scholar
R.W. Sterner, J.J. Elser, Ecological Stoichiometry: the Biology of Elements From Molecules to the Biosphere, Princeton University Press, Princeton, 2002.
Google Scholar
H. Wang, K. Dunning, J.J. Elser, and Y. Kuang, Daphnia species invasion, competitive exclusion, and chaotic coexistence, DCDS-B 12 (2009), pp. 481–493. doi: 10.3934/dcdsb.2009.12.481
Web of Science ®Google Scholar
H. Wang, Y. Kuang, and I. Loladze, Dynamics of a mechanistically derived stoichiometric producer-grazer model, J. Biological Dyn. 2(3) (2008), pp. 286–296. doi: 10.1080/17513750701769881
PubMedGoogle Scholar
H. Wang, Z. Lu, and A. Raghavan, Weak dynamical threshold for the ‘strict homeostasis’ assumption in ecological stoichiometry, Ecol. Model. 384 (2018), pp. 233–240. doi: 10.1016/j.ecolmodel.2018.06.027
Web of Science ®Google Scholar
H. Wang, R.W. Sterner, and J.J. Elser, On the ‘strict homeostasis’ assumption in ecological stoichiometry, Ecological Modelling 243 (2012), pp. 81–88. doi: 10.1016/j.ecolmodel.2012.06.003
Web of Science ®Google Scholar
T. Xie, X. Yang, X. Li, and H. Wang, Complete global and bifurcation analysis of a stoichiometric predator-prey model, J. Dyn. Differ. Equ. 30(2) (2018), pp. 447–472. doi: 10.1007/s10884-016-9551-5
Web of Science ®Google Scholar
X. Yang, X. Li, H. Wang, and Y. Kuang, Stability and bifurcation in a stoichiometric producer-grazer model with knife edge, SIAM J. Applied Dyn. Sys. 15(4) (2016), pp. 2051–2077. doi: 10.1137/15M1023610
Web of Science ®Google Scholar

A. Appendices

A. Appendix 1. stability analysis

A.1. Case 1: $\bar{p} < K q$ and $\bar{p} < θ \bar{x}$ : E₁ = $(\bar{p} / q, 0, \bar{p})$

Boundary equilibrium: $(\frac{d q (\hat{a} + T) - \hat{c} T}{q (d q - \hat{c})}, 0, \frac{d q (\hat{a} + T) - \hat{c} T}{d q - \hat{c}})$ Jacobian: $A_{1} = [\begin{matrix} - r & - \frac{c \bar{p}}{a q + \bar{p}} & \frac{r}{q} \\ 0 & \frac{c \hat{e} \bar{p} q}{θ (a q + \bar{p})} - \hat{d} & 0 \\ \frac{\hat{c} (T - \bar{p})}{\hat{a} + T - \bar{p}} & - \frac{\hat{a} \hat{c} θ \bar{p}}{q (\hat{a} + T - \bar{p})^{2}} - \frac{c \bar{p} q}{a q + \bar{p}} & - \frac{\hat{a} \hat{c} \bar{p}}{q (\hat{a} + T - \bar{p})^{2}} - d \end{matrix}]$ Let the entries of matrix $A_{1}$ be labelled with 1–9 (left to right, then top to bottom). Then we have $\begin{aligned} 1. & = - r \\ 2. & = \frac{- c [d q (\hat{a} + T) - \hat{c} T]}{q [a (d q - \hat{c}) + d (\hat{a} + T)] - \hat{c} T} \\ 3. & = \frac{r}{q} \\ 4. & = 0 \\ 5. & = \frac{c \hat{e} q [d q (\hat{a} + T) - \hat{c} T]}{θ (a q (d q - \hat{c}) + d q (\hat{a} + T) - \hat{c} T)} - \hat{d} \\ 6. & = 0 \\ 7. & = d q \\ 8. & = - \frac{θ (d q (\hat{a} + T) - \hat{c} T) (d q - \hat{c})}{\hat{a} \hat{c} q} - \frac{c q [d q (\hat{a} + T) - \hat{c} T]}{\hat{a} d q + (a q + T) (d q - \hat{c})} \\ 9. & = - \frac{d^{2} q^{2} (\hat{a} + T) + \hat{c} T (\hat{c} - 2 d q)}{\hat{a} \hat{c} q} \end{aligned}$

Then $\begin{aligned} T r (A_{1}) & = 1. + 5. + 9 \\ = - r + \frac{c \hat{e} q [d q (\hat{a} + T) - \hat{c} T]}{θ (a q (d q - \hat{c}) + d q (\hat{a} + T) - \hat{c} T)} - \hat{d} - \frac{d^{2} q^{2} (\hat{a} + T) + \hat{c} T (\hat{c} - 2 d q)}{\hat{a} \hat{c} q} \\ d e t (A_{1}) & = 5. (1. * 9. - 3. * 7.) \\ = (\frac{c \hat{e} q [d q (\hat{a} + T) - \hat{c} T]}{θ (a q (d q - \hat{c}) + d q (\hat{a} + T) - \hat{c} T)} - \hat{d}) (r (\frac{d^{2} q^{2} (\hat{a} + T) + \hat{c} T (\hat{c} - 2 d q)}{\hat{a} \hat{c} q}) - r d) \\ A_{1}^{11} & = 5. * 9. \\ = (\frac{c \hat{e} q [d q (\hat{a} + T) - \hat{c} T]}{θ (a q (d q - \hat{c}) + d q (\hat{a} + T) - \hat{c} T)} - \hat{d}) (- \frac{d^{2} q^{2} (\hat{a} + T) + \hat{c} T (\hat{c} - 2 d q)}{\hat{a} \hat{c} q}) \\ A_{1}^{22} & = 1. * 9. - 3. * 7. \\ = r (\frac{d^{2} q^{2} (\hat{a} + T) + \hat{c} T (\hat{c} - 2 d q)}{\hat{a} \hat{c} q}) - r d \\ A_{1}^{33} & = 1. * 5. \\ = - r * (\frac{c \hat{e} q [d q (\hat{a} + T) - \hat{c} T]}{θ (a q (d q - \hat{c}) + d q (\hat{a} + T) - \hat{c} T)} - \hat{d}) \end{aligned}$

A.2. Appendix 2. Case 2: $\bar{p} < K q$ and $\bar{p} > θ \bar{x}$ : $E_{1}$ = $(\bar{p} / q, 0, \bar{p})$

Boundary equilibrium: $(\frac{d q (\hat{a} + T) - \hat{c} T}{q (d q - \hat{c})}, 0, \frac{d q (\hat{a} + T) - \hat{c} T}{d q - \hat{c}})$ Jacobian: $A_{2} = [\begin{matrix} - r & - \frac{c \bar{p}}{a q + \bar{p}} & \frac{r}{q} \\ 0 & \frac{c \hat{e} \bar{p}}{a q + \bar{p}} - \hat{d} & 0 \\ \frac{\hat{c} (T - \bar{p})}{\hat{a} + T - \bar{p}} & - \frac{\hat{a} \hat{c} θ \bar{p}}{q (\hat{a} + T - \bar{p})^{2}} - \frac{c \bar{p} q}{a q + \bar{p}} & - \frac{\hat{a} \hat{c} \bar{p}}{q (\hat{a} + T - \bar{p})^{2}} - d \end{matrix}]$ Let the entries of matrix $A_{1}$ be labelled with 1–9 (left to right, then top to bottom). Then we have $\begin{aligned} 1. & = - r \\ 2. & = \frac{- c [d q (\hat{a} + T) - \hat{c} T]}{q [a (d q - \hat{c}) + d (\hat{a} + T)] - \hat{c} T} \\ 3. & = \frac{r}{q} \\ 4. & = 0 \\ 5. & = \frac{c \hat{e} [d q (\hat{a} + T) - \hat{c} T]}{a q (d q - \hat{c}) + d q (\hat{a} + T) - \hat{c} T} - \hat{d} \\ 6. & = 0 \\ 7. & = d q \\ 8. & = - \frac{θ (d q (\hat{a} + T) - \hat{c} T) (d q - \hat{c})}{\hat{a} \hat{c} q} - \frac{c q [d q (\hat{a} + T) - \hat{c} T]}{\hat{a} d q + (a q + T) (d q - \hat{c})} \\ 9. & = - \frac{d^{2} q^{2} (\hat{a} + T) + \hat{c} T (\hat{c} - 2 d q)}{\hat{a} \hat{c} q} \end{aligned}$

Then $\begin{aligned} T r (A_{2}) & = 1. + 5. + 9. \\ = - r + \frac{c \hat{e} [d q (\hat{a} + T) - \hat{c} T]}{a q (d q - \hat{c}) + d q (\hat{a} + T) - \hat{c} T} - \hat{d} - \frac{d^{2} q^{2} (\hat{a} + T) + \hat{c} T (\hat{c} - 2 d q)}{\hat{a} \hat{c} q} \\ d e t (A_{2}) & = 5. (1. * 9. - 3. * 7.) \\ = (\frac{c \hat{e} [d q (\hat{a} + T) - \hat{c} T]}{a q (d q - \hat{c}) + d q (\hat{a} + T) - \hat{c} T} - \hat{d}) (r (\frac{d^{2} q^{2} (\hat{a} + T) + \hat{c} T (\hat{c} - 2 d q)}{\hat{a} \hat{c} q}) - r d) \\ A_{2}^{11} & = 5. * 9. \\ = (\frac{c \hat{e} [d q (\hat{a} + T) - \hat{c} T]}{a q (d q - \hat{c}) + d q (\hat{a} + T) - \hat{c} T} - \hat{d}) (- \frac{d^{2} q^{2} (\hat{a} + T) + \hat{c} T (\hat{c} - 2 d q)}{\hat{a} \hat{c} q}) \\ A_{2}^{22} & = 1. * 9. - 3. * 7. \\ = r (\frac{d^{2} q^{2} (\hat{a} + T) + \hat{c} T (\hat{c} - 2 d q)}{\hat{a} \hat{c} q}) - r d \\ A_{2}^{33} & = 1. * 5. \\ = - r * (\frac{c \hat{e} [d q (\hat{a} + T) - \hat{c} T]}{a q (d q - \hat{c}) + d q (\hat{a} + T) - \hat{c} T} - \hat{d}) \end{aligned}$

Contrasting stoichiometric dynamics in terrestrial and aquatic grazer–producer systems

Abstract

1. Introduction

2. Model formulation

Table 1. The parameter (P) values (V) used for simulations [Citation18].