Asymmetrical resource competition in aquatic producers: Constant cell quota versus variable cell quota

1.   Introduction

Asymmetric resource competition is an important form of competition. It describes an unfair resource allocation among individuals of a population or between populations. The reasons for this unfairness include spatial distribution of resources, population characteristics, individual differences, etc. Such competition is ubiquitous in nature. For example, Lawton and Hassell in ^[1] stated that asymmetric resource competition among insects is a more common phenomenon than symmetric resource competition with equal opportunities to compete for resources. Terrestrial plants also exhibit strong asymmetry to resources due to individual differences. Taller and bigger plants always dominate the competition and obtain more resources ^[2].

Aquatic producers are the basis of aquatic communities and influence energy flow and material cycling in aquatic ecosystems. Their growth is limited by light and nutrients ^[3,4,5,6,7]. Light comes from the water surface and changes with the seasons or day and night. Aquatic producers photosynthesize mainly by absorbing blue-violet and red light in the spectrum. Light intensity descends exponentially over water depth since it is absorbed by water and aquatic producers ^[8,9,10,11]. Nutrients in aquatic water bodies come from a variety of sources. For example, natural rainfall usually causes an inflow of nutrients at the surface; industrial or domestic wastewater enters aquatic habitats through underground pipes or tributaries. Each of the nutrient inputs may disrupt the nutrient balance in the aquatic environment. The nutrients we consider here are mainly derived from lake bottom sediments. They are transported into aquatic habitats by water exchange and turbulence ^[12,13,14]. In this situation, the different spatial direction of light and nutrient supply causes asymmetric resource competition among aquatic producers.

Jäger and Diehl in ^[15] modeled the asymmetrical resource competition between aquatic producers for light and nutrients in a shallow aquatic environment. Let $x$ be the water depth coordinate with $x = 0$ at the water surface and $x = x_p+x_b$ at the water bottom. Aquatic producers are divided into two parts: pelagic producers ($P$) and benthic producers ($B$). The former is located in the pelagic habitat ($x\in[0, x_p]$), mainly composed of various phytoplankton. The latter, including submerged macrophytes and benthic algae, lives in the benthic habitat ($x\in[x_p, x_p+x_b]$). Light ($I_p, I_b$) first passes through the pelagic habitat and then reaches the benthic habitat. Following the Lambert-Beer law ^[8], it is given by $I_p(x, P) = I_0\exp\left(-k_0x-k_pPx\right), x\in[0, x_p]$ in the pelagic habitat and $I_b(x, P, B) = I_p(x_p, P)\exp\left(-k_0(x-x_p)-k_bB(x-x_p)\right), x\in[x_p, x_p+x_b]$ in the benthic habitat. Nutrients ($U, V$) are just the opposite, passing through the benthic and pelagic habitat in turn. This process involves two nutrient exchanges. One is the nutrient exchange between the benthic habitats and sediment $((b/x_b)(V_0-V))$. The other is the nutrient exchange between the pelagic and benthic habitat ($(a/x_p)(V-U)$ and $(a/x_b)(V-U)$). This form of spatial supply of resources results in pelagic producers having the priority to use light and benthic producers having the preemptive right for nutrients. This creates unfair allocations of resources between pelagic and benthic producers. In ^[15], Jäger and Diehl also introduced a mathematical model to describe this asymmetric resource competition among aquatic producers. Their model can be simplified as

The biological significance of variables and parameters in the model (1.1) can be found in Table 1. The growth rate of aquatic producers takes the multiplication of Monod functions

This means that the resource-based growth in the model (1.1) is in the Monod form. One of its distinguishing feature is that aquatic producers' cell quota (nutrient to carbon ratio $c_p, c_b$) is constant. The first objective of the present paper is to analyze the dynamic properties of the model (1.1) theoretically and to derive rigorously the basic ecological reproductive indexes for pelagic and benthic producer invasions. These analyses can well validate and complement the numerical simulations of model (1.1) in ^[15].

In aquatic ecosystems, it has long been recognized that the cell quota in aquatic producers is not fixed but constantly changing ^[16]. This variable cell quota characterizes aquatic producer quality and affects the biodiversity of aquatic ecosystems. It has been applied to investigate some ecological mechanisms and elucidate important ecological problems. For example, producer and grazer interactions ^[6,17,18], the coexistence of three species systems ^[19,20,21], and the degradation of organic matter by bacteria ^[22,23,24]. The Droop form is most commonly used to describe the resource-based growth of a population that depends on the variable cell quota. In ^[25], Wang et al. compared the Monod (constant cell quota) and droop (variable cell quota) forms for resource-based population dynamics through theoretical analysis and experimental data. They stated that a population growth model with the two forms has similar dynamics in the closed nutrients, whereas in the case of the low nutrient uptake rate, the transient dynamics are significantly different. Therefore, it is of great significance to investigate the differences and similarities of asymmetric resource competition among aquatic producers under constant and variable cell quotas.

