Dynamics of a stoichiometric phytoplankton-zooplankton model with season-driven light intensity

1.   Introduction

Phytoplankton stands as a pivotal group of primary producers within aquatic ecosystems. Through photosynthesis, phytoplankton generate ample oxygen for marine organisms, while assuming a vital role in energy flow, material cycling, and information dissemination within aquatic environments ^[1,2]. Nonetheless, over the past few decades, incidents of phytoplankton blooms have become increasingly prevalent globally ^[3,4,5,6]. Although only $2\%$ of phytoplankton species are considered harmful ^[7], these blooms disrupt the equilibrium of aquatic ecosystems and pose risks to human health. Concurrently, the mechanisms governing phytoplankton blooms remain largely elusive, and the anticipation and mitigation of such occurrences have posed challenges for contemporary researchers.

Researchers have convincingly demonstrated that phytoplankton blooms result from intricate interactions among various biotic factors ^[8,9,10] and abiotic factors ^{[11,12,13,14]}. Among these factors are nutrients, pivotal components for the sustenance of phytoplankton growth. Many nutrients, encompassing nitrogen and phosphorus, find their way into water bodies through diverse pathways, elevating nutrient levels or inducing eutrophication within aquatic ecosystems ^[11]. Accordingly, aquatic ecosystem eutrophication emerges as a pivotal trigger for phytoplankton blooms. However, it is noteworthy that phytoplankton blooms are not solely confined to eutrophic aquatic environments; occurrences have been documented even in nutrient-poor lake ecosystems ^{[15,16,17,18]} and marine ecosystems ^[19,20,21]. Consequently, the intricate relationship between nutrient concentrations and phytoplankton growth remains an ongoing investigation and discourse area.

Seasonal fluctuations exert a pervasive influence on natural ecosystems. These cyclic variations, omnipresent in nature, impact diverse abiotic facets, such as temperature and light intensity, thereby intricately intertwining with the dynamics of organisms. Empirical studies have underscored the pronounced impact of seasonal oscillations on phytoplankton proliferation. For instance, Thackeray et al. ^[22] investigated the ramifications of long-term climatic alterations on phytoplankton growth through a blend of field experimentation and statistical modeling techniques. Their findings unveiled a profound correlation between the periodic fluctuations in phytoplankton density and the corresponding seasonal periodic variation. Wang et al. ^[23] delved into the sway of climatic seasonality on the proliferation of harmful phytoplankton, elucidating that the rhythmic pulse of seasonal changes holds a pivotal role in shaping the dynamics of detrimental phytoplankton. Huppert ^[24] studied the impact of seasonal periodicity in climate on phytoplankton growth using simple non-autonomous dynamic models. Chen et al. ^[25] examined the effects of temperature and seasonal periodicity by developing a nutrient-phytoplankton dynamic model, which showed good agreement with the measured data from Lake Taihu. Their study indicated that seasonal temperature variations drive periodic oscillations in phytoplankton density, with the intensity of phytoplankton blooms strongly correlated with the seasonal mean temperature, peak temperature variations, and the duration of high temperatures. Researchers looking at the impact of seasonal periodicity on phytoplankton growth mainly focus on the effects of seasonal temperature variations on phytoplankton growth. Notably, phytoplankton utilize solar energy to drive photosynthesis and generate oxygen, underscoring the critical role of light in their life cycle. Given the sway of seasonal cycles, variations in light intensity follow suit. Consequently, it becomes imperative to investigate the repercussions of season-driven light intensity on the growth patterns of phytoplankton.

In aquatic ecosystems, zooplankton play a significant role as mediators of material cycles and conduits for energy transfer, bridging primary producers to higher trophic levels ^{[26,27,28,29,30,31]}. Their influence extends to impacting the geochemical cycles of marine organisms through direct and indirect interactions ^[32,33,34]. This centrality of zooplankton in aquatic systems has spurred significant scientific interest. For instance, Luo ^[31] scrutinized the interplay between seasonally-driven temperature variations and zooplankton density in shaping phytoplankton growth. The resultant theoretical insights showcased the influential regulatory role of zooplankton in curbing excessive phytoplankton proliferation, thereby thwarting the onset of red tide phenomena. Elser and Urabe ^[35] discussed how Daphnia controls the biomass of Chlorophyta through selective grazing in freshwater ecosystems, which significantly influences nutrient cycling and the structure of phytoplankton communities. Irigoien et al. ^[36] investigated how Copepod grazing affects Diatom populations and influences the overall structure of phytoplankton communities in marine ecosystems. They highlighted that Copepod predation could create a 'loophole' that allows certain phytoplankton species to bloom, significantly impacting the dynamics and composition of marine ecosystems. Therefore, an inclusive understanding of phytoplankton population dynamics necessitates incorporating zooplankton-driven top-down control mechanisms.

