A mathematical model between keystone species: Bears, salmon and vegetation

1.   Introduction

A keystone species is one on which other species in the ecosystem largely depend. So much so that if removed there are drastic effects on the ecosystem ^[1]. Bears, salmon and vegetation are unique in an ecosystem because all three are keystone species dependent on one another ^[2]. In this ecosystem, salmon and vegetation are both food resources for bears, but the food quality differs significantly.

Salmon are much more nutrient-rich compared to vegetation, i.e., the nitrogen/phosphorous etc. concentration within salmon's bodies is much higher compared to vegetation. However, salmon must maintain homeostasis of element concentration within the body so that they must excrete excessive nutrients into the environment if the element concentration exceeds a threshold. On the other side, vegetation is less nutrient-rich but can absorb the remaining nutrient in the environment so the nutrient level within vegetation may differ significantly over time at an order of different magnitude ^[3,4,5].

Bears share common characteristics with salmon where bears are also nutrient-rich but must maintain a homeostasis of element concentration within their bodies ^[3,4,5]. When consuming poor-quality food where the element concentration is low, the energy conversion from food consumption to the bear biomass cannot be maximized. In the bear-salmon-vegetation ecosystem, salmon is good-quality food, whereas vegetation is poor-quality food. Therefore, when studying the bear-salmon-vegetation ecosystem, it is important to devise a model that depicts not only the energy transfer but also the nutrient flow between trophic levels.

One of the well-known models that study stoichiometric population dynamics is the one proposed by Loladze and Kuang ^[6]. The model builds on classical predator-prey models of one consumer and one food resource but incorporates a single limiting element phosphorus into modeling to distinguish food quantity and food quality. Following ^[6], stoichiometric models have been studied extensively, see ^{[7,8,9,10,11,12,13]} for example. More recently, Phan, Elser and Kuang extended the previous producer-grazer model framework by including multiple shared limiting elements in the modeling ^[14]. The formulated model excludes the non-smoothness in the original stoichiometric model but replicates qualitatively similar dynamics in a wide range of parameters ^[14]. We follow the modeling framework of the aforementioned studies but extend the model to include a specific nutrient-rich food resource, which agrees with the bear-salmon-vegetation ecosystem. In particular, we choose nitrogen as the limiting element of the system because nitrogen is vital for all salmon, bears and vegetation.

We organize the paper in the following. In Section 2, we formulate a stoichiometric model of one consumer and two food resources, where one food resource is of good quality but the other food resource is of poor quality. In Section 3, we analyze the existence and local stability of the steady-state solutions. In Section 4, we conduct numerical simulations that confirm the analytical results. Moreover, numerical simulations also reveal interesting dynamics when certain parameters vary in a range. We end this paper in Section 5 with conclusions and discussions.

2.   Model formulation

In the bear-salmon-vegetation ecosystem, salmon and vegetation are both food resources for bears but bears do not consume certain organs of salmon. Rather, bears leave salmon corpses around vegetation so that the vegetation absorbs any remaining nitrogen that ultimately contributes to the growth of its population. Moreover, bears and salmon are much more nutrient-rich compared to vegetation. However, salmon and bears must maintain homeostasis in nitrogen concentration within their bodies and excrete excessive nitrogen into the environment. Vegetation is less nutrient-rich but can absorb the remaining nitrogen in the environment so the nitrogen concentration in vegetation may vary at an order of different magnitude.

We first formulate a model that includes the energy transfer among trophic levels when nitrogen is abundant among all three species. Denote $S(t), V(t), B(t)$ by the salmon population density (measured in carbon per litre), the vegetation population density (measured in carbon per litre) and the bear population density (measured in carbon per litre) respectively. The model is

where $\alpha$ is the growth rate of the salmon population, $f(S)$ is the consumption rate of salmon by the bear, $\gamma$ is the intrinsic growth rate of the vegetation, $K$ is the carrying capacity of the environment, $g(V)$ represents the consumption rate of vegetation by the bear, $e_1$ and $e_2$ are the energy conversion rate from the consumption of the salmon and the vegetation respectively and $\lambda$ is the natural death rate of the bear. In (2.1), $\epsilon$ represents the rate of contribution to the vegetation growth due to the element recycling from the nearby salmon corpses left by bears. Furthermore, the bears produce excrement after digesting the salmon, which provides another source for vegetation to absorb the nitrogen, which contributes positively to vegetation growth at a rate of $\psi$ ^[2].