Motivated by the above considerations, we propose the following asymmetric resource competition model with the variable cell quota between pelagic and benthic producers:

Here $Q_p$ and $Q_b$ represent the cell quotas ($N:C$) in pelagic and benthic producers, respectively. In model (1.2), the resource-based growth of aquatic producers is the Droop form. It is expressed as $r_pu_p(Q_p)g_p(P)$ and $r_bu_b(Q_b)g_b(P, B)$, where $u_i(Q_i) = (1-Q_{\min, i}/Q_i), Q_{\min, i}\leq Q_i\leq Q_{\max, i}, i = p, b$. The biomass reduction of aquatic producers is $m_pP, m_bB$ owing to respiration, predation and death, and $(v/x_p)P$ due to sinking. The aquatic producer nutrient uptake rate is $h_p(Q_p)f_p(U)$ and $h_b(Q_b)f_b(V)$, where $h_i(Q_i) = \delta_i(Q_{\max, i}-Q_i)/(Q_{\max, i}-Q_{\min, i}), i = p, b$. The dilution rate of cell quota is $r_pu_p(Q_p)g_p(P)$ and $r_bu_b(Q_b)g_b(P, B)$. Nutrients can be released after biomass loss by aquatic producers with proportions $\theta_p, \theta_b$. Another objective of this study is to explore the dynamics of model (1.2) and to compare the differences and similarities in the dynamic properties of the two models under constant and variable cell quotas. These comparisons contain the variation pattern of the basic ecological reproductive indexes and asymmetric resource competition for varying environmental factors.

The structure of this paper is organized as follows. In Section 2, we investigate the dynamics of models (1.1) and (1.2), including the dissipation, existence and stability of equilibria. The basic ecological reproductive indexes for aquatic producer invasion are rigorously derived. The differences and similarities in the dynamical properties of the two models are illustrated. In Section 3, we explore the variation pattern of the basic ecological reproductive indexes and the results of asymmetric resource competition for varying environmental factors via sensitivity analysis and some numerical bifurcation diagrams. A brief discussion and summary are in the last section.

2.   Model dynamics

In this section, we investigate the dynamics of models (1.1) and (1.2). Considering the ecological background of (1.1) and (1.2), we will explore the solutions of (1.1) and (1.2) with the nonnegative initial values

By standard mathematical arguments, (1.2) ((1.1)) has a unique nonnegative global solution for any initial values satisfying (2.1).

2.1.   Dynamics of model (1.2)

This subsection is devoted to studying the dynamic properties of the model (1.2) containing dissipation, the existence and stability of equilibria.

Theorem 2.1. System (1.2) is dissipative.

Proof. Note that $Q_{\min, i}\leq Q_i(0)\leq Q_{\max, i}$ for $i = p, b$. From the $Q_p$ and $Q_b$ equations in (1.2), we have $Q_{\min, i}\leq Q_i(t)\leq Q_{\max, i}, i = p, b$ for all $t\geq 0$. It follows from the $P$ and $B$ equations in (1.2) that

This means that

where $A_1, A_2$ satisfy

By the $U$ and $V$ equations in (1.2), we obtain

for sufficiently large $t$. Consider the following auxiliary systems

It is obvious that (2.2) has a unique positive equilibrium $(h_1^*, h_2^*)$ satisfying

We claim that $(h_1^*, h_2^*)$ is globally asymptotically stable. The Jacobian matrix at $(h_1^*, h_2^*)$ is

This shows that two eigenvalues of $J(h_1^*, h_2^*)$ have negative real parts since

Define the Dulac function $D(h_1, h_2) = h_1^{-1}h_2^{-1}$ in ${\mathbb{R}}_+^2$. A direct calculation gives

Hence, there is no positive periodic orbit for (2.2) in ${\mathbb{R}}_+^2$. This indicates that $(h_1^*, h_2^*)$ is globally asymptotically stable. From the comparison theorem, we have

since the $U$ and $V$ equations are a cooperative system.

We now investigate the existence and stability of equilibria of model (1.2). The four possible equilibria are shown below:

$E_1^v = (0, Q_{p1}, V_0, 0, Q_{b1}, V_0)$, where

$E_2^v = (0, Q_{p2}, U_2^v, B_2^v, Q_{b2}, V_2^v)$, where $Q_{p2}, U_2^v, B_2^v, Q_{b2}, V_2^v$ solve

$E_3^v = (P_3^v, Q_{p3}, U_3^v, 0, Q_{b3}, V_3^v)$, where $P_3^v, Q_{p3}, U_3^v, Q_{b3}, V_3^v$ solve

$E_4^v = (P_4^v, Q_{p4}, U_4^v, B_4^v, Q_{b4}, V_4^v),$ where $P_4^v, Q_{p4}, U_4^v, B_4^v, Q_{b4}, V_4^v$ solve

To explore asymmetrical resource competition among aquatic producers, we define the basic ecological reproductive indexes with variable cell quota for pelagic and benthic producers by

These indexes describe the average number of new aquatic producers produced by aquatic producers during one life cycle. This means that they represent the reproductive capacity of aquatic producers.

Theorem 2.2. $E_1^v$ always exists, and it is locally asymptotically stable if $\max\{R_{0}^{p, v}, R_{0}^{b, v}\} < 1$, while $E_1^v$ is unstable if $\max\{R_{0}^{p, v}, R_{0}^{b, v}\} > 1$. Furthermore, if

then $E_1^v$ is globally asymptotically stable.

Proof. It is obvious that $E_1^v\equiv(0, Q_{p1}, V_0, 0, Q_{b1}, V_0)$. The Jacobian matrix at $E_1^v$ is

where

Note that

and if $\max\{R_{0}^{p, v}, R_{0}^{b, v}\} < 1$, then $a_{11}, a_{44} < 0$. By the Routh-Hurwitz criterion, all the eigenvalues of $J(E_1^v)$ have negative real parts. This suggests that $E_1$ is locally asymptotically stable if $\max\{R_{0}^{p, v}, R_{0}^{b, v}\} < 1$. Conversely, if $\max\{R_{0}^{p, v}, R_{0}^{b, v}\} > 1$, then $E_1$ is unstable.

From the $P$ and $B$ equations in (1.2), we have

since $Q_{\min, i}\leq Q_i(t)\leq Q_{\max, i}$, $i = p, b$ for any $t\geq0$. Then

if (2.5) holds. From the theory of asymptotical autonomous systems ^[26], the $U$ and $V$ equations in (1.2) reduce to

Following the similar arguments as those in (1.2), we obtain

Thus, the $Q_p$ and $Q_b$ equations reduce to

which imply that

Hence, $E_1^v$ is globally attractive, and then it is globally asymptotically stable.

Theorem 2.3. $E_2^v$ exists if and only if $R_{0}^{b, v} > 1$. Moreover, if $R_{1}^{p, v} < 1$, then $E_2^v$ is locally asymptotically stable, while $E_2^v$ is unstable if $R_{1}^{p, v} > 1$.

Proof. From Theorem 2 in ^[23], $E^v_2$ exists uniquely if and only if $R_{0}^{b, v} > 1$. The Jacobian matrix at $E^v_2$ is

where

It is clear that $a_{11}, a_{22}$ are the two eigenvalues of $J(E_2^v)$. The remaining four eigenvalues of $J(E_2^v)$ satisfy $\lambda^4+a_1\lambda^3+a_2\lambda^2+a_3\lambda+a_4 = 0$, where

A direct calculation shows $a_{11} < 0$ if $R_{1}^{p, v} < 1$ and

By the Routh-Hurwitz criterion, all eigenvalues of $J(E_2^v)$ have negative real parts. Therefore, $E_2^v$ is locally asymptotically stable. On the contrary, if $R_{1}^{p, v} > 1$, then $E_2^v$ is unstable.

By similar arguments as those in Theorem 2.3, the existence and stability of $E_3^v$ are ensured by the following theorem. The details are omitted here.

Theorem 2.4. $E_3^v$ exists if and only if $R_{0}^{p, v} > 1$. Moreover, if $R_{1}^{b, v} < 1$, then $E_3^v$ is locally asymptotically stable, while $E_3^v$ is unstable if $R_{1}^{b, v} > 1$.

Remark 2.5. From Theorems 2.3 and 2.4, we have the following conclusion. If $R_{0}^{j, v} > 1$ and $R_{1}^{j, v} < 1$, $j = p, b$, then model (1.2) has a bistable structure, where both $E_2^v$ and $E_3^v$ are locally asymptotically stable.