In the real world, organisms require various nutrients to maintain their growth and reproduction ^[37], among which the most three common nutrients are carbon, nitrogen, and phosphorus ^[38]. Nevertheless, it has long been known that the between-species composition of chemical elements is different; moreover, the heterogeneity of chemical elements affects the mass transfer between trophic levels ^[39]. The scientific study of the balance among multiple chemical elements in ecological interactions is referred to as stoichiometry, and it can be applied to explore energy balance, nutrient cycling, and nutrient restrictions. One of the key advantages of using stoichiometry in ecological models is its ability to provide a more realistic representation of nutrient limitations and elemental imbalances, which are crucial for accurately predicting population dynamics and ecosystem processes. Loladze et al. ^[40] used the form of a minimum function to characterize the effect of prey quality on the food production rate of predators and improved the traditional Rosenzweig-McArthur type model to present the LKE model, which reveals that prey quality also plays a crucial role in the density variation of predators. In order to reveal the mechanism of nutrient content on the growth of phytoplankton and zooplankton in aquatic ecosystems, it is realistic to incorporate chemical heterogeneity into the model.

Inspired by the preceding discourse, we endeavor to construct a stoichiometric phytoplankton-zooplankton model that encompasses the influence of general season-driven light intensity. Our fundamental objectives are twofold: 1) To scrutinize the intricate interaction dynamics between phytoplankton and zooplankton in the context of stoichiometric limitations and time-varying parameters. 2) To establish sufficient conditions for phytoplankton proliferation and subsequently propose strategic interventions to mitigate the occurrence of phytoplankton outbreaks.

This article is structured as follows. In Section 2, we introduce the development of a stoichiometric phytoplankton-zooplankton model augmented with the incorporation of season-driven light intensity. In Section 3, the dynamics of the model are delineated, encompassing aspects such as positive invariance, dissipativity, boundary dynamics, and internal dynamics. In Section 4, we showcase the application of numerical simulations to vividly validate our proposed theoretical findings. Furthermore, we explore the impacts of seasonal light intensity, nutrient availability, and zooplankton on phytoplankton growth in this section. Finally, our conclusions are succinctly encapsulated in Section 5.

2.   Model formulation

In this study, we focus on elucidating the influence of two pivotal constituents, carbon (C) and phosphorus (P), while presupposing the copious abundance of all other elements within the ambient environment. The proportional relationship between these fundamental chemical constituents, particularly the phosphorus-to-carbon ratio, stands as an illustrative gauge of producer quality. Furthermore, the assimilation of carbon by the producer is circumscribed by its photosynthetic prowess. In conditions of ample luminosity, it is anticipated that carbon (C) availability constraints are less pronounced compared to phosphorus (P).

In ^[40], Loladze et al. investigated the following stoichiometric LKE model

Here, $x$ and $y$ are the carbon concentration of phytoplankton and zooplankton respectively, $r$ represents the maximal growth rate of phytoplankton, $K$ is the constant carrying capacity of phytoplankton that depends on light intensity, $f(x)$ is the zooplankton's ingestion rate, which is a bounded smooth function, $\hat{e}$ describes the production efficiency in carbon terms for zooplankton, $d$ presents the loss rate of zooplankton ^[41], $P$ is the total phosphorus content, and $\theta$ stands for fixed cell quota of zooplankton. The model is based on the following three assumptions:

Assumption 1. The total phosphorus content $P$ in the system is fixed.

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

Assumption 3. The total phosphorus content $P$ is divided into two pools: nutrient in the phytoplankton and zooplankton.

Remark 1. The zooplankton's ingestion rate $f(x)$ satisfies the following assumptions:

In addition, under these assumptions for $f(x)$, it follows that $f(x)/x$ has the following properties:

Remark 2. Based on ^[42], Assumption 2, often referred to as "strict homeostasis", posits that heterotrophs (such as zooplankton) maintain a constant nutrient composition despite variations in their diet. This is based on the understanding that the stoichiometric variability of heterotrophs is generally much lower than that of autotrophs.

Environmental fluctuations exhibit a characteristic dependence on time, with seasonal periodicity being a prevalent temporal pattern in ecosystems, exerting substantial influence on the dynamics of organisms. Consequently, an authentic model should duly consider this seasonal effect. Building upon the Lotka-Volterra type LKE model, we incorporate $K(t)$ as the carrying capacity of phytoplankton, which is contingent upon the season-driven light intensity. Therefore, the augmented realism in our model, encompassing the influence of season-driven light intensity, is delineated by

where $K(t)$ is assumed to be a periodic function with period one year. The parameters and their biological significance are listed in Table 1.

Remark 3. ^[40] There are two minimum functions in (2.2). The first minimum function is

which expresses that the maximum biomass of phytoplankton is limited by two factors, light limitation (represented by $K(t)$, the carrying capacity influenced by season-driven light intensity) and phosphorus limitation (specifically, the availability of phosphorus in the environment). The second minimum function is

which indicates that the quality of their food influences the conversion efficiency of zooplankton.

In order to facilitate the discussion below, we introduce the following transformation