The above model (2.1) assumes a linear growth rate for the salmon and a logistic growth for the vegetation. Having spent most of their lives feeding and growing at sea, salmon return to spawn and die. Salmons return to their natal streams carrying marine-derived nutrients in their body tissues. In particular, adult salmon are rich in nitrogen ^[2]. Salmon do not compete for resources which implies it is safe to assume a linear growth rate of $\gamma.$ Vegetation, on the other hand, competes for resources limited by solar energy and follows a logistic growth rate. However, the vegetation population is not only limited by light but is also regulated by nitrogen concentration.

Because of the homeostasis of bears and salmon, we can safely assume that the N: C (nitrogen: carbon) ratio within bears and salmon are $\theta_1$ and $\theta_2$ respectively. It follows that the remaining nitrogen in the environment is $N-\theta_1 B(t)- \theta_2 S(t).$ Moreover, following ^[6], we assume that the N: C (nitrogen: carbon) ratio within the vegetation never falls below a level $N_V.$ Combining the assumptions, we obtain another carrying capacity of the vegetation limited by nitrogen as

Together with the carrying capacity limited by solar energy, we obtain an improved carrying capacity for the vegetation population

Following ^[6], we assume that all of the nitrogen within the system gets recycled immediately and there is no pool of free nitrogen in the environment.

It follows that the N: C ratio in vegetation is

Vegetation serves as good-quality food for bears if

or equivalently

and is poor-quality food otherwise. If vegetation is good-quality food, the biomass conversion from the vegetation to the bear follows the maximum energy intake $e_2.$ However, if vegetation is poor-quality food, the biomass conversion is reduced by a ratio

Hence, biomass conversion efficiency is not a constant but depends on both energetic and nutrient limitations. Note that we can assume that salmon is a good-quality food source for bears since salmon is nutrient-rich. Following ^[6], we formulate the biomass conversion from vegetation consumption as

Taking all the aforementioned evidence into consideration, we arrive at the model below

where the parameters are listed in Table 1.

In the following analysis, we assume a linear functional response for both predation of the salmon and the vegetation, i.e., $f(S) = \beta S$ and $g(V) = \delta V.$ We analyze a special case where the element recycling for the vegetation is minimum and can be ignored, i.e., $\epsilon = 0$ and $\psi = 0.$

In the following analysis, we restrict the positive invariant set $\Omega$ to be

where the population densities are bounded by the total element concentration $N$ in the ecosystem.

3.   Mathematical analysis

3.1.   Steady state solutions

We first analyze steady-state solutions of (2.2) when $\epsilon = 0, \psi = 0.$ The steady-state solutions are determined by

Here, (3.1) implies that $S = 0$ or $B = \alpha/\beta.$ If salmon is at extinction, i.e., $S = 0,$ then (3.3) implies that

If the bear population is at extinction, i.e., $B = 0,$ then (3.2) simplifies to

Therefore, we obtain an extinction equilibrium $E_0(0, 0, 0)$ and a vegetation-only equilibrium $E_1(0, \min\{K, N/N_V\}, 0)$ that always exist.

If

then we have

or

When (3.4) holds, (3.2) reduces to

If $K < (N-\theta_1 B)/N_V,$ then (3.6) leads to

Therefore, a boundary equilibrium $E_2(0, V^{*}, B^{*})$ exists if

If $K > (N-\theta_1 B)/N_V,$ then (3.6) simplifies to

which is equivalent to

where $\Omega_1 = \theta_1 \delta V^*, \Omega_2 = -V^*(\theta_1 \gamma +\delta N), \Omega_3 = \gamma V^* (N-N_V V^*).$ By the invariant set $\Omega,$ we have $B < N/\theta_1.$ We solve for $B$ and reject the larger root to obtain

Therefore, a boundary equilibrium $E_3(0, V^*, B^{**})$ exists if

Finally, when (3.5) holds, (3.2) reduces to