Next, we use $m_p$ as the bifurcation parameter to study the existence of $E_4^v$. Let

and denote $\hat m_p^v$ by $m_b = r_bu_b(Q_{b3}(\hat m_p))g_b(P^v_3(\hat m_p), 0)$. We consider the coexistence equilibrium $E_4^v$ bifurcating from $\Pi = \{(m_p, 0, Q_{p2}, U^v_2, B^v_2, Q_{b2}, V^v_2): m_p > 0\}$ at $m_p = \bar m_p^v$ and meeting $m_p = 0$ or $\Gamma = \{(m_p, P^v_3, Q_{p3}, U^v_3, 0, Q_{b3}, V^v_3): m_p\in(0, m_p^v)\}$ at $m_p = \hat m_p^v$.

Theorem 2.6. Assume that $R_{1}^{p, v} > 1$ and $R_{1}^{b, v} > 1$ hold. Then for each fixed $m_b\in(0, m_b^v)$, $E_4^v$ exists if $m_p\in(\max\{0, \hat m_p^v\}, \bar m_{p}^v)$.

Proof. This proof is divided into two parts. The first part is to explore the existence of $E_4^v$ near $(\hat m_p^v, 0, Q_{p2}, U^v_2, B^v_2, Q_{b2}, V^v_2)$ by the local bifurcation theory (see Theorem 1.7 in ^[27]). The second part to prove that (1.2) has at least one $E_4^v$ for $m_p\in(\max\{0, \hat m_p^v\}, \bar m_{p}^v)$ by applying the global bifurcation theory (see Theorem 3.3 and Remark 3.4 in ^[28]).

(ⅰ) Local bifurcation. Define $H: {\mathbb{R}}^7_+ \rightarrow {\mathbb{R}}^6$ as

Obviously, $H(m_p, 0, Q_{p2}, U^v_2, B^v_2, Q_{b2}, V^v_2) = 0$. Let $F: = H_{(P, Q_p, U, B, Q_b, V)}(\bar m^v_{p}, 0, Q_{p2}, U^v_2, B^v_2, Q_{b2}, V^v_2)$. It follows that

for any $(\eta_1, \eta_2, \eta_3, \eta_4, \eta_5, \eta_6)\in {\mathbb{R}}^6_+$, where

For $(\eta_1, \eta_2, \eta_3, \eta_4, \eta_5, \eta_6)\in \ker F$, one can obtain

Let $\eta_1 = 1$. Note that (2.7) is a five-dimensional homogeneous linear equation, and its coefficient determinant is not zero. Hence, there exists a unique solution ($1, \hat{\eta}_2, \hat{\eta}_3, \hat{\eta}_4, \hat{\eta}_5, \hat{\eta}_6$) satisfying (2.7). This means that $\dim\ker F = 1$ and $\ker F = \text{span}\{1, \hat{\eta}_2, \hat{\eta}_3, \hat{\eta}_4, \hat{\eta}_5, \hat{\eta}_6\}$. It is obvious that

and codim range $F = 1$. A direct calculation gives

which does not belong to $\text{range }F$.

According to the Crandall-Rabinowitz bifurcation theorem (see Theorem 1.7 in ^[27]), the smooth curve $\Upsilon = \{(\mu_{p}(s), P^v_4(s), Q_{p4}(s), U^v_4(s), B^v_4(s), Q_{b4}(s), V^v_4(s)):0 < s < \varepsilon\}$ for some $\varepsilon > 0$ near $(\bar\mu_p^v, 0, Q_{p2}, U^v_2, B^v_2, Q_{b2}, V^v_2)$ contains all positive coexistence equilibria of (1.2) with the form

(ⅱ) Global bifurcation. Let $\Lambda$ be the set of all positive coexistence equilibria of (1.2). From Theorem 3.3 and Remark 3.4 in ^[28], there exists a connected set $\Lambda^+$ in $\Lambda$ such that $\Lambda^+$ connects to $\Pi$ and contains $\Upsilon$ and its closure includes $(\bar\mu_{p}^v, 0, Q_{p2}, U^v_2, B^v_2, Q_{b2}, V^v_2)$. Furthermore, $\Lambda^+$ satisfies one of the following: 1) $\Lambda^+$ is not compact in ${\mathbb{R}}^7_+$; 2) $\Lambda^+$ meets another bifurcation point $(\tilde m_p, 0, Q_{p2}, U^v_2, B^v_2, Q_{b2}, V^v_2)$ with $\tilde m_p\not = \bar m_p^v$; 3) $\Lambda^+$ contains $(m_p, \hat P^v_4, Q_{p2}+\hat Q_{p4}, U^v_2+\hat U^v_4, B^v_2+\hat B^v_4, Q_{b2}+\hat Q_{b4}, V^v_2+\hat V^v_4)$ with $0\neq(\hat P^v_4, \hat Q_{p4}, \hat U^v_4, \hat B^v_4, \hat Q_{b4}, \hat V^v_4)\in {\mathbb{X}}$, where ${\mathbb{X}}$ is a closed complement of $\ker F = \text{span}\{1, \hat{\eta}_2, \hat{\eta}_3, \hat{\eta}_4, \hat{\eta}_5, \hat{\eta}_6\}$ in ${\mathbb{R}}^6_+$.