to scale (2.2) into more concise form

where

3.   Dynamics of (2.3)

In this section, the dynamics of (2.3) are analyzed, including boundary dynamics and internal dynamics. We introduce the following notations,

Theorem 1. The set

is positively invariant with respect to (2.3), where $k = \min\{K_{\max}, P/q\}$.

Proof. To prove this theorem, we need only examine that the solution $(x(t), y(t))$ of (2.3) with initial conditions $(x(0), y(0))\in \Omega$ remains in $\Omega$ for $t\geq 0$.

Assume that there exists a $t_1 > 0$ at which $(x(t_1), y(t_1))$ first crosses the boundary of $\Omega$.

Case 1. On the left boundary $\{(x, y): x = 0, \; 0 < y < P\}$, let $y_1 = \max_{t\in[0, t_1]}y(t) < P$, and assume that $x(t_1) = 0$, then we have

where $\alpha$ is a constant. Thus, $x(t)\geq x(0)\exp\{\alpha t_1\} > 0$, which contradicts $x(t_1) = 0$. Hence, no trajectory starting in $\Omega$ can across the left boundary of $\Omega$.

Case 2. On the right boundary $\{(x, y): x = k, \; 0 < qx+y < P\}$, we assume $x(t_1) = k$, then

A standard comparison argument reveals that $x(t) < k, \; t\in[0, t_1]$ and hence, no trajectory touches the right boundary of $\Omega$.

Case 3. On the bottom boundary $\{(x, y):0 < x < k, \; y = 0\}$, assume that $y(t_1) = 0$, then

thus $y(t)\geq y(0)\exp\{-dt\} > 0$, which contradicts $y(t_1) = 0$. Therefore, all orbits starting in $\Omega$ cannot escape {from} the bottom boundary of $\Omega$.

Case 4. On the top boundary $\{(x, y):0 < x < k, \; qx+y = P\}$ of $\Omega$, suppose that

Note that $t_1 > 0$ is the first time that the trajectory touched the top boundary, thus for any $t\in[0, t_1)$, $qx(t)+y(t) < P$, and then, it yields that

By (3.1), we obtain $x(t_1) = (P-y(t_1))/q$, and substitute it into (2.3), it follows that

For (3.1), it follows that $q = (P-y(t_1))/x(t_1)$. Then

From (3.3) and (3.4), it follows that

which contradicts (3.2). The proof is complete. □

Theorem 2. Let $(x(t), y(t))$ be the solution of (2.3) through $x(t_0) > 0$ and $y(t_0) > 0$, then

that is, (2.3) is dissipative.

Proof. By a standard comparison argument, it is trivial to obtain that $\liminf_{t\rightarrow \infty}x(t)\geq 0$, $\limsup_{t\rightarrow \infty}x(t)\leq k$, and $\liminf_{t\rightarrow \infty}y(t)\geq 0$.

By (2.3), note that

then we have $\limsup_{t\rightarrow \infty}y(t)\leq P$. Therefore, (2.3) is dissipative. □

Theorem 3. The extinction equilibrium $E_0 = (0, 0)$ always exists, and when $P\leq q K_{\min}$, there exists a phytoplankton-only equilibrium $E_1 = (k, 0)$. Furthermore, if $P\leq q K_{\min}$, then $E_0$ is an unstable saddle. If $P\leq q K_{\min}$ and $\min\{\hat{e}, q\}f(P/q) < d$, then $E_1$ is locally asymptotically stable. Moreover, if $\hat{e}f(P/q) < d$, then $E_1$ is globally asymptotically stable.

Proof. When $P\leq q K_{\min}$, one has $k = \min\{K_{\max}, P/q\} = P/q$, and (2.3) can be transformed into

where

In order to derive the equilibria, we need to solve

Obviously, the boundary equilibria are $E_0 = (0, 0)$ and $E_1 = (k, 0)$, see Figure 1.

Consider the Jacobian matrix of (3.5),

Then

It demonstrates that $E_0$ is an unstable saddle.

The Jacobian matrix at $E_1$ reads

Consequently, when

the eigenvalues of $J(E_1)$ are both negative, that's to say $E_1$ is locally asymptotically stable.

Based on (3.5), we deduce

When $\hat{e}f(P/q) < d$, according to comparison theorem, one can obtain that $y(t)\rightarrow 0$ as $t \rightarrow \infty$. According to the theory of asymptotically autonomous systems ^[45], (3.5) can be simplified to following limiting system

This implies that when $t\rightarrow \infty$, $x(t)\rightarrow P/q$. As a result, $E_1$ exhibits global attractivity, indicating $E_1$ is globally asymptotically stable. □

Theorem 4. Assume that

then internal equilibrium $E_2 = (x_2, y_2)$ exists. If

then $E_2$ is locally asymptotically stable.

Proof.

The internal equilibrium solves

By solving $F(x, y) = 0$, it produces that

Based on the expression of $y = g(x)$, it can be characterized as following two cases.