Solving for $V$ and rejecting the trivial solution leads to

Therefore, a boundary equilibrium $E_4(0, \tilde{V}, \tilde{B})$ exists if

Next, we analyze the existence of steady-state solutions when the bear population exists, i.e., $\hat{B} = \alpha/\beta.$ Substituting $\hat{B}$ into (3.3) gives

If $\delta \hat{B} V > [\delta \hat{B} (N- \theta_1 \hat{B} -\theta_2 S)]/\theta_1,$ solving for $S$ in (3.11) leads to

We further substitute $\bar{S}$ into (3.2) and obtain

By solving for $V$ and denoting the positive solution by $\bar{V},$ we obtain

Therefore, a positive equilibrium $E_5 (\bar{S}, \bar{V}, \hat{B})$ exists if

If $\delta \hat{B} V < [\delta \hat{B} (N- \theta_1 \hat{B} -\theta_2 S)]/\theta_1,$ (3.3) simplifies to

Solving for $S$ gives

We substitute (3.14) into (3.2) and obtain

If $K < (N-\theta_1 \hat{B})/N_V-[\theta_2(\lambda-e_2\delta V)]/(e_1\beta N_V),$ or equivalently $[e_1 \beta K N_V-(N-\theta_1 \hat{B})e_1 \beta+\lambda \theta_2]/[e_2 \delta \theta_2] < V,$ then (3.15) reduces to

Solving for $V$ and denoting the positive solution by $\hat{V}$ gives

Therefore, a positive equilibrium $E_6(\hat{S}, \hat{V}, \hat{B})$ exists if

If $K > (N-\theta_1 \hat{B})/N_V-[\theta_2(\lambda-e_2\delta V)]/(e_1\beta N_V),$ or equivalently $V < [e_1 \beta K N_V-(N-\theta_1 \hat{B})e_1 \beta+\lambda \theta_2]/[e_2 \delta \theta_2],$ then (3.15) reduces to

Solving for $V$ and denoting the positive solution by $V^+$ leads to

Therefore, a positive equilibrium $E_7(S^+, V^+, \hat{B})$ exists if

3.2.   Stability of the steady-state solutions

Next, we analyze the local stability of the steady-state solutions. Direct calculations lead to the Jacobian matrix

where

Then, the characteristic equation follows as $\vert \mu I-J \vert = 0.$ By substituting each equilibrium into the characteristic equation, we are able to obtain the following theorems that state the stability results. The theorem below shows the stability of the bear-extinction equilibria.

Theorem 3.1. The trivial equilibrium $E_0(0, 0, 0)$ and the vegetation-only equilibrium $E_1(0, \; \min\{K, N/N_V \}, 0)$ are unstable.

Proof. By evaluating (3.18) at the trivial equilibrium $E_0$ and solving the characteristic equation $\vert \mu I-J \vert = 0,$ we obtain the characteristic roots

which demonstrates that $E_0$ is unstable.

Similar calculations give the characteristic roots at $E_1$

which shows that $E_1$ is unstable.

The following theorem shows the stability results of the salmon-extinction equilibria $E_2, E_3, E_4.$

Theorem 3.2. The salmon-extinction equilibrium $E_2$ is locally asymptotically stable if $\alpha/\beta < B^*$ and is unstable if otherwise. The salmon-extinction equilibrium $E_3$ is locally asymptotically stable if $\alpha/\beta < B^{**} < \gamma/\delta$ and is unstable if otherwise. The salmon-extinction equilibrium $E_4$ is locally asymptotically stable if $\alpha/\beta < \tilde{B} < \gamma/\delta$ and is unstable if otherwise.

Proof. Direct calculations lead to the characteristic equation at $E_2$

where $J_{22} = \gamma-(2 \gamma V^*)/K-\delta B^*, J_{23} = -\delta V^*, J_{32} = e_2 \delta B^*.$ By using (3.7), we can simplify $J_{22} = \gamma-(2 \gamma V^*)/K-\delta B^* = \gamma -(2 \gamma V^*)/K-\gamma(1-V^*/K) = -(\gamma V^*)/K < 0.$ It follows that $\mu_1 = \alpha-\beta B^*, \mu_2+\mu_3 = -(\gamma V^*)/K < 0, \mu_2 \mu_3 = \delta V^*e_2 \delta B^* > 0.$ Therefore, $E_2$ is locally asymptotically stable if $\mu_1 < 0,$ i.e., $(\alpha/\beta) < B^*$ and is unstable if otherwise.