If 3) holds, then $\hat{P}^v_4 = 0$. It is a contradiction to $\hat{P}^v_4 > 0$. If 2) holds, then there is a coexistence equilibria sequence $\left\{\left(m^{i}_{p}, (P^v_4)^{i}, Q^{i}_{p4}, (U^{v}_4)^{i}, (B^{v}_4)^{i}, Q^{i}_{b4}, (V^{v}_4)^{i}\right)\right\}$ such that

as $i\rightarrow \infty$. It follows from the first equality in (2.3) that

Letting $i \rightarrow \infty$ gives

which indicates that $\tilde{m}_p = \bar m_p^v$. It is a contradiction, and then 2) does not hold.

According to the above discussion, 1) must hold, and then $\Lambda^+$ is not compact in ${\mathbb{R}}^7_+$. From the first equality in (2.3), one can see $0 < m_p < \bar m_p^v$ if (1.2) has positive coexistence equilibria. By Theorem 2.1, we have

for any $m_p\in(0, \bar m_p^v)$. This means that $\Lambda^+$ must meet the boundary of $(0, \bar m_p^v)\times {\mathbb{R}}^6_+$. Note that $\Lambda^+$ connects to $\Pi$ as $m_p \rightarrow\bar m_p^{v, -}$ and $\Lambda^+$ cannot meet $(m_p, 0, Q_{p1}, V_0, 0, Q_{b1}, V_0)$ for any $m_p\in(0, \bar m_p^v)$. Therefore, one of the following two alternatives must happen. The first alternative is $m_p \rightarrow0$ for some fixed $m_b\in(0, m_b^v)$, which indicates that the projection of $\Lambda^+$ on the $m_p$-axis contains the interval $(0, \bar m_p^v)$. The second alternative is that $\Lambda^+$ meets $\Gamma$ at $m_p = \hat m_p^v$. Thus, $E_4^v$ exists on $(\hat m_p^v, \bar m_p^v)$. The proof is complete.

To facilitate an understanding of the dynamics of model (1.2), we use the loss rates $m_p$ and $m_b$ as parameters to describe the attractive region of the above equilibria. From (2.6) and Theorems 2.2–2.6, we let

Figure 1 displays the districts of pelagic and benthic producers from extinction to survival in the $(m_p, m_b)$-plane. In $\Delta_1^v$, the solutions of model (1.2) converge to $E_1^v$. It indicates that both pelagic and benthic producers are extirpated. In $\Delta_2^v$, benthic producers win asymmetric resource competition while pelagic producers go extinct. In this region, $E_2^v$ attracts all solutions. Correspondingly, pelagic producers dominate aquatic ecosystems while benthic producers disappear in $\Delta_3^v$. Then the solutions of model (1.2) converge to $E_3^v$. Pelagic and benthic producers can coexist in the region $\Delta_4^v$, and $E_4^v$ is an attractor. $\Delta_2^v\cap \Delta_3^v$ is a bistable region where the solutions converge to $E_2^v$ or $E_3^v$ for different initial values.

2.2.   Dynamics of model (1.1)

We investigate the dynamic properties of model (1.1). The four possible equilibria are shown below: $E_1^c\equiv(0, V_0, 0, V_0)$, $E_2^c = (0, U_2^c, B_2^c, V_2^c)$, where $U_2^c, B_2^c, V_2^c$ solve

$E_3^c = (P_3^c, U_3^c, 0, V_3^c)$, where $P_3^c, U_3^c, V_3^c$ solve

$E_4^c = (P_4^c, U_4^c, B_4^c, V_4^c),$ where $P_4^c, U_4^c, B_4^c, V_4^c$ solve

We define the basic ecological reproductive indexes with the constant cell quota as

Let