● $r < f'(0)P$. It is trivial that $g(0) = r/f'(0)$ and $g(k) = 0$. Then the $y = g(x)$ is a smooth curve connecting the initial point $(0, r/f'(0))$ and the ending point $(k, 0)$.

● $r\geq f'(0)P$. Similarly, one gets $g(0) = P$ and $g(k) = 0$. Thus the $y = g(x)$ is a smooth curve which connects the initial point $(0, P)$ and the ending point $(k, 0)$.

We first consider the case that $\hat{e} > q$, at which $G(x, y) = 0$ means

Based on Figure 1 and the analysis above, it is obvious and intuitional that internal equilibrium $E_2$ exists.

If $E_2$ lies below the line $\hat{e}x+y = P$, then $x_2 = f^{-1}(d/\hat{e}), \; y_2 = g(x_2)$, and the Jacobian matrix of $E_2$ is

The determinant and trace are

If (3.6) holds, then $\text{Tr} J(E_2) < 0$, therefore, both eigenvalues of $J(E_2)$ have negative real parts, and $E_2$ is locally asymptotically stable.

If $E_2$ lies above the line $\hat{e}x+y = P$, then $x_2, \; y_2$ satisfy

In this case, the Jacobian matrix at $E_2$ is

Note that (3.6), direct calculation yields that $\text{Tr} J(E_2) < 0$. Consequently, $E_2$ is locally asymptotically stable.

Next, we examine the scenario when $\hat{e}\leq q$. In this particular case, $E_2$ always lies below the line $\hat{e}x+y = P$. By employing similar arguments as preseneted above, it can be concluded that $E_2$ is locally asymptotically stable. □

In order to explore the existence and global asymptotically stability of the positive periodic solution of (2.3), we introduce some concepts and results from ^[46] that are basic for the following discussion.

Let $X, \; Z$ be normed vector spaces, $L:\text{Dom} L\subset X\rightarrow Z$ be a linear mapping, and $N:X\rightarrow Z$ be a continuous mapping. The mapping $L$ is called a Fredholm mapping of index zero if $\text{dim Ker}L = \text{codim Im} L < +\infty$ and $\text{Im} L$ is closed in $Z$. If $L$ is a Fredholm mapping of index zero, then there exist continuous projectors $R:X\rightarrow X$ and $Q:Z\rightarrow Z$ such that $\text{Ker} R = \text{Ker} L$, $\text{Im} L = \text{Ker} Q = \text{Im}(I-Q)$. It follows that $L|\text{Dom} L\cap \text{Ker} R:(I-R) X\rightarrow \text{Im}L$ is invertible. We denote the inverse of that map by $K_R$. If $\Omega$ is an open bounded subset of $X$, the mapping $N$ will be called $L$-compact on $\bar{\Omega}$ if $QN(\bar{\Omega})$ is bounded and $K_R(I-Q)N:\bar{\Omega}\rightarrow X$ is compact. Since $\text{Im} Q$ is isomorphic to $\text{Ker} L$, there exist isomorphisms $J:\text{Im} Q\rightarrow \text{Ker} L$.

Lemma 1. (^[46]) Let $L$ be a Fredholm mapping of index zero and let $N$ be $L$-compact on $\bar{\Omega}$. Suppose

(a) for each $\lambda\in(0, 1)$, every solution $x$ of $Lx = \lambda Nx$ is such that $x\notin\partial\Omega$;

(b) $QNx\neq 0$ for each $x\in \partial\Omega\cap \text{Ker} L$ and $\text{deg}\{JQN, \Omega\cap\text{Ker} L, 0\}\neq0$. Then the equation $Lx = Nx$ has at least one solution lying in $\text{Dom} L\cap \Omega$.

Theorem 5. If $\exp\{2r\omega\}P\leq \hat{e} f^{-1} (d/\hat{e})$ and $d < Pf'(0) < r+d$, then (2.3) has at least one strictly positive $\omega$-periodic solution.

Proof. Considering the biological significance of (2.3), we set $x(0) > 0$ and $y(0) > 0$. Let $x(t) = \exp\{x_1(t)\}$ and $y(t) = \exp\{x_2(t)\}$, then (2.3) turns into

Define

and $||s|| = ||(x_1, x_2)^{\text{T}}|| = \max_{t\in[0, \omega]}|x_1(t)|+\max_{t\in[0, \omega]}|x_2(t)|$ for any $s\in X$ (or $Z$). It is not difficult to verify that $X$ and $Z$ are Banach spaces with norm $||\cdot||$. Let

Then, $\text{Ker} L = \mathbb{R}^2$, $\text{Im} L = \left\{z\in Z:\int_0^\omega z(t)\text{d} t = 0\right\}$ is a closed subset of $Z$, and $\text{dim Ker} L = 2 = \text{codim Im} L$. In addition, $R$ and $Q$ are continuous mappings such that $\text{Im} R = \text{Ker} L$, $\text{Ker} Q = \text{Im} L = \text{Im} (I-Q)$. As a result, the inverse mapping (to $L$) exists, that is, $K_R:\text{Im}L\rightarrow \text{Dom} L\cap \text{Ker} R$ and

From Arzela-Ascoli theorem, it follows that $N$ is $L$-compact on $\bar{\Omega}$, where $\Omega\in X$ is an open bounded set.

In order to apply Lemma 1, we need to find an appropriate open bounded set $\Omega$. Consider the operator equation $Ls = \lambda Ns, \; \lambda\in(0, 1)$, that is,

Set $s = s(t)\in X$ be the solution of (3.8) for some $\lambda \in(0, 1)$, then integrating (3.8) from $0$ to $\omega$ yields that

and

Moreover, we derive

From (3.9)–(3.12), one can conclude

Note that $s(t) = (x_1(t), x_2(t))^{\text{T}}$, then there are $\xi_i$, $\eta_i$ such that

By (3.11), we have

which means $x_1(\xi_1)\leq\ln[K_{\max}]$, then

By utilizing (3.11) again, we derive

then $x_2(\xi_2)\leq\ln [r\exp\{K_1\}/f(\exp\{K_1\})]$, and

In addition,

which illustrates $x_1(\eta_1)\geq \ln[f^{-1}(d/\hat{e})]$, then

Therefore, $\max_{t\in[0, \omega]}|x_1(t)|\leq \max\{|K_1|, |K_3|\}: = K_5$. Furthermore, note that $\exp\{2r\omega\}P\leq \hat{e} f^{-1} (d/\hat{e})$, we have $P\leq \hat{e}\exp\{x_1(\xi_1)\}$, then

Consequently, $x_2(\eta_2)\geq\ln[P-d/f'(0)]$, and

Therefore, $\max_{t\in[0, \omega]}|x_2(t)|\leq \max\{|K_2|, |K_4|\}: = K_6$. Obviously, $K_i(i = 1, \cdots, 6)$ are independent of $\lambda$.

Next, consider the following algebraic equations

Based on the assumption $P > \hat{e}K_{\max}$, the solutions of these algebraic equations lie in following cases:

● If $\hat{e}v_1+v_2\leq P$, then there exists at least one solution. Denote by $(v_1^*, v_2^*)$ the largest solution, and it satisfies $v_1^* = f^{-1}(d/\hat{e})$, $0 < v_2^* < P$.

● If $\hat{e}v_1+v_2 > P$, $P < qK_{\min}$, then there exists at least one solution. Denote by $(v_1^{**}, v_2^{**})$ the largest solution, which satisfies $0 < v_1^{**} < P/q$, $0 < v_2^{**} < rP/(r+d)$.

● If $\hat{e}v_1+v_2 > P$, $P\geq qK_{\min}$, then there exists at least one solution $(\hat{v}_1, \hat{v}_2)$, satisfying $0 < \hat{v}_2 < P,$

Therefore, for the solution $(\tilde{v}_1, \tilde{v}_2)$ of the algebraic equations, there exists a sufficiently large number $K_7$ such that $||(\ln\tilde{v}_1, \ln\tilde{v}_2)^{\text{T}}|| = |\ln\tilde{v}_1|+|\ln\tilde{v}_2|\leq K_7$. Define $\Omega = \{s(t) = (x_1(t), x_2(t))^{\text{T}}\in X:||s|| < K\}$, where $K = K_5+K_6+K_7$, then $\Omega$ satisfies condition (a) in Lemma 1. When $s\in\partial\Omega\cap\text{Ker}L = \partial\Omega\cap \mathbb{R}^2$, $s$ is a constant vector in $\mathbb{R}^2$ with $||s|| = K$. Then $QNs\neq 0$.

For $v\in[0, 1]$, construct homotopic maps

and

Obviously, $H_v(x_1, x_2)$ satisfies that $H_0(x_1, x_2) = G(x_1, x_2)$, $H_1(x_1, x_2) = QN(x_1, x_2)$, $0\notin H_v(\partial\Omega\cap \text{Ker}L)$. In addition, under the condition in Theorem 5, we have

where $J$ is an identity map due to $\text{Im}P = \text{Ker}L$. Hence, $\Omega$ satisfies all the conditions in Lemma 1, and (3.7) has at least one solution $s^*(t) = (x_1^*(t), x_2^*(t))$ in $\text{Dom}L\cap \bar{\Omega}$. Set $x^*(t) = \exp\{x_1^*(t)\}$ and $y^*(t) = \exp\{x_2^*(t)\}$, then $(x^*(t), y^*(t))^\text{T}$ is a $\omega$-periodic positive solution of (2.3). The proof is complete. □

Lemma 2. (^[47]) Let $h$ be a real number and $f$ be a nonnegative function defined on $[h, +\infty)$ such that $f$ is integrable on $[h, +\infty)$ and is uniformly continuous on $[h, +\infty)$. Then $\lim_{t\rightarrow +\infty}f(t) = 0$.

Theorem 6. Let $(x^*(t), y^*(t))$ be the positive $\omega$-periodic solution of (2.3). If

then $(x^*(t), y^*(t))$ is globally asymptotically stable.

Proof. Let $(x(t), y(t))$ be the solution of (2.3) with positive initial values $x(t_0) > 0$ and $y(t_0) > 0$. We have proved that the set $\Gamma = \{(x, y):0\leq x\leq k, 0\leq y \leq P\}$ is an ultimately bounded region of (2.3). Then there exists a $T_1 > 0$ such that $(x(t), y(t))\in \Gamma$ and $(x^*(t), y^*(t))\in \Gamma$ for $t\geq t_0+T_1$. Consider the Lyapunov function defined by

Along the trajectory of (2.3), calculating the right derivative of $V(t)$ yields that

where $\eta$ is a positive constant. Integrating the above differential inequality from $t_0+T_1$ to $t$ produces

Then, it follows that

which means $|x(t)-x^*(t)|+|y(t)-y^*(t)|\in L^1([t_0+T_1, +\infty))$.

Note that $x^*(t)$, $y^*(t)$ are bounded, and $x(t)$, $y(t)$ are ultimate bounded in $\Gamma$, then $x(t), \; y(t), \; x^*(t),$ $y^*(t)$ all have bounded derivatives for $t\geq t_0+T_1$. Hence, $|x(t)-x^*(t)|+|y(t)-y^*(t)|$ is uniformly continuous on $[t_0+T_1, +\infty)$. From Lemma 2, it follows that

The proof is complete. □

The following theorem is established for the boundary periodic solution of Eq (2.3).

Theorem 7. If $P > qK_{\min}$, then (2.3) has a boundary periodic solution $(x_*(t), 0)$, where

Furthermore, if

then $(x_*(t), 0)$ is globally asymptotically stable.

Proof. Let $(x(t), y(t))$ be any solution of (2.3) with positive initial values $x(t_0) > 0$ and $y(t_0) > 0$. Since $\Gamma = \{(x, y):0\leq x\leq k, 0\leq y \leq P\}$ is an ultimately bounded region of (2.3), there exists a $T_1 > 0$ such that $(x(t), y(t))\in \Gamma$ and $(x_*(t), 0)\in \Gamma$ for $t\geq t_0+T_1$. From the discussion in Theorem 6, to illustrate $(x_*(t), 0)$ is globally asymptotically stable, we only need to prove that

Consider the following Lyapunov function

Calculating the right derivative of $V(t)$ along the solution of (2.3) yields that

where $\gamma$ is a positive constant. The remaining part of the proof follows the same steps as those discussed in Theorem 6. For the sake of brevity, we omit it here. □

4.   Numerical simulations

In this section, we present several numerical simulations based on the theoretical findings discussed in the preceding sections. These simulations validate the theoretical results and investigate the influence of seasonal light intensity, phosphorus availability, and zooplankton loss rate on the system dynamics. In this context, we select the consumption rate function $f(x)$ and phytoplankton carrying capacity $K(t)$ as the form

then $K(t)$ is a season-dependent function, attaining its maximum during summer and minimum during winter (Figure 2). The parameter values used in our simulations are given in Table 1, unless explicitly stated otherwise.

Initially, we conduct simulations to validate the sufficient theoretical criteria from the preceding section. As illustrated in Figure 3(a), when Theorem 3 is satisfied, the phytoplankton persists, and the zooplankton becomes extinct. In Figure 3(b), both phytoplankton and zooplankton exhibit persistence in constant eventually. Additionally, Figure 3(c) reveals the existence of a strictly positive $\omega$-periodic solution $(x^*(t), y^*(t))$ for (2.3), which is globally asymptotically stable. Furthermore, Figure 3(d) demonstrates the presence of a globally asymptotically stable boundary $\omega$-periodic solution $(x_*(t), 0)$ for (2.3), indicating periodic oscillation in phytoplankton density while zooplankton remains extinct.

The intensity of light and the availability of phosphorus ($P$) can modulate the quality of phytoplankton, as indicated in the preceding analysis. Particularly, the phosphorus availability of phytoplankton exerts a substantial influence on the predation dynamics of zooplankton and is crucial for the overall behavior of (2.3). Additionally, zooplankton impose top-down control on phytoplankton in aquatic ecosystems, thereby playing a pivotal role in regulating phytoplankton outbreaks. The loss rate of zooplankton, denoted by $d$, holds vital significance for the persistence of zooplankton population. Consequently, we select the phosphorus availability $P$ and the loss rate of zooplankton $d$ as bifurcation parameters to investigate the dynamic responses of (2.3) under continuous variations of these parameters. Figure 4 exhibits the bifurcation surfaces, illustrating the densities of all populations concerning the implicit variation of the loss rate and phosphorus availability, represented by $d$ and $P$, respectively. From Figure 4, one observes that phytoplankton grow rapidly with enhancing phosphorus and loss rate of zooplankton initially. As phosphorus availability continues to increase, the density of phytoplankton undergoes pronounced oscillations when the loss rate is low, while it displays regular periodic oscillations when the loss rate is high (see Figure 4(a)). Furthermore, zooplankton cannot sustain at significantly low phosphorus availability and high loss rate. When the loss rate is low, increasing phosphorus availability allows zooplankton to persist stably or cyclically (see Figure 4(b)).

To gain a comprehensive understanding of the bifurcation characteristics and to further investigate total phosphorus availability ($P$) and loss rate ($d$), we conduct in-depth analyses by examining specific cross-sections along the gradients of $P$ from $0$ to $0.12$ and $d$ from $0$ to $0.6$. In each bifurcation diagram, the blue curve represents the maximum biomass values, and the red curve represents the minimum biomass values.

In scenarios where the loss rate of zooplankton is low ($d = 0.05$), both phytoplankton and zooplankton coexist in a stable or cyclic state. That is, when $0 < P < 0.012$, phytoplankton growth is constrained exclusively by phosphorus availability. During this interval, phytoplankton and zooplankton can coexist within a stable steady-state. On the other hand, as the phosphorus availability spans from $0.012$ to $0.12$, the growth of phytoplankton becomes constrained by the season-driven light intensity $K(t)$. The populations can coexist via cyclic oscillation, as Figure 5 depicts. More specifically, in the $0 < P < 0.012$ range, the zooplankton biomass exhibits an increase with the rise in phosphorus availability, reflecting the positive impact of enhanced food quality. Moreover, the phytoplankton biomass experiences an initial increase followed by a decrease attributed to the predation by zooplankton. Figure 5 demonstrates that at $P = 0.012$, the internal equilibrium loses its stability and undergoes a supercritical Hopf bifurcation, leading to the coexistence of all populations through oscillations as $P$ increases from $0.012$ to $0.017$. However, as $P$ continues to expand further (e.g., $0.017 < P < 0.12$ in Figure 5), abundant phosphorus availability induces erratic oscillatory behavior in all populations.

From Figures 6 and 7, it is evident that the irregular oscillations gradually diminish as the loss rate of zooplankton increases. When $P < 0.012$ in Figure 6(b), the low food quality leads to the extinction of zooplankton. Furthermore, in Figure 6(a), the phytoplankton biomass increases with the rise in phosphorus availability. As phosphorus availability further increases, food quality improves, enabling zooplankton survive in the cyclic oscillation form. Combining the observations from Figure 6(a) and (b), the system undergoes a dynamic trajectory with the following pattern: $E_1$-existence equilibrium $(0 < P < 0.003)$ $\rightarrow$ $(x_*(t), 0)$ boundary periodic solution $(0.003 < P < 0.012)$ $\rightarrow$ $(x^*(t), y^*(t))$ positive periodic solution $(0.012 < P < 0.12)$, which is similar to the one observed when $d = 0.5$ (Figure 7).

Next, we delve into the impact of varying the loss rate $d$ on the dynamics of (2.3), focusing on high and low levels of phosphorus availability denoted by $P = 0.12$ and $P = 0.0015$, respectively. For $P = 0.12$, the phosphorus availability is at a high level, and the abundance of high-quality food benefits the survival of zooplankton. In this regime, phytoplankton and zooplankton coexist in irregular oscillation when $d < 0.51$. However, as the loss rate surpasses $0.51$, food quality becomes a significant limitation for zooplankton survival, leading to extinction. Consequently, the system transitions to a boundary periodic solution, and the biomass of phytoplankton continues to increase while maintaining regular oscillations.

From Figure 9, it is evident that when the phosphorus availability is at a low level ($P = 0.0015$), due to the low food quality, both phytoplankton and zooplankton coexist at low densities in the range of $0 < d < 0.15$. As $d$ further increases, the zooplankton biomass declines, while the phytoplankton density increases due to the absence of zooplankton predation. For $d > 0.3$, zooplankton go extinct due to the low food quality. When considering the combined information from Figures 8 and 9, it becomes evident that irregular oscillation disappears with decreasing phosphorus availability, and the threshold value of $d$ for zooplankton extinction declines. This observation is attributed to the low food quality not adequately supporting the growth and survival of zooplankton. Based on the above bifurcation diagrams concerning $P$ and $d$, Figure 10 presents time-series solutions of (2.3) under various parameter settings. These variations demonstrate that (2.3) can reflect the diverse bloom patterns seen in natural ecosystems. In LKE model, when $0.56 < K < 0.98$, the system exhibits periodic oscillatory behavior ^[40]. However, this oscillatory behavior is limited to a single type and has a very short period, making it unsuitable for describing the periodic blooms of phytoplankton, as shown in Figure 11. Therefore, it is necessary to incorporate season-driven light intensity in ecological modeling.

5.   Conclusions and discussion

Ubiquitous seasonal periodic variations in the natural environment profoundly influence numerous abiotic factors, including temperature and light intensity parameters. Light intensity significantly impacts phytoplankton, influencing their growth and dynamics. For example, light intensity affects the species composition of the phytoplankton community, the vertical distribution of phytoplankton, and their uptake of nutrients. One of the most important effects is that increased light intensity can stimulate photosynthesis, enhancing phytoplankton biomass. We comprehensively introduce a seasonal-driven light intensity factor to apprehend the ramifications of periodic fluctuations on organisms and proposes a phytoplankton-zooplankton model represented as (2.2). This model emerges as a nonautonomous food chain model, characterized by stoichiometric constraints, thus encapsulating the intricate phytoplankton heterogeneity, the connotation of food quality, and the dynamic nature of time-varying parameters. Employing qualitative analysis, the plausible attractors inherent in (2.3) encompass boundary equilibrium, internal equilibrium, boundary periodic solutions, and positive periodic solutions (as demonstrated in Figure 3). The ensuing section delves into a systematic exploration of the effects of phosphorus availability and the loss rate of zooplankton on the intricate dynamics governed by (2.3), elucidated through comprehensive numerical simulations.

The simulation outcomes have illustrated how seasonal-driven light intensity and nutrient-limited grazing combined to create blooms. In Figure 3(c), phytoplankton biomass peaks around the 200th day, while zooplankton biomass reaches its lowest point simultaneously due to low-quality food. Then, we carry out some bifurcation diagrams to demonstrate the pivotal significance of phosphorus availability and loss rate of zooplankton in shaping the growth dynamics and coexistence patterns of phytoplankton and zooplankton.} In scenarios where the loss rate $d$ remains constant (Figures 5–7), elevated phosphorus availability notably fosters the concurrent presence of phytoplankton and zooplankton. Specifically, when the zooplankton loss rate is set at a low value ($d = 0.05$), an initial stable equilibrium arises, marked by the coexistence of both phytoplankton and zooplankton. This state primarily results from phytoplankton growth being delimited by phosphorus availability. Subsequently, as phosphorus availability escalates, the development of phytoplankton becomes influenced by the seasonally fluctuating light intensity. Consequently, both phytoplankton and zooplankton manifest cyclic oscillations. Notably, when phosphorus availability surpasses a threshold, say $P > 0.017$, the proliferation of abundant phosphorus precipitates erratic oscillatory behavior across all population categories.

Insights derived from Figures 6 and 7 distinctly illustrate a noticeable trend: Irregular oscillations gradually diminish as the loss rate of zooplankton increases.} Conversely, when phosphorus availability remains consistently high, the initial coexistence of phytoplankton and zooplankton is characterized by irregular oscillations. However, as the loss rate surges, zooplankton face extinction. In the context of $P = 0.0015$, the densities of phytoplankton and zooplankton remain persistently low, owing to the restricted availability of phosphorus. Upon intensifying the loss rate, zooplankton density undergoes further decline, while phytoplankton density increases due to the absence of zooplankton predation. Despite ample food, zooplankton biomass fails to increase significantly due to poor sustenance quality. In LKE model, when $0.56 < K < 0.98$, the system exhibits periodic oscillatory behavior ^[40]. However, this oscillatory behavior is limited to a single type and has a very short period, making it unsuitable for describing the periodic blooms of phytoplankton, as shown in Figure 11. When introducing season-driven light intensity, the model can produce diverse bloom patterns, as seen in Figure 10. Hence, it is significant to incorporate $K(t)$ as the carrying capacity of phytoplankton, which is contingent upon the season-driven light intensity.