Similarly, by evaluating (3.18) at $E_3,$ we obtain the characteristic equation at $E_3$

where

This implies that $\mu_1 = \alpha-\beta B^{**}, \mu_2+\mu_3 = J_{22}, \mu_2 \mu_3 = -J_{23}J_{32} > 0.$ By using (3.9), we can further simplify $J_{22} = \gamma-2(\gamma -\delta B^{**})-\delta B^{**} = \delta B^{**}-\gamma.$ Therefore, $E_3$ is locally asymptotically stable if $\mu_1 < 0$ and $\mu_2+\mu_3 < 0,$ which are equivalent to $\alpha/\beta < B^{**} < \gamma/\delta.$

Finally, direct calculations give the characteristic equation at $E_4$

where

We can further simplify $J_{22}$ by using (3.10) to obtain $J_{22} = \delta \tilde{B}-\gamma.$ Moreover, $J_{33} < 0$ is satisfied automatically when $E_4$ exists. Therefore, $E_4$ is locally asymptotically stable if $\alpha/\beta < \tilde{B} < \gamma/\delta$ and is unstable otherwise.

Finally, the following theorem shows the local stabilities of the positive equilibria $E_5, E_6, E_7.$

Theorem 3.3. The positive equilibrium $E_5(\bar{S}, \bar{V}, \hat{B})$ is locally asymptotically stable if $e_1\beta \theta_1 > e_2\delta \theta_2$ and is unstable if otherwise. The positive equilibrium $E_6(\hat{S}, \hat{V}, \hat{B})$ is locally asymptotically stable if

and is unstable if otherwise. The positive equilibrium $E_7(S^+, V^+, \hat{B})$ is always unstable.

Proof. Direct calculations give the characteristic equation at $E_5$

where

By using (3.12), we can simplify $J_{33}$ as $J_{33} = -e_2 \delta \hat{B} < 0.$ Moreover, by using (3.13), we can simplify $J_{22}$ to

The above analyses show that $E_5$ is locally asymptotically stable if $J_{31} > 0$ or equivalently $e_1 \beta \theta_1 > e_2 \delta \theta_2.$

Similarly, by evaluating (3.18) at $E_6,$ we obtain the characteristic equation

where

By the Routh-Hurwitz stability criterion, $E_6$ is locally asymptotically stable if

Here we can simplify $J_{22}$ and $J_{33}$ by using (3.16) and (3.14) respectively to obtain

which further reduce (3.19) to

Therefore, $E_6$ is locally asymptotically stable if $J_{23}J_{32}+J_{13}J_{31} < 0$ or equivalently

Finally, similar to the above analyses, we obtain the characteristic equation at $E_7$

where

By the Routh-Hurwitz criterion, $E_7$ is locally asymptotically stable if

Here it is obvious that $-(J_{23}J_{32}+J_{13}J_{31}) > 0.$ Next, by substituting (3.14) and (3.17) into $J_{22},$ we can simplify

which is satisfied when $E_7$ exists. Moreover, $J_{22}(J_{23}J_{32}+J_{13}J_{31}) > J_{13}(J_{31}J_{22}-J_{21}J_{32})$ is equivalent to $J_{22}J_{23}+J_{13}J_{21} > 0.$ It remains to verify that $J_{13}(J_{31}J_{22}-J_{21}J_{32}) > 0$ or equivalently $J_{31}J_{22}-J_{21}J_{32} < 0.$ Here $J_{31}J_{22}-J_{21}J_{32} < 0$ can be simplified to $0 < -\gamma \theta_2 N_V e_2 \delta V^2$ by substituting (3.14) and (3.17), which leads to the contradiction. Therefore, $E_7$ is always unstable when exists. Thus, completes the proof.

4.   Numerical simulations

4.1.   Simulation of (2.2) with a linear predation