and denote $\hat m_p^c$ as $m_b = r_bf_b(V^c_3(\hat m_p^c))g_b(P^c_3(\hat m_p^c), 0)$.

Carrying out similar arguments to those in Theorems 2.1–2.6, we obtain the following theorem. The details of the proof are omitted here.

Theorem 2.7. $\rm{{(i)}}$ System (1.1) is dissipative.

$\rm{{(ii)}}$ $E_1^c$ always exists, and it is locally asymptotically stable if $\max\{R_{0}^{p, c}, R_{0}^{b, c}\} < 1$, while $E_1^c$ is unstable if $\max\{R_{0}^{p, c}, R_{0}^{b, c}\} > 1$. Furthermore, if $m_{p} > r_pg_p(0)-v/x_p, \; m_b > r_bg_b(0, 0)$, then $E_1^c$ is globally asymptotically stable.

$\rm{{(iii)}}$ $E_2^c$ exists if and only if $R_{0}^{b, c} > 1$. Moreover, if $R_{1}^{p, c} < 1$, then $E_2^c$ is locally asymptotically stable, while $E_2^c$ is unstable if $R_{1}^{p, c} > 1$.

$\rm{{(iv)}}$ $E_3^c$ exists if and only if $R_{0}^{p, c} > 1$. Moreover, if $R_{1}^{b, c} < 1$, then $E_3^c$ is locally asymptotically stable, while $E_3^c$ is unstable if $R_{1}^{b, c} > 1$.

$\rm{{(v)}}$ Assume that $R_{1}^{p, c} > 1$ and $R_{1}^{b, c} > 1$ hold. Then for each fixed $m_b\in(0, m_b^c)$, $E_4^c$ exists if $m_p\in(\max\{0, \hat m_p^c\}, \bar m_{p}^c)$.

Remark 2.8. The above theoretical analysis results give the threshold conditions for pelagic and benthic producers to invade aquatic ecosystems, respectively, and the criterion for their coexistence. These findings explain and complement the results of the numerical analysis in ^[15]. It follows from (ⅲ) and (ⅳ) in Theorem 2.7 that both $E_2^c$ and $E_3^c$ are locally asymptotically stable if $R_{0}^{j, c} > 1$ and $R_{1}^{j, c} < 1$, $j = p, b$.

In order to compare with model (1.2), we also take $(m_p, m_b)$ as the parameters and define

From Figure 2, one can see that the solutions of model (1.1) converge to $E_i^c$ in each region $\Delta_i^c$. $\Delta_2^c\cap \Delta_3^c$ is also a bistable region of $E_2^c$ and $E_3^c$. The ecological interpretation of the corresponding region is the same as in Figure 1. It can be seen from Figures 1 and 2 that the dynamic behavior of models (1.1) and (1.2) are similar, mainly including equilibria, bistability, and no oscillation. This suggests that there is no essential difference in the dynamics of models for constant and variable cell quotas. However, the ranges of $\Delta_i^v$, $\Delta_i^c$ and the bistable region are not the same. Therefore, the results of the asymmetric competition are not identical within certain parameter ranges for two different types of cell quotas.

3.   Constant versus variable cell quotas

Models (1.1) and (1.2) have different forms of cell quota. This difference would bring some changes to the asymmetric resource competition among aquatic producers. In the following, we will compare the asymmetric resource competition and the effects of environmental factors between the constant and variable cell quotas. These comparisons contain the basic ecological reproductive indexes (2.4) and (2.8), the results of asymmetric resource competition, and the evolution trend of the biomass densities of pelagic and benthic producers with environmental factors. In Table 2, we list the ecologically reasonable parameter values applied in the numerical analysis.