We introduce season-driven light intensity into the model, exerting a notable influence on the light-dependent carrying capacity denoted as $K(t)$. Light intensity and various other biotic and abiotic factors, including phosphorus availability, predation rates, and loss rates, are frequently subject to seasonal variations or temporal fluctuations. Consequently, the model's complexity escalates, rendering the exploration of seasonal effects more intricate. Additionally, $K(t)$ was used in a specific form in the numerical simulations. We did not explore the effects of factors such as amplitude and period on the system. However, discussing the intensity and period of season-driven is also an interesting question. What is more, the authenticity of the zooplankton's ingestion rate can be enhanced by calibrating it against empirical data obtained from natural settings. This enhancement would enable the model to address real-world complexities more effectively.

Furthermore, the linearity attributed to the loss rate of zooplankton is an assumption herein. However, it is imperative to acknowledge that this loss rate is inherently influenced by and varies in response to interactions at the subsequent trophic level. Consequently, formulating a more detailed food chain model could more explicitly elucidate the mechanisms underlying phytoplankton blooms. A comprehensive elucidation of these concerns has the potential to offer novel perspectives on the subject under investigation. However, these intricate discussions are deferred to future research endeavors.

Use of AI tools declaration

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

Acknowledgments

This research is supported by the National Natural Science Foundation of China (No. 12071068).

Conflict of interest

The authors declare there is no conflict of interest.