Now we explore the dynamics of (2.2) numerically. For the numerical simulation, we will be assuming that the salmon-bear-vegetation ecological system lies within the riparian forests of Alaska. More specifically, we will be looking at a $100\text{km}^2$ region of Lynx Creek, a tributary of the Wood River Lakes system in the Bristol Bay region of southwestern Alaska, USA $(59^{\circ}29' N, 158^{\circ}55' W)$ ^[2]. This is a well-documented region where we will be able to get accurate values for our parameters. Based on data from ^[2], we explore the following scenarios by using biologically realistic parameters.

Figure 1 shows that the salmon-extinction boundary equilibrium is locally asymptotically stable. In particular, 1(a) demonstrates that $E_2$ is locally asymptotically stable, whereas both $E_3, E_4$ do not exist under the parameter set. Figure 1(b) indicates that $E_3$ is locally asymptotically stable, whereas both $E_2, E_4$ do not exist under the parameter set. Moreover, Figure 1(c) shows that $E_4$ is locally asymptotically stable, whereas both $E_2, E_3$ do not exist under the parameter set. Via our extensive numerical experiments, we find that $E_2$ cannot coexist with either/both $E_3, E_4.$ Figure 2 demonstrates that all salmon, vegetation and bear populations coexist together. Figure 2(a) shows that the positive equilibrium $E_5$ is locally asymptotically stable, whereas $E_6$ does not exist under the parameter set. Figure 2(b) shows that the positive equilibrium $E_6$ is locally asymptotically stable, whereas $E_5$ does not exist under the parameter set.

Next, we explore how the vegetation growth rate impacts the long-term dynamics of (2.2). We generate Figure 3 by using the same set of parameter as Figure 2(b) but varying $\gamma.$ Figure 3(a) demonstrates that a small vegetation growth rate drives the vegetation population to extinction, whereas salmon coexist with bears, but the populations oscillate periodically. On the other hand, the salmon population goes to extinction, but the vegetation coexists with bears at a steady state when $\gamma$ is large. The results are reasonable biologically because a small vegetation growth rate cannot sustain the persistence of the vegetation population. However, a large vegetation growth rate facilitates the nitrogen recycling of the vegetation, which supports the bear population via consumption but drives the salmon population to extinction due to the scarcity of nitrogen.

Finally, Figure 4(a) demonstrates another scenario where the salmon, bear and vegetation coexist but the populations oscillate periodically. However, Figure 4(b) shows that a larger vegetation growth rate may stabilize the oscillating populations to a steady state.

4.2.   Simulation of (2.2) with a logistic growth of the salmon population and Holling type Ⅱ predation

Next, we explore the extended model of (2.2) where the salmon population follows a logistic growth

We consider that the predation of the salmon and the vegetation now follow the Holling type Ⅱ functional response ^[15]

where $h_1, h_2$ represent the half-saturating constants respectively.

Figure 5(a) demonstrates the bifurcation diagram of the salmon population with respect to the salmon carrying capacity $K_1.$ Moreover, Figure 5(b) demonstrates the bifurcation diagram of the vegetation population with respect to $K_1.$ The bifurcation diagram of the bear population with respect to $K_1$ is similar to Figure 5(a) and is thus omitted. Figure 5 shows that the salmon, vegetation and bear populations coexist at a steady state if the salmon carrying capacity is relatively small. However, when $K_1$ increases and passes $4.1,$ the vegetation population goes to extinction, whereas the salmon population coexists with the bear population at a steady state. Moreover, when $K_1$ further increases and passes $5.7,$ the vegetation population remains at extinction, but the salmon and bear populations oscillate periodically. The results are biologically reasonable because if the salmon carrying capacity is relatively small, the salmon, vegetation and bear compete for the nutrients but coexist due to the scarcity of the salmon population. However, if the salmon carrying capacity is at the intermediate range, the competition for nutrients drives the vegetation population to extinction because the salmon and bear are more nutrient-rich. Finally, a large salmon carrying capacity destabilizes the coexistence steady-state of the salmon and bear due to the enrichment of salmon.