In view of the model analysis in Section 2, $R_{i}^{j, v} = 1$ and $R_{i}^{j, c} = 1$, $i = 0, 1, j = p, b$ are the critical thresholds of aquatic producers from extinction to survival. From (2.4) and (2.8), one can observe that the basic ecological reproductive indexes in models (1.1) and (1.2) are not the same. A significant difference is that the indexes in (2.4) depend on variable cell quotas of pelagic and benthic producers. Figure 3 shows the changing trend of the basic ecological reproductive indexes for varying sediment nutrient concentration $V_0$ and light intensity $I_0$. The following phenomenons can be seen: 1) $R_1^{b, v} < 1$ if $V_0 > 141.5$ (Figure 3(a)); 2) $R_i^{b, c} < 1$, $i = 0, 1$ when $V_0$ tends to $0.5$ and $R_1^{b, c} < 1$ if $V_0 > 367.5$ (Figure 3(b)); 3) $R_{i}^{j, v} > 1$, $R_{i}^{b, c} > 1$, $i = 0, 1, j = p, b$ if $I_0\in(75,846)$ and $R_1^{p, v} < 1$ if $I_0 > 846$ (Figure 3(c)); 4) $R_{i}^{j, v} > 1$, $R_{i}^{b, c} > 1$, $i = 0, 1, j = p, b$ if $I_0\in(196, 1200)$. Phenomenons 1) and 2) illustrate that pelagic and benthic producers with variable cell quota are more likely to coexist in nutrient-poor environments, while in eutrophic environments, pelagic producers win asymmetric resource competition and dominate aquatic ecosystems. The reason is that variable cell quotas can well offset the adverse effects of nutrient deprivation. High nutrient input breaks the balance of resource supply, making pelagic producers dominant in asymmetric competition. Phenomenons 3) and 4) indicate that low light is detrimental for aquatic producers to coexist, while high light causes benthic producers with variable cell quota to win asymmetric resource competition. These studies suggest that the basic ecological reproductive indexes are not consistent for constant and varying cell quotas, and thus the coexistence region of pelagic and benthic producers are also very different for varying $V_0$ and $I_0$.

Based on the Latin Hypercube Sampling (LHS) and Partial Rank Correlation Coefficients (PRCCs) analyses, we explore the correlation of basic ecological reproductive indexes to model parameters. Due to the lack of available data on the distribution function, it is reasonable to assume that all input parameters conform to a normal distribution according to previous studies. The mean is the parameter value in Table 1, and the standard deviation is 0.1 times the mean. From Figure 4, one can see the correlation between $R_{i}^{j, v}$ and $R_{i}^{b, c}$, $i = 0, 1, j = p, b$ on model parameters that are linked to environmental factors. The value of the histogram represents the degree of correlation between them. High values indicate that they are more correlated, indicating that this parameter has a more significant influence on $R_{i}^{j, v}$ and $R_{i}^{b, c}$. The overall observation shows that parameters $I_0, k_0, x_p, r_b, r_p, m_b, m_p$ have a relatively large influence on the reproductive indexes, while others are not too significant for some indexes. Sensitivity analysis of basic ecological reproductive indexes also reveals differences between constant and variable cell quotas. For example, in Figure 4($c_1$), ($c_2$), the indexes $R_1^{p, v}$, $R_1^{b, v}$ show the correlations for the parameters.

We now consider the influences of environmental factors on the biomass densities of pelagic and benthic producers. These environmental factors have a necessarily close connection with resource supply, including parameters $I_0, k_0$ related to light, parameters $V_0, a, b$ related to nutrients. The following numerical bifurcation diagrams reveal the evolution trend of pelagic and benthic producer biomass densities for varying $I_0, k_0, V_0, a, b$.

Figure 5(a) shows that very low water surface light intensity is harmful to both pelagic and benthic producers. As $I_0$ increases, pelagic producers first invade aquatic ecosystems. A sharp regime shift follows, with benthic producers invading aquatic habitats and rapidly increasing biomass, while pelagic producer biomass rapidly declines. This suggests that the low light intensity is beneficial for pelagic producers, and the high light intensity allows benthic producers to win in asymmetric resource competition. For varying background light attenuation coefficient $k_0$, a sharp regime shift also occurs from benthic to pelagic producer dominance. Especially if the water is very turbid, both pelagic and benthic producers become extinct (see Figure 5(b)). For constant and variable cell quotas, pelagic and benthic producer biomass show similar evolutionary trends and regime transitions in Figure 5 (see solid and dashed lines). Two differences can also be observed. One is that when only one type of producer is present, its biomass is higher for the variable cell quota. This suggests that the variable cell quota facilitates the increase in producer biomass. The other is that the coexistence range of pelagic and benthic producers is relatively small when the cell quota is varied relative to the constant cell quota. This is because the variable cell quota reduces the dependence of pelagic producers on nutrients, making it easier to win in asymmetric resource competition.

From Figure 6, one can see the changes in pelagic and benthic producer biomass for varying sediment nutrient concentration $V_0$ and nutrient exchange rates $a, b$. Low nutrient concentrations or exchange rates allow benthic producers to win the asymmetric competition, while high ones are beneficial for pelagic producers to dominate aquatic habitats. During this process, there is a clear regime switch from benthic to pelagic producers. Similar to light-related environmental factors, pelagic and benthic producers have higher biomass and small coexistence areas for the variable cell quota when nutrient-related factors change. These findings indicate that the variable cell quota can influence asymmetric resource competition among aquatic producers and exhibit properties that differ from the constant cell quota.

