© Wuhan University 2026

0 Introduction

Zooplankton have a significant impact on the structure and function of the lake ecosystem by acting as both predators of bacteria and phytoplankton and prey for fish. The primary energy source for zooplankton is believed to be organic carbon derived mainly from endogenous sources within the lake. However, advancements in stable isotope technology have demonstrated that carbon within aquatic food webs also originates from external inputs such as terrestrial organic carbon (TOC)^[1-2]. Recent studies indicate that external sources, including terrestrial dissolved organic carbon (T-DOC), terrestrial particulate organic carbon (T-POC), and terrestrial prey (T-prey), along with surface runoff and other pathways, significantly contribute to allochthonous carbon in lakes^[3-4]. Research has focused on how zooplankton acquire and utilize TOC^[4-6]. Findings from previous studies suggest that TOC plays a vital role in supporting zooplankton growth and reproduction^[7-9]. Additionally, studies suggest that TOC may account for 22% to 75% of zooplankton diets in nutrient-poor lakes^[10]. Zooplankton acquire these subsidies by consuming organic detritus and feeding on heterotrophic bacteria that feed on TOC^[11-13].

Due to the challenges of measuring plankton biomass directly, mathematical models are essential for understanding the physical and biological dynamics of plankton ecosystems^[14-15]. Previous research has primarily focused on zooplankton-phytoplankton systems, emphasizing dynamic behavior, equilibrium point stability, Hopf branches, global stability assessments, and global Hopf branches^[16-20]. In the context of the Three Gorges Reservoir area, a novel nutrient-plankton dynamics model was developed to explore the relationship between global stability and nutrients. The findings revealed that environmental pollution, nutrient cycling, and overfishing have the potential to affect water quality and the collapse of aquatic ecosystems^[21]. Furthermore, Zhang and Wang^[22] introduced Holling type I functional responses for phytoplankton-nutrient interactions and Holling type Ⅱ for phytoplankton-zooplankton interactions. Their analysis of the model’s global bifurcation and equilibrium point stability revealed that phytoplankton blooms can occur even at low nutrient input rates.

TOC plays a significant role in sustaining the planktonic food web in lakes and is essential for consumer nutrition. However, research on plankton systems with TOC remains limited. Inspired by reviews in Refs. [23-24], we propose a model of the carbon-phytoplankton-zooplankton (CPZ) system, which includes TOC, phytoplankton (P), and zooplankton (Z). The development of this CPZ dynamical system will enhance the understanding of planktonic ecosystem interactions. Through the analysis in this study, we aim to gain valuable insights into the overall functioning of aquatic ecosystems. Ultimately, this research seeks to contribute to a more comprehensive understanding of how TOC maintains stability in aquatic systems.

The paper is structured as follows. Section 1 formulates the CPZ model and presents two key theorems: one establishing the existence of a positive invariant set that guarantees nonnegative solutions, and another proving the existence of equilibrium points under biologically feasible parameter conditions. And then three theorems are rigorously established in Section 2 to characterize the instability of the plankton-free equilibrium and the global asymptotic stability of both the phytoplankton-free equilibrium and the coexistence equilibrium. Section 3 presents selected numerical simulations to validate the theoretical analyses and illustrate the dynamical behaviors of the CPZ model. Finally, Section 4 concludes the paper with a brief discussion.

1 CPZ Model with TOC

The CPZ dynamical model with Holling type Ⅱ functional responses for TOC-zooplankton interactions can be characterized by the following differential system:

$\{\begin{array}{l}\frac{\mathrm{d}C}{\mathrm{d}t}=\delta \left(C_{\mathrm{0}}-C\right)-\frac{a_{\mathrm{1}}C}{a_{\mathrm{2}}+C}Z,\\ \frac{\mathrm{d}P}{\mathrm{d}t}=\gamma P-\beta PZ-dP,\\ \frac{\mathrm{d}Z}{\mathrm{d}t}=k_{\mathrm{1}}\beta PZ+k_{\mathrm{2}}\frac{a_{\mathrm{1}}C}{a_{\mathrm{2}}+C}Z-\mu Z.\end{array}$(1)

For simplicity, we only consider TOC, phytoplankton, and zooplankton in the aquatic system. Let $C=C\left(t\right)$ represent the concentration of TOC available for zooplankton uptake at time $t$, while $P=P\left(t\right)$ and $Z=Z\left(t\right)$ denote the densities of phytoplankton and zooplankton, respectively. The parameter $\delta$ is the dilution rate, and $C_{\mathrm{0}}$ is the constant input concentration of TOC. Parameters $a_{\mathrm{1}}$ and $a_{\mathrm{2}}$ denote the maximum specific uptake rate and half-saturation constant of zooplankton, respectively. Additionally, $\gamma$ denotes the intrinsic growth rate of phytoplankton, whereas $d$ represents their death rate. It is reasonable to assume that $\gamma >d$. The parameter $\beta$ describes the predation coefficient of zooplankton, while parameters $k_{\mathrm{1}}$ and $k_{\mathrm{2}}$ represent the conversion rates from phytoplankton to zooplankton and from TOC to zooplankton respectively. Finally, $\mu$ indicates the death rate of zooplankton. All parameters are assumed to be nonnegative.

In a biological framework, the densities of phytoplankton and zooplankton are nonnegative for all $t>\mathrm{0}$. Similarly, the input concentration of TOC remains nonnegative for all $t>\mathrm{0}$. Consequently, system (1) is governed by the following initial conditions,

$C\left(\mathrm{0}\right)>\mathrm{0},P\left(\mathrm{0}\right)>\mathrm{0},Z\left(\mathrm{0}\right)>\mathrm{0}.$(2)

Next, we mainly focus on the positively invariant region of the CPZ dynamical system, which underpins the existence of several equilibria.

Theorem 1   The set $R_{\mathrm{3}}^{+}=\{\left(C,P,Z\right)|C,P,Z\ge \mathrm{0}\}$ is a positively invariant region for system (1).

Proof   When $C=\mathrm{0}$, the first equation of (1) reduces to $\frac{\mathrm{d}C}{\mathrm{d}t}=\delta C_{\mathrm{0}}$, it implies

$C\left(t\right)=\delta C_{\mathrm{0}}t+C\left(\mathrm{0}\right).$(3)

Integrating both sides of the last two equations of (1) with respect to $t$, we obtain the expressions for $P$ and $Z$, where

$P\left(t\right)=P\left(\mathrm{0}\right){\mathrm{e}}^{{\int }_{\mathrm{0}}^t\left(\gamma -\beta Z\left(\xi \right)-d\right)\mathrm{d}\xi },Z\left(t\right)=Z\left(\mathrm{0}\right){\mathrm{e}}^{{\int }_{\mathrm{0}}^t\left(k_{\mathrm{1}}\beta P\left(\xi \right)+k_{\mathrm{2}}\frac{a_{\mathrm{1}}C\left(\xi \right)}{a_{\mathrm{2}}+C\left(\xi \right)}-\mu \right)\mathrm{d}\xi }.$(4)

Using the initial value conditions (2)-(4), it is evident that $C,P,Z\ge \mathrm{0}$ hold for all $t\ge \mathrm{0}$. This confirms that $R_{\mathrm{3}}^{+}==\{\left(C,P,Z\right)|C,P,Z\ge \mathrm{0}\}$ is indeed the positively invariant region for system (1).

The study subsequently focuses on positive constant equilibrium solutions of system (1), which are derived through the solution of the following algebraic equation system,

$\{\begin{array}{l}\delta \left(C_{\mathrm{0}}-C\right)-\frac{a_{\mathrm{1}}C}{a_{\mathrm{2}}+C}Z=\mathrm{0},\\ \gamma P-\beta PZ-dP=\mathrm{0},\\ k_{\mathrm{1}}\beta PZ+k_{\mathrm{2}}\frac{a_{\mathrm{1}}C}{a_{\mathrm{2}}+C}Z-\mu Z=\mathrm{0}.\end{array}$(5)

The solution to the third equation of (5) yields $Z=\mathrm{0}$ or

$P=\frac{\mu }{k_{\mathrm{1}}\beta }-\frac{k_{\mathrm{2}}{a}_{\mathrm{1}}C}{k_{\mathrm{1}}\beta \left(a_{\mathrm{2}}+C\right)}\triangleq {\phi }_{\mathrm{1}}\left(C\right).$(6)

Similarly, the second equation of (5) implies in $P=\mathrm{0}$ or

$Z=\frac{\gamma -d}{\beta }.$(7)

When $Z=\mathrm{0}$, solving the first two equations of (5) yields $C=C_{\mathrm{0}}$, and $P=\mathrm{0}$.

Therefore, the plankton-free equilibrium $E_{\mathrm{0}}\left(C_{\mathrm{0}},\mathrm{0,0}\right)$ always exists. This equilibrium corresponds to a state in which only TOC persists, with no viable plankton populations.

Next, we turn to analyzing interior equilibria, which require solving the full system of the algebraic equations (5).

The third equation of (5) can be rewritten as

$\frac{a_{\mathrm{1}}C}{a_{\mathrm{2}}+C}Z=\frac{\mu }{k_{\mathrm{2}}}Z-\frac{k_{\mathrm{1}}\beta P}{k_{\mathrm{2}}}Z.$

Substituting the left-hand side expression with its right-hand side into the first equation of (5), we derive the following algebraic equation:

$\delta \left(C_{\mathrm{0}}-C\right)-\frac{\mu }{k_{\mathrm{2}}}Z+\frac{k_{\mathrm{1}}\beta P}{k_{\mathrm{2}}}Z=\mathrm{0}.$(8)

When $P=\mathrm{0}$, the solutions to equations (6) and (8) are given by

$C_{\mathrm{1}}=\frac{\mu a_{\mathrm{2}}}{k_{\mathrm{2}}{a}_{\mathrm{1}}-\mu },Z_{\mathrm{1}}=\frac{k_{\mathrm{2}}\delta \left(C_{\mathrm{0}}-C_{\mathrm{1}}\right)}{\mu }.$(9)

Consequently, when $k_{\mathrm{2}}{a}_{\mathrm{1}}>\mu$ and $C_{\mathrm{0}}>C_{\mathrm{1}}$, the phytoplankton-free equilibrium $E_{\mathrm{1}}\left(C_{\mathrm{1}},\mathrm{0},Z_{\mathrm{1}}\right)$ exists, which involves TOC and zooplankton. This equilibrium characterizes a stable state where TOC and zooplankton coexist in the absence of phytoplankton.

To investigate the coexistence equilibrium, we first substitute the zooplankton density from (7) into the carbon balance equation (8), a key step in linking phytoplankton and TOC dynamics. Subsequently, by performing algebraic manipulation (i.e., multiplying both sides by $k_{\mathrm{2}}\beta$), we obtain the following algebraic equation:

$k_{\mathrm{2}}\beta \delta \left(C_{\mathrm{0}}-C\right)-\mu \left(\gamma -d\right)+k_{\mathrm{1}}\beta \left(\gamma -d\right)P=\mathrm{0}.$(10)

This equation describes the balance between TOC input and planktonic nutrient cycling. Solving equation (10) leads to the expression

$P=\frac{\mu }{k_{\mathrm{1}}\beta }-\frac{k_{\mathrm{2}}\delta \left(C_{\mathrm{0}}-C\right)}{k_{\mathrm{1}}\left(\gamma -d\right)}\triangleq {\phi }_{\mathrm{2}}\left(C\right).$(11)

From (6) it is easy to see that $\phi {\text{'}}_{\mathrm{1}}\left(C\right)=-\frac{k_{\mathrm{2}}{a}_{\mathrm{1}}{a}_{\mathrm{2}}}{k_{\mathrm{1}}\beta {\left(a_{\mathrm{2}}+C\right)}^{\mathrm{2}}}<\mathrm{0}$. Thus, ${\phi }_{\mathrm{1}}\left(C\right)$ is a strictly decreasing and continuous function of $C$. It should also be noted that ${\phi }_{\mathrm{1}}\left(\mathrm{0}\right)>\mathrm{0},{\phi }_{\mathrm{1}}\left(C_{\mathrm{1}}\right)=\mathrm{0}$, and $\underset{C\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{\phi }_{\mathrm{1}}\left(C\right)=\frac{\mu -k_{\mathrm{2}}{a}_{\mathrm{1}}}{k_{\mathrm{1}}\beta }<\mathrm{0}$ if $k_{\mathrm{2}}{a}_{\mathrm{1}}>\mu$.

Similarly, from (11) we know that ${{\phi }^{\text{'}}}_{\mathrm{2}}\left(C\right)=\frac{k_{\mathrm{2}}\delta }{k_{\mathrm{1}}\left(\gamma -d\right)}>\mathrm{0}$. Given that $\gamma >d,$ this derivative indicates that${\phi }_{\mathrm{2}}\left(C\right)$ is a strictly increasing and continuous function of $C$. It is straightforward to see that ${\phi }_{\mathrm{2}}\left(C_{\mathrm{0}}\right)>\mathrm{0}.$

Let $\overline{C}=\frac{\mu \left(\gamma -d\right)}{k_{\mathrm{2}}\beta \delta }$. Then ${\phi }_{\mathrm{2}}\left(\mathrm{0}\right)>\mathrm{0}$ if $C_{\mathrm{0}}<\overline{C}$, and ${\phi }_{\mathrm{2}}\left(C_{\mathrm{0}}-\overline{C}\right)=\mathrm{0}$.

Thus, the equation ${\phi }_{\mathrm{1}}\left(C\right)={\phi }_{\mathrm{2}}\left(C\right)$ admits a unique positive solution if and only if $C_{\mathrm{0}}-\overline{C}<C_{\mathrm{1}}$, i.e., $C_{\mathrm{0}}<C_{\mathrm{1}}+\overline{C}$. Consequently, the unique coexistence equilibrium $E_{\mathrm{2}}\left(C_{\mathrm{2}},P_{\mathrm{2}},Z_{\mathrm{2}}\right)$ exists when $C_{\mathrm{0}}<C_{\mathrm{1}}+\overline{C}$. This equilibrium includes TOC, phytoplankton, and zooplankton, with $Z_{\mathrm{2}}=\frac{\gamma -d}{\beta }$, $P_{\mathrm{2}}={\phi }_{\mathrm{1}}\left(C_{\mathrm{2}}\right)$, and $C_{\mathrm{2}}$ being the solution to the equation ${\phi }_{\mathrm{1}}\left(C\right)={\phi }_{\mathrm{2}}\left(C\right).$

The previous analysis can be summarized as follows.

Theorem 2   For system (1), the plankton-free equilibrium $E_{\mathrm{0}}\left(C_{\mathrm{0}},\mathrm{0,0}\right)$ always exists. When $k_{\mathrm{2}}{a}_{\mathrm{1}}>\mu$, the phytoplankton-free equilibrium $E_{\mathrm{1}}\left(C_{\mathrm{1}},\mathrm{0},Z_{\mathrm{1}}\right)$ exists for $C_{\mathrm{0}}>C_{\mathrm{1}}$, where $C_{\mathrm{1}}=\frac{\mu a_{\mathrm{2}}}{k_{\mathrm{2}}{a}_{\mathrm{1}}-\mu },Z_{\mathrm{1}}=\frac{k_{\mathrm{2}}\delta \left(C_{\mathrm{0}}-C_{\mathrm{1}}\right)}{\mu }.$ Meanwhile, the unique coexistence equilibrium $E_{\mathrm{2}}\left(C_{\mathrm{2}},P_{\mathrm{2}},Z_{\mathrm{2}}\right)$ exists if and only if $C_{\mathrm{0}}<C_{\mathrm{1}}+\overline{C}$, where $C_{\mathrm{2}}$ is the unique solution to ${\phi }_{\mathrm{1}}\left(C\right)={\phi }_{\mathrm{2}}\left(C\right),P_{\mathrm{2}}={\phi }_{\mathrm{1}}\left(C_{\mathrm{2}}\right)$, and $Z_{\mathrm{2}}=\frac{\gamma -d}{\beta }$.

2 Stability of Equilibria

In this section, we will discuss the asymptotic stability of equilibria $E_{\mathrm{0}}$, $E_{\mathrm{1}}$, and $E_{\mathrm{2}}$ for system (1).

First, we investigate the local stability of equilibrium points, with a primary emphasis on the eigenvalues of the Jacobian matrices calculated at each specific equilibrium. Formally, an equilibrium point is locally asymptotically stable if every eigenvalue of its Jacobian matrix has a negative real part. In particular, the real eigenvalues must be strictly negative, and complex eigenvalues must have negative real parts.

The Jacobian matrix $\mathit{J}\left(C,P,Z\right)$ of system (1) at any point $\left(C,P,Z\right)$ in the positive quadrant $R_{\mathrm{3}}^{+}$ is given by

$\mathit{J}=\left(\begin{array}{ccc}-\delta -\frac{a_{\mathrm{1}}{a}_{\mathrm{2}}}{{\left(a_{\mathrm{2}}+C\right)}^{\mathrm{2}}}Z& \mathrm{0}& -\frac{a_{\mathrm{1}}C}{a_{\mathrm{2}}+C}\\ \mathrm{0}& \gamma -\beta Z-d& -\beta P\\ \frac{a_{\mathrm{1}}{a}_{\mathrm{2}}}{{\left(a_{\mathrm{2}}+C\right)}^{\mathrm{2}}}{k}_{\mathrm{2}}Z& k_{\mathrm{1}}\beta Z& k_{\mathrm{1}}\beta P+k_{\mathrm{2}}\frac{a_{\mathrm{1}}C}{a_{\mathrm{2}}+C}-\mu \end{array}\right).$(12)

The subsequent theorem characterizes the local asymptotic stability of equilibria for system (1).

Theorem 3   The plankton-free equilibrium $E_{\mathrm{0}}\left(C_{\mathrm{0}},\mathrm{0,0}\right)$ is unstable because the Jacobian matrix $\mathit{J}\left(E_{\mathrm{0}}\right)$ has at least one eigenvalue with a positive real part$.$ The phytoplankton-free equilibrium $E_{\mathrm{1}}\left(C_{\mathrm{1}},\mathrm{0},Z_{\mathrm{1}}\right)$ is locally asymptotically stable if and only if $C_{\mathrm{0}}>C_{\mathrm{1}}+\overline{C}$. Conversely, the unique coexistence equilibrium $E_{\mathrm{2}}\left(C_{\mathrm{2}},P_{\mathrm{2}},Z_{\mathrm{2}}\right)$ is locally asymptotically stable when $C_{\mathrm{0}}<C_{\mathrm{1}}+\overline{C}$.

Proof   The Jacobian matrix for system (1) at the plankton-free equilibrium $E_{\mathrm{0}}$ is determined by

$\mathit{J}\left(E_{\mathrm{0}}\right)=\left(\begin{array}{ccc}-\delta & \mathrm{0}& -\frac{a_{\mathrm{1}}{C}_{\mathrm{0}}}{a_{\mathrm{2}}+C_{\mathrm{0}}}\\ \mathrm{0}& \gamma -d& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& k_{\mathrm{2}}\frac{a_{\mathrm{1}}{C}_{\mathrm{0}}}{a_{\mathrm{2}}+C_{\mathrm{0}}}-\mu \end{array}\right).$(13)

The characteristic equation corresponding to $E_{\mathrm{0}}$ is derived as

$\left(\lambda +\delta \right)\left(\lambda -\left(\gamma -d\right)\right)\left(\lambda -\left(k_{\mathrm{2}}\frac{a_{\mathrm{1}}{C}_{\mathrm{0}}}{a_{\mathrm{2}}+C_{\mathrm{0}}}-\mu \right)\right)=\mathrm{0}.$(14)

Solving (14), the eigenvalues are obtained as ${\lambda }_{\mathrm{1}}=-\delta$, ${\lambda }_{\mathrm{2}}=\gamma -d$, $\mathrm{a}\mathrm{n}\mathrm{d}{\lambda }_{\mathrm{3}}=k_{\mathrm{2}}\frac{a_{\mathrm{1}}{C}_{\mathrm{0}}}{a_{\mathrm{2}}+C_{\mathrm{0}}}-\mu$.$$Since $\gamma >d$, it follows that ${\lambda }_{\mathrm{2}}=\gamma -d>\mathrm{0}$.

The plankton-free equilibrium $E_{\mathrm{0}}\left(C_{\mathrm{0}},\mathrm{0,0}\right)$ is unstable because the eigenvalue ${\lambda }_{\mathrm{2}}$ is positive.

When phytoplankton-free equilibrium $E_{\mathrm{1}}$ exists, its Jacobian matrix is

$\mathit{J}\left(E_{\mathrm{1}}\right)=\left(\begin{array}{ccc}-\delta -\frac{a_{\mathrm{1}}{a}_{\mathrm{2}}}{{\left(a_{\mathrm{2}}+C_{\mathrm{1}}\right)}^{\mathrm{2}}}{Z}_{\mathrm{1}}& \mathrm{0}& -\frac{\mu }{k_{\mathrm{2}}}\\ \mathrm{0}& \gamma -\beta Z_{\mathrm{1}}-d& \mathrm{0}\\ \frac{a_{\mathrm{1}}{a}_{\mathrm{2}}}{{\left(a_{\mathrm{2}}+C_{\mathrm{1}}\right)}^{\mathrm{2}}}{k}_{\mathrm{2}}{Z}_{\mathrm{1}}& k_{\mathrm{1}}\beta Z_{\mathrm{1}}& \mathrm{0}\end{array}\right).$(15)

One of the eigenvalues of (15) is ${\lambda }_{\mathrm{1}}=\gamma -\beta Z_{\mathrm{1}}-d$. Substituting $Z_{\mathrm{1}}=\frac{k_{\mathrm{2}}\delta \left(C_{\mathrm{0}}-C_{\mathrm{1}}\right)}{\mu }$ and $\overline{C}=\frac{\mu \left(\gamma -d\right)}{k_{\mathrm{2}}\beta \delta }$ into this expression, we have ${\lambda }_{\mathrm{1}}=\gamma -k_{\mathrm{2}}\beta \delta \frac{C_{\mathrm{0}}-C_{\mathrm{1}}}{\mu }-d=\frac{k_{\mathrm{2}}\beta \delta }{\mu }\left(\frac{\mu \left(\gamma -d\right)}{k_{\mathrm{2}}\beta \delta }-\left(C_{\mathrm{0}}-C_{\mathrm{1}}\right)\right)=\frac{k_{\mathrm{2}}\beta \delta }{\mu }\left(\overline{C}-\left(C_{\mathrm{0}}-C_{\mathrm{1}}\right)\right)$.

It follows that

${\lambda }_{\mathrm{1}}=\gamma -\beta Z_{\mathrm{1}}-d<\mathrm{0}\mathrm{h}\mathrm{o}\mathrm{l}\mathrm{d}\mathrm{s}\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{n}{C}_{\mathrm{0}}>C_{\mathrm{1}}+\overline{C}.$(16)

The other eigenvalues ${\lambda }_{\mathrm{2}}$ and ${\lambda }_{\mathrm{3}}$ are determined by the equation

${\lambda }^{\mathrm{2}}+\left(-\delta -\frac{a_{\mathrm{1}}{a}_{\mathrm{2}}}{{\left(a_{\mathrm{2}}+C_{\mathrm{1}}\right)}^{\mathrm{2}}}{Z}_{\mathrm{1}}\right)\lambda +\frac{a_{\mathrm{1}}{a}_{\mathrm{2}}\mu }{{\left(a_{\mathrm{2}}+C_{\mathrm{1}}\right)}^{\mathrm{2}}}{Z}_{\mathrm{1}}=\mathrm{0},$

where ${\lambda }_{\mathrm{2}}+{\lambda }_{\mathrm{3}}=-\delta -\frac{a_{\mathrm{1}}{a}_{\mathrm{2}}}{{\left(a_{\mathrm{2}}+C_{\mathrm{1}}\right)}^{\mathrm{2}}}{Z}_{\mathrm{1}}<\mathrm{0},{\lambda }_{\mathrm{2}}{\lambda }_{\mathrm{3}}=\frac{a_{\mathrm{1}}{a}_{\mathrm{2}}\mu }{{\left(a_{\mathrm{2}}+C_{\mathrm{1}}\right)}^{\mathrm{2}}}{Z}_{\mathrm{1}}>\mathrm{0}.$Thus, ${\lambda }_{\mathrm{2}}<\mathrm{0},\mathrm{a}\mathrm{n}\mathrm{d}{\lambda }_{\mathrm{3}}<\mathrm{0}.$

Consequently, all eigenvalues of the Jacobian matrix (15) exhibit negative real parts under the condition $C_{\mathrm{0}}>C_{\mathrm{1}}+\overline{C}$. Therefore, the phytoplankton-free equilibrium $E_{\mathrm{1}}\left(C_{\mathrm{1}},\mathrm{0},Z_{\mathrm{1}}\right)$ is locally asymptotically stable when it satisfies $C_{\mathrm{0}}>C_{\mathrm{1}}+\overline{C}$.

When coexistent equilibrium $E_{\mathrm{2}}$ exists, the Jacobian matrix at $E_{\mathrm{2}}$ is

$\mathit{J}\left(E_{\mathrm{2}}\right)=\left(\begin{array}{ccc}-\delta -\frac{a_{\mathrm{1}}{a}_{\mathrm{2}}}{{\left(a_{\mathrm{2}}+C_{\mathrm{2}}\right)}^{\mathrm{2}}}{Z}_{\mathrm{2}}& \mathrm{0}& -\frac{\mu -k_{\mathrm{1}}\beta P_{\mathrm{2}}}{k_{\mathrm{2}}}\\ \mathrm{0}& \mathrm{0}& -\beta P_{\mathrm{2}}\\ \frac{a_{\mathrm{1}}{a}_{\mathrm{2}}}{{\left(a_{\mathrm{2}}+C_{\mathrm{2}}\right)}^{\mathrm{2}}}{k}_{\mathrm{2}}{Z}_{\mathrm{2}}& k_{\mathrm{1}}\beta Z_{\mathrm{2}}& \mathrm{0}\end{array}\right).$(17)

The characteristic equation for the coexistence equilibrium $E_{\mathrm{2}}$ can be written as

${\lambda }^{\mathrm{3}}+A_{\mathrm{1}}{\lambda }^{\mathrm{2}}+A_{\mathrm{2}}\lambda +A_{\mathrm{3}}=\mathrm{0},$(18)

where $A_{\mathrm{1}}=\delta +\frac{a_{\mathrm{1}}{a}_{\mathrm{2}}}{{\left(a_{\mathrm{2}}+C_{\mathrm{2}}\right)}^{\mathrm{2}}}{Z}_{\mathrm{2}}>\mathrm{0},A_{\mathrm{2}}=k_{\mathrm{1}}{\beta }^{\mathrm{2}}{P}_{\mathrm{2}}{Z}_{\mathrm{2}}+\frac{a_{\mathrm{1}}{a}_{\mathrm{2}}}{{\left(a_{\mathrm{2}}+C_{\mathrm{2}}\right)}^{\mathrm{2}}}\left(\mu -k_{\mathrm{1}}\beta P_{\mathrm{2}}\right)Z_{\mathrm{2}}=\left[k_{\mathrm{1}}{\beta }^{\mathrm{2}}{P}_{\mathrm{2}}+\frac{a_{\mathrm{1}}{a}_{\mathrm{2}}}{{\left(a_{\mathrm{2}}+C_{\mathrm{2}}\right)}^{\mathrm{2}}}\left(\mu -k_{\mathrm{1}}\beta P_{\mathrm{2}}\right)\right]Z_{\mathrm{2}}$, and $A_{\mathrm{3}}=\left(\delta +\frac{a_{\mathrm{1}}{a}_{\mathrm{2}}}{{\left(a_{\mathrm{2}}+C_{\mathrm{2}}\right)}^{\mathrm{2}}}{Z}_{\mathrm{2}}\right)k_{\mathrm{1}}{\beta }^{\mathrm{2}}{P}_{\mathrm{2}}{Z}_{\mathrm{2}}>\mathrm{0}$.

From equation (6), we readily deduce that $\mu -k_{\mathrm{1}}\beta P_{\mathrm{2}}=\frac{k_{\mathrm{2}}{a}_{\mathrm{1}}{C}_{\mathrm{2}}}{a_{\mathrm{2}}+C_{\mathrm{2}}}$. Substituting this expression into $A_{\mathrm{2}}$, we rewrite it as $A_{\mathrm{2}}=\left[k_{\mathrm{1}}{\beta }^{\mathrm{2}}{P}_{\mathrm{2}}+\frac{a_{\mathrm{1}}^{\mathrm{2}}{a}_{\mathrm{2}}{k}_{\mathrm{2}}{C}_{\mathrm{2}}}{{\left(a_{\mathrm{2}}+C_{\mathrm{2}}\right)}^{\mathrm{3}}}\right]Z_{\mathrm{2}}$. As stated in Theorem 2, when $C_{\mathrm{0}}<C_{\mathrm{1}}+\overline{C}$, system (1) has a unique positive equilibrium $E_{\mathrm{2}}\left(C_{\mathrm{2}},P_{\mathrm{2}},Z_{\mathrm{2}}\right)$, with $C_{\mathrm{2}}>\mathrm{0},P_{\mathrm{2}}>\mathrm{0},\mathrm{a}\mathrm{n}\mathrm{d}{Z}_{\mathrm{2}}>\mathrm{0}$. All terms within the expression for $A_{\mathrm{2}}$ are positive when $C_{\mathrm{0}}<C_{\mathrm{1}}+\overline{C}$. This directly verifies that $A_{\mathrm{2}}>\mathrm{0}$ when the coexistence equilibrium exists.

These results indicate that

${\Delta }_{\mathrm{1}}=A_{\mathrm{1}}>\mathrm{0},{\Delta }_{\mathrm{2}}=A_{\mathrm{1}}{A}_{\mathrm{2}}-A_{\mathrm{3}}=k_{\mathrm{2}}{a}_{\mathrm{1}}^{\mathrm{2}}{a}_{\mathrm{2}}\left(\delta +\frac{a_{\mathrm{1}}{a}_{\mathrm{2}}{Z}_{\mathrm{2}}}{{\left(a_{\mathrm{2}}+C_{\mathrm{2}}\right)}^{\mathrm{2}}}\right)\frac{C_{\mathrm{2}}{Z}_{\mathrm{2}}}{{\left(a_{\mathrm{2}}+C_{\mathrm{2}}\right)}^{\mathrm{3}}}>\mathrm{0},\mathrm{a}\mathrm{n}\mathrm{d}{\Delta }_{\mathrm{3}}=A_{\mathrm{3}}{\Delta }_{\mathrm{2}}>\mathrm{0}.$

By the Routh-Hurwitz stability criterion^[25], it follows that all eigenvalues of equation (18) have negative real parts. Consequently, the unique coexistence equilibrium $E_{\mathrm{2}}\left(C_{\mathrm{2}},P_{\mathrm{2}},Z_{\mathrm{2}}\right)$ is locally asymptotically stable if and only if $C_{\mathrm{0}}<C_{\mathrm{1}}+\overline{C}$. This completes the proof of Theorem 3.

Next, we will establish the global stability of the equilibria $E_{\mathrm{1}}$ and $E_{\mathrm{2}}$ using the Lyapunov-LaSalle theorem^[26].

Theorem 4   The phytoplankton-free equilibrium $E_{\mathrm{1}}\left(C_{\mathrm{1}},\mathrm{0},Z_{\mathrm{1}}\right)$ is globally asymptotically stable if and only if $C_{\mathrm{0}}>C_{\mathrm{1}}+\overline{C}$.

Proof   Let $f\left(C\right)=\frac{a_{\mathrm{1}}C}{a_{\mathrm{2}}+C}$. Define a Lyapunov function as

$V(t)={\int }_{C_{\mathrm{1}}}^C\frac{f\left(\xi \right)-f\left(C_{\mathrm{1}}\right)}{f\left(\xi \right)}\mathrm{d}\xi +\frac{k_{\mathrm{1}}}{k_{\mathrm{2}}}P+\frac{\mathrm{1}}{k_{\mathrm{2}}}{\int }_{Z_{\mathrm{1}}}^Z\frac{\xi -Z_{\mathrm{1}}}{\xi }\mathrm{d}\xi .$(19)

Along the trajectories of system (1), computing the time derivative of $V(t)$ yields

$\begin{array}{l}{\frac{\mathrm{d}V}{\mathrm{d}t}|}_{(\mathrm{1})}=\frac{f\left(C\right)-f\left(C_{\mathrm{1}}\right)}{f\left(C\right)}\{\delta \left(C_{\mathrm{0}}-C\right)-f\left(C\right)Z+f\left(C\right)Z_{\mathrm{1}}-\left[\delta \left(C_{\mathrm{0}}-C_{\mathrm{1}}\right)-f\left(C_{\mathrm{1}}\right)Z_{\mathrm{1}}\right]-f\left(C\right)Z_{\mathrm{1}}\}+\frac{k_{\mathrm{1}}}{k_{\mathrm{2}}}\left(\gamma P-\beta PZ-dP+\beta P{Z}_{\mathrm{1}}-\beta P{Z}_{\mathrm{1}}\right)+\frac{\mathrm{1}}{k_{\mathrm{2}}}\frac{Z-Z_{\mathrm{1}}}{Z}\left(k_{\mathrm{1}}\beta PZ+k_{\mathrm{2}}f\left(C\right)Z-\mu Z-k_{\mathrm{2}}f\left(C_{\mathrm{1}}\right)Z+k_{\mathrm{2}}f\left(C_{\mathrm{1}}\right)Z\right)\\ =-\delta \frac{f\left(C\right)-f\left(C_{\mathrm{1}}\right)}{f\left(C\right)}\left(C-C_{\mathrm{1}}\right)-\frac{Z_{\mathrm{1}}}{f\left(C\right)}{\left(f\left(C\right)-f\left(C_{\mathrm{1}}\right)\right)}^{\mathrm{2}}-\frac{k_{\mathrm{1}}}{k_{\mathrm{2}}}P\left(\beta Z_{\mathrm{1}}-\left(\gamma -d\right)\right)\end{array}$(20)

Substituting $f\left(C\right)=\frac{a_{\mathrm{1}}C}{a_{\mathrm{2}}+C}$ and $f\left(C_{\mathrm{1}}\right)=\frac{a_{\mathrm{1}}{C}_{\mathrm{1}}}{a_{\mathrm{2}}+C_{\mathrm{1}}}$ into (20), we obtain

${\frac{\mathrm{d}V}{\mathrm{d}t}|}_{(\mathrm{1})}=-\delta \frac{a_{\mathrm{1}}{a}_{\mathrm{2}}{\left(C-C_{\mathrm{1}}\right)}^{\mathrm{2}}}{a_{\mathrm{1}}C\left(a_{\mathrm{2}}+C_{\mathrm{1}}\right)}-\frac{Z_{\mathrm{1}}}{f\left(C\right)}{\left(f\left(C\right)-f\left(C_{\mathrm{1}}\right)\right)}^{\mathrm{2}}-\frac{k_{\mathrm{1}}}{k_{\mathrm{2}}}P\left(\beta Z_{\mathrm{1}}-\left(\gamma -d\right)\right).$(21)

It follows from equation (16) that $\beta Z_{\mathrm{1}}-\left(\gamma -d\right)>\mathrm{0}$ when $C_{\mathrm{0}}>C_{\mathrm{1}}+\overline{C}$. Note that $f\left(C\right)$ is a positive increasing function, and solutions of system (1) remain nonnegative. Thus, all terms in (21) are non-positive, implying ${\frac{\mathrm{d}V}{\mathrm{d}t}|}_{(\mathrm{1})}\le \mathrm{0}$. Furthermore ${\frac{\mathrm{d}V}{\mathrm{d}t}|}_{(\mathrm{1})}=\mathrm{0}$ if and only if $C=C_{\mathrm{1}}$, $P=\mathrm{0}$, and $Z=Z_{\mathrm{1}}$. Therefore, $M=\{\left(C_{\mathrm{1}},\mathrm{0},Z_{\mathrm{1}}\right)\}$ is the largest invariant subset of system (1). By Lyapunov-LaSalle invariant principle and Theorem 3^[26], $E_{\mathrm{1}}\left(C_{\mathrm{1}},\mathrm{0},Z_{\mathrm{1}}\right)$ is globally asymptotically stable when $C_{\mathrm{0}}>C_{\mathrm{1}}+\overline{C}$. This completes the proof.

Theorem 5   The unique coexistence equilibrium $E_{\mathrm{2}}\left(C_{\mathrm{2}},P_{\mathrm{2}},Z_{\mathrm{2}}\right)$ is globally asymptotically stable if and only if $C_{\mathrm{0}}<C_{\mathrm{1}}+\overline{C}$.

Proof   We make the same assumption $f\left(C\right)=\frac{a_{\mathrm{1}}C}{a_{\mathrm{2}}+C}$, and define a new Lyapunov function

$V={\int }_{C_{\mathrm{2}}}^C\frac{f\left(\xi \right)-f\left(C_{\mathrm{2}}\right)}{f\left(\xi \right)}\mathrm{d}\xi +\frac{k_{\mathrm{1}}}{k_{\mathrm{2}}}{\int }_{P_{\mathrm{2}}}^P\frac{\xi -P_{\mathrm{2}}}{\xi }\mathrm{d}\xi +\frac{\mathrm{1}}{k_{\mathrm{2}}}{\int }_{Z_{\mathrm{2}}}^Z\frac{\xi -Z_{\mathrm{2}}}{\xi }\mathrm{d}\xi .$(22)

By calculating the time derivative of $V$ along the trajectories of system (1), we obtain

${\frac{\mathrm{d}V}{\mathrm{d}t}|}_{(\mathrm{1})}=\frac{f\left(C\right)-f\left(C_{\mathrm{2}}\right)}{f\left(C\right)}\left\{\delta \left(C_{\mathrm{0}}-C\right)-f\left(C\right)Z+f\left(C\right)Z_{\mathrm{2}}-\left[\delta \left(C_{\mathrm{0}}-C_{\mathrm{2}}\right)-f\left(C_{\mathrm{2}}\right)Z_{\mathrm{2}}\right]-f\left(C\right)Z_{\mathrm{2}}\right\}+\frac{k_{\mathrm{1}}}{k_{\mathrm{2}}}\frac{P-P_{\mathrm{2}}}{P}\left(\gamma P-\beta PZ-dP+\beta P{Z}_{\mathrm{2}}-\beta P{Z}_{\mathrm{2}}\right)+\frac{\mathrm{1}}{k_{\mathrm{2}}}\frac{Z-Z_{\mathrm{2}}}{Z}\left(k_{\mathrm{1}}\beta PZ+k_{\mathrm{2}}f\left(C\right)Z-\mu Z+k_{\mathrm{2}}f\left(C_{\mathrm{2}}\right)Z-k_{\mathrm{2}}f\left(C_{\mathrm{2}}\right)Z+k_{\mathrm{1}}\beta P_{\mathrm{2}}Z-k_{\mathrm{1}}\beta P_{\mathrm{2}}Z\right)=-\delta \frac{f\left(C\right)-f\left(C_{\mathrm{2}}\right)}{f\left(C\right)}\left(C-C_{\mathrm{2}}\right)-\frac{Z_{\mathrm{2}}}{f\left(C\right)}{\left(f\left(C\right)-f\left(C_{\mathrm{2}}\right)\right)}^{\mathrm{2}}$(23)

Following a similar approach to analyze the global asymptotic stability of $E_{\mathrm{1}}$, substituting $f\left(C\right)=\frac{a_{\mathrm{1}}C}{a_{\mathrm{2}}+C}$ and $f\left(C_{\mathrm{2}}\right)=\frac{a_{\mathrm{1}}{C}_{\mathrm{2}}}{a_{\mathrm{2}}+C_{\mathrm{2}}}$ into the Lyapunov function yields

${\frac{\mathrm{d}V}{\mathrm{d}t}|}_{(\mathrm{1})}=-\delta \frac{a_{\mathrm{1}}{a}_{\mathrm{2}}}{\left(a_{\mathrm{2}}+C\right)\left(a_{\mathrm{2}}+C_{\mathrm{2}}\right)f\left(C\right)}{\left(C-C_{\mathrm{2}}\right)}^{\mathrm{2}}-\frac{Z_{\mathrm{2}}}{f\left(C\right)}{\left(f\left(C\right)-f\left(C_{\mathrm{2}}\right)\right)}^{\mathrm{2}}.$(24)

Clearly, ${\frac{\mathrm{d}V}{\mathrm{d}t}|}_{(\mathrm{1})}\le \mathrm{0}$ holds when $C_{\mathrm{2}}>\mathrm{0},Z_{\mathrm{2}}>\mathrm{0}$, i.e., $C_{\mathrm{0}}<C_{\mathrm{1}}+\overline{C}$. Furthermore, ${\frac{\mathrm{d}V}{\mathrm{d}t}|}_{(\mathrm{1})}=\mathrm{0}$ if and only if $C=C_{\mathrm{2}}$, $P=P_{\mathrm{2}}$, and $Z=Z_{\mathrm{2}}$.

Thus, any solution of system (1) tends to $M=\{\left(C_{\mathrm{2}},P_{\mathrm{2}},Z_{\mathrm{2}}\right)\}$, which is the maximal invariant subset of system (1). By applying the Lyapunov⁃LaSalle theorem in conjunction with Theorem 3, the global asymptotic stability of $E_{\mathrm{2}}\left(C_{\mathrm{2}},P_{\mathrm{2}},Z_{\mathrm{2}}\right)$ is established. The proof is completed.

3 Simulation Results

In this section, we use numerical simulations to investigate the global dynamics of model (1) for the input concentration threshold value $C_{\mathrm{0}}$.

The values of these parameters are set as follows.

$C_{\mathrm{0}}=\mathrm{2.5},\delta =\mathrm{0.4},a_{\mathrm{1}}=\mathrm{1.35},a_{\mathrm{2}}=\mathrm{5},\gamma =\mathrm{0.48},\beta =\mathrm{0.13},d=\mathrm{0.1},\mu =\mathrm{0.2},k_{\mathrm{1}}=\mathrm{0.6},k_{\mathrm{2}}=\mathrm{0.5}.$

From the specified parameter values, it is evident that $C_{\mathrm{1}}\approx \mathrm{2.1}$, $\overline{C}\approx \mathrm{2.9}$, which satisfy the condition $C_{\mathrm{0}}<C_{\mathrm{1}}+\overline{C}$ in Theorem 5. The globally asymptotically stable coexistence equilibrium $E_{\mathrm{2}}\left(C_{\mathrm{2}},P_{\mathrm{2}},Z_{\mathrm{2}}\right)$ implies that phytoplankton and zooplankton populations will converge to a stable coexistence state under low TOC inputs. This theoretical prediction is validated by the numerical simulation results shown in Fig. 1. Figure 1(a) illustrates that TOC-phytoplankton-zooplankton dynamics exhibit limit cycle behavior. Figure 1(b) reveals that the amplitude of phytoplankton density trajectories surpasses that of TOC concentration trajectories, thereby indicating that the fluctuations in phytoplankton density are more pronounced than those in TOC concentration. Figure 1 suggests that phytoplankton play a critical role in promoting zooplankton growth and highlights zooplankton’s feeding preference for phytoplankton. By contrast, even minor TOC additions dampen oscillation amplitudes, leading to asymptotic stabilization of phytoplankton density. This suggests that lentic ecosystems are highly susceptible to trophic destabilization from low-level allochthonous subsidies.

Fig. 1 Stable behavior of ${\mathit{E}}_{\mathrm{2}}$ for ${\mathit{C}}_{\mathrm{0}}=\mathrm{2.5}$

The influx of allochthonous organic matter varies among lakes. Therefore, we explore the impact of threshold $C_{\mathrm{0}}$ variations on ecosystem dynamics in Fig. 2. Increasing $C_{\mathrm{0}}$ to 15 while maintaining other parameters constant, we obtain $C_{\mathrm{0}}>C_{\mathrm{1}}+\overline{C}$, and $k_{\mathrm{2}}{a}_{\mathrm{1}}=\mathrm{0.675}>\mu =\mathrm{0.2}$. According to Theorem 4, the equilibrium $E_{\mathrm{1}}\left(C_{\mathrm{1}},\mathrm{0},Z_{\mathrm{1}}\right)$ is globally asymptotically stable under sufficiently high input TOC concentration.

Fig. 2 Stable behavior of ${\mathit{E}}_{\mathrm{1}}$ for ${\mathit{C}}_{\mathrm{0}}=\mathrm{15}$

When the concentration of the input TOC is increased to $C_{\mathrm{0}}=\mathrm{15}$, Fig. 2(a) confirms that $E_{\mathrm{1}}$ is asymptotically stable. As depicted in Fig. 2(b), the amplitude of allochthonous organic matter concentration trajectories exceeds that of phytoplankton density trajectories. These results imply that fluctuations in allochthonous organic matter are more pronounced than those in phytoplankton. Simulation in Fig. 2 demonstrates that phytoplankton density declines to zero, signifying the extinction of this population. At the point $E_{\mathrm{1}}\left(C_{\mathrm{1}},\mathrm{0},Z_{\mathrm{1}}\right)$, when $C_{\mathrm{0}}>C_{\mathrm{1}}+\overline{C}$, the condition $\beta Z_{\mathrm{1}}>\gamma -d$ holds. This imbalance results in a negative net growth rate for phytoplankton, hence insufficient to sustain positive population growth. The predation pressure $\beta Z_{\mathrm{1}}$ exceeding the phytoplankton survival threshold triggers the extinction cascade. Moreover, Fig. 2 shows a corresponding increase in zooplankton density concurrent with phytoplankton decline, implying that zooplankton growth is influenced by allochthonous organic matter, TOC exerts inhibitory effects on phytoplankton through the term $\frac{a_{\mathrm{1}}C}{a_{\mathrm{2}}+C}Z$. Specifically, zooplankton exhibits a preference for external resources. A high allochthonous organic carbon influx may suppress phytoplankton growth, potentially causing extinction. This results from zooplankton consuming terrestrial organic carbon, which indirectly inhibits phytoplankton via trophic cascades. TOC boosts zooplankton biomass, thereby increasing predation pressure on phytoplankton through higher grazing rates. This paper emphasizes an important conclusion that the input TOC concentration influences equilibrium states and stability.

4 Conclusion

In order to explore the ecological role of TOC in lake ecosystems, a novel terrestrial organic carbon-phytoplankton-zooplankton model is proposed, where zooplankton consumption of TOC follows a Holling type Ⅱ functional response. Theoretical analysis and numerical simulations reveal that the input TOC concentration $C_{\mathrm{0}}$ is a key factor which influences the dynamics of model (1).

Theorem 4   formalizes that zooplankton can persist indefinitely even in the absence of phytoplankton, provided TOC influx exceeds a critical threshold, especially in nutrient-poor lakes. Moreover, fluctuations in TOC input concentrations have been demonstrated to govern the observed ecological dynamics. Specifically, results derived from system (1) provide robust corroboration for the hypotheses postulated in Ref. [24]. This work extends the classic phytoplankton-zooplankton modeling paradigm by demonstrating that TOC not only serves as an alternative carbon source but also restructures trophic interactions.

Notably, the model underscores the importance of terrestrial-aquatic carbon coupling for maintaining ecosystem resilience. Mechanistically, increased TOC inputs may inhibit phytoplankton blooms via resource dilution, leading to stable zooplankton densities. These findings advocate for regulating TOC input as a viable management approach to control harmful algal blooms, particularly in bloom-prone lakes.

The paper’s primary shortcoming resides in the model’s restricted validation scope, which is confined to scenarios with low filter-feeder biomass (e.g., early-stage eutrophication). New modules are imperative for high-predation environments such as shellfish farms. For intensive predation scenarios (e.g., commercial shellfish aquaculture), the future research will implement ecological complexity integration via modular design. The proposed extension involves incorporating filter-feeder functional groups (e.g., fish), thereby expanding the system into a four-variable dynamic framework.

References

All Figures

	Fig. 1 Stable behavior of ${\mathit{E}}_{\mathrm{2}}$ for ${\mathit{C}}_{\mathrm{0}}=\mathrm{2.5}$
In the text

	Fig. 2 Stable behavior of ${\mathit{E}}_{\mathrm{1}}$ for ${\mathit{C}}_{\mathrm{0}}=\mathrm{15}$
In the text