Next, to compare the dynamics of the salmon-vegetation-bear model with the grazer-consumer model in ^[6], we run the simulation by using the same set of parameters as in ^[6], except the salmon carrying capacity $K_1$ and the half-saturating constant of the salmon consumption $h_1.$ Figure 6 shows the bifurcation diagram of the salmon population/vegetation population/bear population with respect to the vegetation carrying capacity $K$ respectively. If $K$ is relatively small, Figure 6 shows that the salmon, bear and vegetation populations coexist at a steady state. However, when $K$ increases and passes $0.24,$ the salmon population goes to extinction, whereas the vegetation and bear populations coexist at a steady state. If $K$ further increases and passes $0.55,$ the vegetation and bear coexist but the populations oscillate periodically. If $K$ further increases and passes $0.97,$ the vegetation and bear populations coexist but return to the steady state status. Finally, if $K$ becomes relatively large and passes $1.5,$ the salmon population is no longer extinct but coexists with the vegetation and bear populations at a steady state.

Overall, when $K$ is within the intermediate range, i.e., the salmon population is at extinction, Figures 6(a), 6(b) and 6(c) demonstrate similar dynamics with the producer-grazer model in ^[6]. The vegetation-bear coexistent steady state loses stability due to energy enrichment. However, a larger carrying capacity limited by the light energy drives the oscillating vegetation and bear populations to a steady state because of the competition for nutrients. Different from the results in ^[6], the bear population does not go to extinction but coexists with either/both salmon or/and vegetation populations when $K$ is relatively small or large. The results suggest that a low or high carrying capacity limited by light energy facilitates the persistence of all salmon, vegetation and bear populations. Moreover, an intermediate carrying capacity limited by light energy drives the salmon population to extinction.

5.   Conclusions and discussion

In this paper, we study an ecosystem of bears, salmon and vegetation, where bears consume salmon and vegetation for survival. Because salmon return to their natal streams to spawn and carry marine-derived nutrients, it is important to characterize the nutritional level difference between salmon, bears, and vegetation. In general, salmon and bears are more nutrient-rich compared to vegetation but must maintain homeostasis of the nutritional level within their bodies. On the other side, vegetation is less nutrient-rich but can recycle the remaining nutrient in the environment and the nutritional level within the body may differ significantly.

We propose a stoichiometric predator-prey model that tracks both the energy flow and the nutrient recycling from one trophic level to another. Analytical results show that boundary equilibria $E_0, E_1$ where bears are extinct exist but are always unstable. Moreover, boundary equilibria $E_2, E_3, E_4$ where salmon are extinct but bears persist may exist and remain locally asymptotically stable if certain conditions are satisfied. Positive equilibria $E_5, E_6, E_7$ where salmon, bears and vegetation coexist may exist if certain conditions are satisfied. Analyses show that $E_5$ and $E_6$ may remain locally asymptotically stable under certain conditions, but $E_7$ is always unstable.

Numerical simulations demonstrate that a small vegetation growth rate may drive the vegetation population to extinction where the salmon population and the bear population coexist in the periodic setting. Moreover, a large vegetation growth rate may drive the salmon population to extinction but the vegetation coexists with the bears at a steady state. Alternatively, the salmon, bears and vegetation populations may coexist periodically. In this scenario, a large vegetation growth rate may stabilize the system and drive the salmon, bear and vegetation to coexist at a steady state.

In this paper, for the analytical analysis, we adopt the linear functional response for the predation of the salmon and the predation of the vegetation in (2.2) for simplicity. However, a linear functional response has its limitation and is more suitable for an ecosystem with a sparse population density. A Holling type Ⅱ functional response has a saturating effect when population density becomes large and therefore may be suitable for a broader regime ^[15]. Figures 5 and 6 in the simulation also confirm that rich dynamics may occur if the Holling type Ⅱ functional response is adopted. Because of the nonlinearity of the Holling type Ⅱ functional response, we expect that the analytical analysis becomes more challenging but on the other hand may deepen our understanding of the ecosystem of keystone species, which leaves as future works.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

CM is grateful for funding from the NSERC USRA. XW is grateful for funding from the NSERC of Canada (RGPIN-2020-06825 and DGECR-2020-00369).

Conflict of interest

The authors declare there is no conflict of interest.