4.   Discussion

Asymmetric competition is widespread in aquatic ecosystems due to the asymmetric supply of resources such as light and nutrients. Jäger and Diehl in ^[15] stated that asymmetric competition between pelagic and benthic producers might have different competition outcomes compared to classical resource competition theories based on numerical simulations. A significant difference is that pelagic and benthic producers can coexist even when one of them is at a disadvantage in terms of both light and nutrient uptake. In contrast, in the classical theory of resource competition, the conditions for the coexistence of two populations are that their utilization of resources must be significantly different. This means that asymmetric competition is more beneficial to the coexistence of pelagic and benthic producers.

In this study, we investigate the dynamic properties of the model (1.1), which explain and complement the numerical analysis results in ^[15]. Model (1.2) is proposed to describe asymmetric resource competition among aquatic producers with the variable cell quota. We also explore the dynamics of model (1.2) and compare the similarities and differences under constant and variable cell quotas. It should be emphasized that models (1.1) and (1.2) are only suitable for describing shallow aquatic environments but not all aquatic habitats.

The basic ecological reproductive indexes $R_{i}^{j, v} = 1$ and $R_{i}^{j, c} = 1$, $i = 0, 1, j = p, b$ for aquatic producer invasions are rigorously derived. If $\max\{R_{0}^{p, v}, R_{0}^{b, v}\} < 1$ ($\max\{R_{0}^{p, c}, R_{0}^{b, c}\} < 1$), then the extinction of aquatic producers is inevitable. If $R_{0}^{b, v} > 1$ and $R_{1}^{p, v} < 1$ ($R_{0}^{b, c} > 1$ and $R_{1}^{p, c} < 1$), then benthic producers win the asymmetric competition and dominate aquatic habitats. Correspondingly, if $R_{0}^{p, v} > 1$ and $R_{1}^{b, v} < 1$ ($R_{0}^{p, c} > 1$ and $R_{1}^{b, c} < 1$), then pelagic producers win the competition and dominate aquatic ecosystems. If $R_{1}^{p, v} > 1$ and $R_{1}^{b, v} > 1$ ($R_{1}^{p, c} > 1$ and $R_{1}^{b, c} > 1$), then pelagic and benthic producers can coexist. Finally, if $R_{1}^{p, v} < 1$ and $R_{1}^{b, v} < 1$ ($R_{1}^{p, c} < 1$ and $R_{1}^{b, c} < 1$), then the models have bistability, where either pelagic or benthic producers may win asymmetric competition.

The constant cell quota (Monod forms) and variable cell quota (Droop forms) have been widely used in aquatic ecological models. The former indicates that nutrient consumption and growth/cell division in aquatic producer cells occur simultaneously, while the latter indicates that the two processes are considered separately ^[25]. The existing studies show that the variable cell quota model describes the data more accurately, while the constant cell quota model is more applicable due to its simple form ^[25]. In view of the importance and wide applicability of constant and variable cell quotas in aquatic ecosystems, elucidating the similarities and differences between the two types of cell quotas can facilitate the further development of aquatic ecological models.

Here we attempt to explore the similarities and differences between constant and variable cell quota models under asymmetric resource competition. Theoretical analysis reveals the similarity of the dynamics of model (1.1) with the constant cell quota and model (1.2) with the variable cell quota. They both have four equilibria and bistable structures. However, asymmetric competition results in models (1.1) and (1.2) are not consistent for different parameter values. Sensitive analysis and bifurcation diagrams show that if there is only one aquatic producer, the aquatic producer biomass is higher when the cell quota is changed, which is beneficial to its survival. However, the variable cell quota reduces the dependence of pelagic producers on nutrients, thus enabling them to win an advantage in the asymmetric resource competition. This causes a reduction in the coexistence of pelagic and benthic producers.

Compared with the research work of Wang et al. in ^[25], our model is composed of two populations and includes the effect of light. Numerical analysis shows that light can bring some differences to some basic ecological reproductive indexes for two different types of cell quotas (see Figure 3(c), (d)). There are still some biological problems that deserve further discussion. For example, the roles of zooplankton ^[21,29] and toxins ^[30,31] in asymmetric competition.

Acknowledgments

This work is supported by NSFC-12271144, 11901140 and Heilongjiang Provincial Natural Science Foundation of China-LH2022A015.

Conflict of interest

The authors declare there is no conflict of interest.