Modeling the role of fish population in mitigating algal bloom

1.   Introduction

Blue-green algae, or bloom-forming cyanobacteria, poses a continual ecological menace to aquatic environments, exerting substantial repercussions on the fragile balance of aquatic ecosystems. This adversely affects water quality, biodiversity, and the overall environmental health ^[1]. The proliferation of algae not only disrupts the visual clarity of the water but also gives rise to an intricate cascade of events leading to anoxic conditions in the deeper realms of water bodies ^[2].

Algae issues fuel-up, when there is an excess of nutrients, primarily nitrogen and phosphorus, entering into the water reservoir. These nutrients typically come from run-off sources, such as fertilized lawns, agricultural fields, pastures, feedlots, septic tanks, etc. ^[3]. Often nutrients accumulate in lakes and ponds and increase their susceptibility to algal blooms, which can have severe consequences for aquatic ecosystems. After the death of algal cells, their decomposition consumes large amounts of dissolved oxygen in the water, leading to oxygen depletion and the subsequent death of aerobic organisms ^[3,4]. This depletion of oxygen is particularly detrimental to fish and other aquatic life that rely on oxygen-rich environments to survive. Furthermore, algal blooms often result in producing foul odors, and impairing the sensory indicators and scenic values of water bodies ^[5]. These changes not only affect the aesthetic and recreational value of water bodies but also disrupt the natural balance of aquatic ecosystems.

Additionally, the death of certain algal cells can release substantial amounts of algal toxins into the water. These toxins pose serious threats to the quality of drinking water and can adversely affect human and animal health, ^[6]. The presence of these toxins in water sources necessitates more rigorous and costly water treatment processes to ensure safety. Moreover, algal blooms can impede water treatment processes by clogging filters and generating disinfection by-products, as described in ^[7]. These complications increase operational costs for water treatment plants, further highlighting the extensive impact of algal blooms on both the environment and human life ^[8].

To effectively address the challenges posed by algal blooms, researchers have developed various methods for algae control in aquatic environments, including physical, chemical, biological, and integrated approaches ^[9]. These methods, however, often come with their own set of limitations and drawbacks. For instance, Chen et al. ^[10] used modified clay to treat algae in Sancha Lake in Sichuan, China, achieving a 98.5% removal efficiency by the third day through simulated sedimentation. However, this method resulted in the accumulation of microcystin toxins in the lake, which negatively impacted the aquatic ecosystem. Liu et al. ^[11] demonstrates the efficacy of three coagulants (polyaluminum chloride, ferric chloride, and cationic starch) in mitigating severe cyanobacterial blooms. However, the significant accumulation of nutrients in sediments following cyanobacterial sedimentation may lead to severe internal phosphorus pollution, jeopardizing environmental restoration through submerged plant cultivation. Chemical agents, such as algaecides and herbicides exhibit some algal control efficacy but pose significant risks to other aquatic organisms, contributing to secondary pollution concerns and limiting their widespread application ^[12]. Additionally, biological methods involving microorganisms, aquatic plants, and aquatic animals for algal removal entail high environmental impacts and ecosystem disturbances ^[13,14].

Biomanipulation techniques, particularly the use of filter-feeding fish and mollusks, have gained attention as a more ecological approach to managing algal blooms ^[15,16]. The utilization of filter-feeding fish, such as silver carp (Hypophthalmichthys molitrix), and mollusks, exemplified by Dreissena rostriformis bugensis, has been extensively documented for inducing a top-down effect on phytoplankton communities ^[15,17]. As algae thrive on nutrient-rich environments, the interplay between fish and these photosynthetic organisms involves complex ecological relationships. Fish act as natural regulators, exerting control through predation on algae and disrupting the conditions conducive to their overgrowth. Furthermore, in the intricate web of aquatic ecosystems, fish play a vital role in the decomposition of detritus, which is a crucial process for sustainability of water clarity and nutrient cycling. Certain fish species, such as tilapia and grass carp are known for their voracious appetite for algae. By grazing on algal biomass, these fish help to control algal biomass and limit the development of blooms, which helps to maintain a balance between nutrient levels and reduces the likelihood of excessive algal growth.

In this study, we investigate how nutrients, algae, detritus, and fish are interconnected, using a mathematical modeling approach. We aim to understand the relationships and interactions between these variables in an aquatic ecosystem. By developing a mathematical framework, we try to gain insights into how changes in one variable affect the others, and ultimately, how their dynamics impact the overall health and sustainability of the aquatic ecosystem.

This paper is structured as follows. Section 2, presents the mathematical model formulation and provides a comprehensive description of its components, including the underlying assumptions, parameters, and variables that characterize the system under investigation. In Section 3, we conduct a thorough model analysis, focusing on the existence of equilibrium points and their local stability properties. A detailed mathematical treatment is provided to establish the conditions under which these equilibria exist and to characterize their stability behavior. Section 4, extends our analysis to explore the bifurcation dynamics of the proposed mathematical model, examining how qualitative changes in the system's behavior occur as key parameters vary. Section 5, complements our theoretical findings with extensive numerical simulations, providing visual verification and deeper insights into the analytical results obtained in the previous sections. Finally, Section 6 synthesizes the key findings of our study, discussing their implications and potential applications while suggesting directions for future research.

2.   The mathematical model

Some mathematical models have aimed to quantify and comprehend the intricate dynamics underlying the relationship between nutrient levels and the proliferation of algae, offering valuable insights into the ecological ramifications within aquatic environments^[18,19]. Also, through mathematical modeling, researchers have sought to unravel the interconnected processes that drive the dynamics of nutrient-algae relationships and shed light on the intricate web of factors that influence the aquatic ecosystem's health ^{[20,21,22,23]}.

Yan et al. ^[24] introduced a hybrid model known as rough set and multidimensional cloud model to predict the trophic and nutrient status values of the water bodies. The model was tested on 24 major lakes, and the experimental findings demonstrate that this hybrid approach yields more precise assessment results compared to other widely used models. Zhang et al. ^[25] developed a structurally dynamic model for Lake Mogan, illustrating the hysteresis response of vegetation and water quality to increasing phosphorous concentration. The model encompassed nine state variables, including phosphorous in phytoplankton, zooplankton, sediment, sediment pore water, submerged plants, epiphytes, detritus, soluble reactive phosphate $({PO_4}^{2-})$, and planktivorous fish in the lake. The findings revealed that within a phosphorous concentration range of 0.16 to 0.25 $mg$ TP $l^{-1}$ (TP stands for total phosphorous), the water state transitioned from a fresh to turbid state, inducing significant alterations in submerged plants. The model accurately predicted shifts from submerged vegetation to phytoplankton at approximately 0.25 $mg$ TP $l^{-1}$ phosphorous concentration. Their observation was in accordance with the observation of ^[19]. The results led to the conclusion that the restoration of shallow lakes occurs at a much slower pace than eutrophication, and beyond the threshold concentration of phosphorous $(0.25$ $mg$ TP $l^{-1})$, the restoration of submerged plants may not be possible.

Misra ^[21] has proposed a dynamical model for a eutrophic water body, incorporating variables, such as nutrient concentration, algal population density, zooplankton population, detritus, and dissolved oxygen concentration. Unlike the approach of Voinov and Tonkikh ^[18], Misra considered nutrient input through water run-off from agricultural fields, not limiting it to detritus as the sole external nutrient source. Their model successfully simulated the eutrophication process, establishing relationships among the variables. Obtained results indicated that an increase in nutrient supply led to higher densities of algae. Additionally, a decrease in dissolved oxygen concentration was observed with an increase in detritus density. Shukla et al. ^[23] presented a mathematical model for a eutrophied water body affected by organic pollutants. Model analysis revealed that the simultaneous impact of water pollution and eutrophication resulted in a more rapid decrease in dissolved oxygen levels compared to the presence of only one phenomenon.

In the context of model formulation, we focus on an aquatic ecosystem which experiences eutrophication due to excessive algae growth triggered by the discharge of nutrients from sources, such as domestic drainage and agricultural run-off. Our model is centered on four dynamical variables:

– $N(t)$: signifies the cumulative concentration of nutrients (phosphorous, nitrogen, etc.),

– $A(t)$: denotes the density of algae,

– $S(t)$: represents the density of detritus (formed due to the death of algae),

– $F(t)$: signifies the population of fish (both herbivorous and carnivorous), at time $t > 0$.

For biological relevance, we assume that the initial values of $N$, $A$, $S$, and $F$ are nonnegative. Therefore, the state space of our proposed model is confined in a subset of $\mathbb{R}_+^4$.

Algae found in water bodies comprises a diverse assemblage of some major taxonomic groups. Various forms within this assemblage have distinct physiological requirements and respond differently to factors, such as light, temperature, and nutrient concentration. For the model formulation, we omit the consideration of the effects of light and temperature on algal growth. Further, we assume that the algal biomass is solely influenced by the availability of nutrients in the water reservoir. Here, we assume a continuous discharge of nutrients into the water body from various sources, such as domestic drainage and agricultural run-off, occurring at a constant rate $q$, and these nutrients naturally deplete over time at a rate $\alpha_0$. We further contemplate the active absorption of these nutrients by the algae to facilitate and sustain their growth. The uptake rate of nutrients by algae follows a saturating curve, meaning that as nutrient concentration increases, the rate of uptake initially rises rapidly but then levels off as the algae's nutrient absorption capacity reaches its maximum. This behavior is described mathematically by the Michaelis-Menten type interaction $k_1\frac{NA}{k_{12}+k_{11}N}$. The parameter $k_1$ is the maximum uptake rate of nutrients by the algae when the nutrient concentration is high, while $k_{11}$ and $k_{12}$ are constants related to the affinity of the algae for the nutrients and the half-saturation constant, respectively. The half-saturation constant $k_{12}$ (when $k_{11} = 1$) is the nutrient concentration at which the uptake rate is half of the maximum rate, indicating the point at which the algae's uptake mechanism is significantly efficient but not yet saturated. This implies that at low nutrient concentrations, the uptake rate is limited by the availability of nutrients, and the algae can absorb them efficiently. As the nutrient concentration increases, the algae's absorption mechanisms become saturated, and the rate of uptake approaches its maximum capacity. This saturation effect is due to the finite number of nutrient transport sites on the algal cells, which can become fully occupied at high nutrient levels, preventing further increase in the uptake rate. As we considered, the density of algae is postulated to rely entirely on the concentration of nutrients present in the aquatic ecosystem. In this context, the growth of algae is directly linked to the same Michaelis-Menten type interaction rate at which nutrients are absorbed by the algae. This relationship is characterized by a proportionality constant $\lambda_1$. Thus, the dynamics of nutrients and algae can be governed by the following set of differential equations.

Moreover, when algae complete their life-cycle, they naturally diminish, sinking to the bottom of the water body and accumulating as detritus. This organic material serves as a substrate for bacteria residing at the lake's bottom. These bacteria are pivotal in the ecosystem, as they decompose the detritus, breaking it into simpler compounds. Through this decomposition process, nutrients are released back into the water, participating in the nutrient cycle of the water body. Ultimately, the activities of these bacteria profoundly influence the ecological equilibrium and health of the water body. Thus, the dynamics of nutrients and algae, when algae is converted into detritus and ultimately into nutrients, is governed by the following equations:

Here, the constants $\alpha_1$ and $\delta$, respectively, signify the rates at which algae and detritus naturally undergo depletion. Furthermore, the transformation of algae into detritus hinges on the natural depletion of algae, while the conversion of detritus into nutrients relies on the natural decline of detritus. Hence, $\pi_1\alpha_1A$ denotes the rate of algae conversion into detritus, and $\pi_2\delta S$ represents the rate of detritus conversion into nutrients. From the second equation of model system (2.2), it is important to note that the following condition must hold to ensure the model's feasibility:

We further consider that fish exhibit growth at a rate $r_1$ due to other sources. Fish also consume algae at a rate $r_2$ and this contributes to their growth. Moreover, the growth of fish is proportionally related to algae consumption rate, which is represented by constant $\lambda_2$. The density of detritus in an aquatic ecosystem influences oxygen levels, consequently impacting the growth rate of the fish population. As detritus accumulates, its decomposition consumes oxygen, potentially leading to decreased oxygen availability for fish. This, in turn, negatively affects the growth and overall health of the fish population within the ecosystem. This impact of detritus on growth of the fish population is considered by the term $\frac{1}{1+mS}$, where $m$ measures the adverse effect of detritus on the growth of fish population. The parameter $m$ plays a crucial role in determining the severity of this effect. A larger value of $m$ indicates a more pronounced negative impact of detritus on fish growth, while a smaller value suggests a more resilient fish population that can better withstand higher detritus levels. The constant $r_0$ represents the mortality rate of fish attributed to overcrowding within the aquatic ecosystem. Following the above considerations, the dynamics of nutrient and algae along with the fish population in the considered aquatic ecosystem can be described by the following differential equations:

The schematic diagram for model system (2.4) is depicted in Figure 1, and a detailed description of parameters and their units are provided in Table 1.

The region of attraction, which encompasses all solutions starting within the positive orthant, is contained in

where $P_m = \min\left\{\alpha_0, \alpha_1(1-\pi_1), \delta(1-\pi_2)\right\}$.

3.   Model analysis

Due to the intrinsic nonlinearity embedded in the system described by (2.4), direct analysis of the model encounters significant challenges. Consequently, in this section, we adopt the qualitative analysis of model system (2.4). Our strategy involves utilizing the stability theory of differential equations. To achieve this objective, we initially showcase the feasibility of equilibrium points and subsequently explore their stability properties.

3.1.   Equilibrium analysis

Model system (2.4) exhibits four nonnegative equilibria, which are stated as follows:

(i) Algae and fish-free equilibrium point $E_0\left(\frac{q}{\alpha_0}, 0, 0, 0\right)$. In this scenario, the aquatic environment contains nutrients, while fish, algae and detritus are conspicuously absent.

(ii) Algae-free equilibrium point $E^*\left(\frac{q}{\alpha_0}, 0, 0, \frac{r_1}{r_0}\right)$. In this state of equilibrium point, the ecosystem encompasses the simultaneous presence of nutrients and fish population with the absence of both algae and detritus.

(iii) Fish-free equilibrium point $E^{**}\left(N^{**}, A^{**}, S^{**}, 0\right)$, where $N^{**} = \frac{k_{12}\alpha_1}{\lambda_1 k_{1}-k_{11}\alpha_1}$, $A^{**} = \frac{\alpha_0\left(\frac{q}{\alpha_0}-N^{**}\right)\left(k_{12}+k_{11}N^{**}\right)}{\left(k_1 N^{**}-\pi_1\pi_2\alpha_1\right)}$, and

$S^{**} = \frac{\pi_1\alpha_1A^{**}}{\delta}$, this equilibrium point is feasible provided $k_1 N^{**} > \pi_1\pi_2\alpha_1$. This equilibrium point exemplifies a situation in the aquatic system where algae and, consequently, detritus are present, while the population of fish is notably absent.

(iv) Coexisting equilibrium point $E_i^*(N_i^*, A_i^*, S_i^*, F_i^*)$. This equilibrium point unfolds the scenario, where nutrients, algae, detritus, and fish population coexist within the considered aquatic ecosystem.

The existence of equilibria $E_0$, $E^*$, and $E^{**}$ are apparent, thus we have not delved into a discussion regarding their existence here. Moreover, we establish the viability of the coexisting equilibrium point $E_i^*$ in the following discussion. Equilibrium point $E_i^*$ can be obtained by analyzing the following equations:

Substituting $S = \frac{\pi_1\alpha_1A}{\delta}$ from Eq (3.3) in Eq (3.1), we get an equation in variable $N$ and $A$, i.e.,

From Eq (3.5), it is apparent that

(i) $N = \frac{\pi_1\pi_2k_{12}\alpha_1}{k_1-\pi_1\pi_2k_{11}\alpha_1}$ is an asymptote,

(ii) if $A = 0$, we have $N = \frac{q}{\alpha_0}$,

(iii) if $N = 0$, we have $A = -\frac{q}{\pi_1\pi_2\alpha_1}$, and $\frac{dN}{dA} = -\alpha_0\left(\frac{q}{\alpha_0}- N\right)\bigg/ \left[A\left(\alpha_0 +\frac{k_1k_{12}A}{(k_{12}+k_{11}N)^2}\right)\right] < 0$ as $N < \frac{q}{\alpha_0}$.

Moreover, substituting $S = \frac{\pi_1\alpha_1A}{\delta}$ from Eq (3.3) in Eq (3.4), we obtain

Substituting the value of $F$ from Eq (3.6) in Eq (3.2), we obtain another equation in variable $N$ and $A$, i.e.,

From above Eq (3.7), it is apparent that

(i) $N = \frac{\delta\left(k_1r_0\alpha_1+r_2^2\lambda_2\right)}{r_0 m\pi_1\alpha_1 \left(\lambda_1 k_1-k_{11}\alpha_1\right)}$ is an asymptote,

(ii) $\frac{dN}{dA} = \frac{r_2}{r_0}\cdot\frac{\left(k_{12}+k_{11}N\right)^2}{\lambda_1 k_1 k_{12}}\cdot \frac{\delta\left(\lambda_2r_2\delta-r_1m\pi_1\alpha_1\right)}{\left(\delta+m\pi_1\alpha_1A\right)^2}$,

(iii) if $A = 0$, we have $N = \frac{r_0\alpha_1k_{12}+r_1r_2}{r_0\left(\lambda_1k_1-\alpha_1k_{11}\right)} > 0$,

(iv) for $N = 0$, we have $A = -\frac{\delta\left(r_0\alpha_1+ r_1r_2\right)}{r_2^2\lambda_2\delta+r_0 m \pi_1\alpha_1^2}$.

Hence, based on the preceding analysis, it can be inferred that the isoclines (3.5) and (3.7) intersect at a single point when $\lambda_2 > \frac{r_1m\pi_1\alpha_1}{r_2\delta} = \lambda_2^*$ (Figure 2(a)). In contrast, these two isoclines may intersect at two, one, or no point(s) if $\lambda_2 < \frac{r_1m\pi_1\alpha_1}{r_2\delta} = \lambda_2^*$ within the positive quadrant (Figures 2(b)–(d)). As a result, model system (2.4) may exhibit one, two, or no coexisting equilibrium point(s) depending on the value of $\lambda_2$.

3.2.   Local stability analysis

In this section, we delve into a comprehensive analysis of the stability behavior of obtained equilibria. This exploration is crucial for providing insights into the dynamical behavior of the model system (2.4). The Jacobian matrix of the proposed model system (2.4) is represented as follows:

Stability of algae and fish-free equilibrium point $E_0\left(\frac{q}{\alpha_0}, 0, 0, 0\right)$: The Jacobian matrix for model system (2.4) at $E_0$ is

with eigenvalues $\Phi_1 = -\alpha_0$, $\Phi_2 = -\delta$, $\Phi_3 = \frac{\lambda_1 k_1 q}{\alpha_0 k_{12}+q k_{11}}-\alpha_1$ and $\Phi_4 = r_1 > 0$. Here, eigenvalues $\Phi_1$ and $\Phi_2$ consistently exhibit negativity, while $\Phi_3$ may be positive or negative depending on parameter values. In contrast, $\Phi_4$ consistently retains a positive sign. Consequently, it can be deduced that algae and fish-free equilibrium point $E_0$ is invariably unstable, whenever it exists.

Stability of algae-free equilibrium point $E^*\left(\frac{q}{\alpha_0}, 0, 0, \frac{r_1}{r_0}\right)$: The Jacobian matrix for model system (2.4) around the equilibrium point $E^*$ is

with eigenvalues $\Phi_1^* = -\alpha_0$, $\Phi_2^* = -\delta$, $\Phi_3^* = -r_1$ and $\Phi_4^* = \frac{\lambda_1k_1 q}{\alpha_0k_{12}+k_{11}q}-\alpha_1-\frac{r_1r_2}{r_0}$. Here, the first three eigenvalues consistently exhibit negativity and the fourth eigenvalue is negative when $\frac{\lambda_1k_1 q}{\alpha_0k_{12}+k_{11}q} < \alpha_1+\frac{r_1r_2}{r_0}$ and positive when $\frac{\lambda_1k_1 q}{\alpha_0k_{12}+k_{11}q} > \alpha_1+\frac{r_1r_2}{r_0}$. Therefore, it can be concluded that the algae-free equilibrium point $E^*$ is unstable when $\frac{\lambda_1k_1 q}{\alpha_0k_{12}+k_{11}q} > \alpha_1+\frac{r_1r_2}{r_0}$ and stable when $\frac{\lambda_1k_1 q}{\alpha_0k_{12}+k_{11}q} < \alpha_1+\frac{r_1r_2}{r_0}$.

Stability of fish-free equilibrium point $E^{**}\left(N^{**}, A^{**}, S^{**}, 0\right)$: The Jacobian of model system (2.4) around the equilibrium point $E^{**}$ is

The matrix $J^{**}$ has one eigenvalue $\Phi^{**}_1 = \frac{r_1+\lambda_2r_2 A^{**}}{(1+mS^{**})}$ and the rest of the three eigenvalues can be obtained by solving the following cubic equation

where

Here, it is noteworthy that eigenvalue $\Phi^{**}_1$ consistently maintains a positive sign. Consequently, it can be asserted that the fish-free equilibrium point, whenever it exists is unstable.

Stability of coexisting equilibrium point $E_i^*\left(N_i^*, A_i^*, S_i^*, F_i^*\right)$: The Jacobian matrix for model system (2.4) at equilibrium point $E_i^*$ can be written as

The characteristic equation of matrix $J^*_i$ is obtained as follows:

where

Here, $\mathcal{C}_1 > 0$ and further using the Routh Hurwitz criterion, we can say that equilibrium point $E^*_i$ is stable if $\mathcal{C}_3 > 0,$ $\mathcal{C}_4 > 0$ and $\mathcal{C}_3\left(\mathcal{C}_1\mathcal{C}_2-\mathcal{C}_3\right)-\mathcal{C}_1^2\mathcal{C}_4 > 0$.

4.   Bifurcation analysis

In this section, we derive the conditions that determine whether the model system (2.4) undergoes a saddle-node bifurcation or a transcritical bifurcation. Through this exploration, we aim to elucidate the intricate dynamics and critical thresholds that govern the system's behavior during these two distinct types of bifurcations.

4.1.   Saddle-node bifurcation

The equilibrium point analysis reveals the existence of two coexisting equilibria for model system (2.4), depending on the parameter values. This suggests the possibility of a saddle-node bifurcation occurring around equilibrium point $E^*_i$. Here, we assume that there exists a $\lambda_1 = {\lambda_1}_b$, such that $\mathcal{C}_4\big\vert_{(\lambda_1 = {\lambda_1}_b)} = 0$. Consequently, Jacobian matrix $J^*_i$ has eigenvalue 0 with algebraic multiplicity 1. Let $\mathfrak{U} = \left[u_1\; u_2\; u_3\; u_4\right]^\top$ and $\mathfrak{W} = \left[w_1\; w_2\; w_3\; w_4\right]$, sequentially represent the right and left eigenvectors of matrix $J^*_i$ corresponding to the $0$ eigenvalue, where

Suppose, $\mathfrak{Q} = \left[q_1, q_2, q_3, q_4\right]^{\top}$, where $q_1$, $q_2$, $q_3$ and $q_4$ are sequentially the righthand side of $dN/dt$, $dA/dt$, $dS/dt$, and $dF/dt$ in model system (2.4). Thus,

Thus, according to the Sotomayor's theorem, the conditions required for the existence of a saddle-node bifurcation are satisfied when $\mathfrak{B}_2\neq 0$. Consequently, we can state the following theorem concerning the existence of saddle-node bifurcation for model system (2.4).

Theorem 1. For the existence of saddle-node bifurcation around the coexisting equilibrium point $E^{*}_i$, there exists a $\lambda_1 = {\lambda_1}_b$, such that $\mathfrak{B}_2\neq 0$.

Remark 1. When growth rate of algae, driven by nutrient uptake (i.e., $\lambda_1$) surpasses the threshold ${\lambda_1}_b$, such that $\mathfrak{B}_2\neq 0$, model system (2.4) undergoes a qualitative change in its dynamical behavior. Specifically, a saddle-node bifurcation indicates the point at which two equilibrium points of the system collide and annihilate each other. Ecologically, this means that if the nutrient uptake rate by algae increases to a critical value, the algae density can rapidly shift from low to high density, potentially leading to an algal bloom. This bifurcation serves as a warning signal for environmental management, indicating that a small increase in uptake rate of nutrients by algae can lead to a dramatic change in algal growth, thereby necessitating proactive measures to prevent the onset of algal blooms.

4.2.   Transcritical bifurcation

Through the stability analysis of the obtained equilibria, it is evident that Jacobian matrix $J^*$ exhibits one eigenvalue that can be positive, negative, or zero depending on the specific values of the model parameters and the remaining three eigenvalues are negative. This suggests the possible occurrence of a transcritical bifurcation around the equilibrium point $E^*$. Thus, in this section, we discuss the existence of a transcritical bifurcation for model system (2.4) around the equilibrium point $E^*$. To investigate this, we designate $\lambda_1$ as the bifurcation parameter. Suppose there exists a $\lambda_1 = {\lambda_1}_p$, such that eigenvalue $\Phi^*_4 = 0$, which gives ${\lambda_1}_p = \left(\frac{\alpha_0k_{12}+k_{11}q}{k_1 q}\right)\left(\alpha_1+\frac{r_1r_2}{r_0}\right)$. Thus, we apply the center manifold theorem ^[26] at $E^*$. We first introduce $N = y_1$, $A = y_2$, $S = y_3$ and $F = y_4$, thus model system (2.4) can be written as

The linearized matrix for system (2.4) around the equilibrium point $E^*$ at $\lambda_1 = {\lambda_1}_p$ is

Thus, the right eigenvector $(\mathcal{U})$ and left eigenvector $(\mathcal{W})$ associated with eigenvalue 0 are represented as follows

The coefficients $a$ and $b$ of Theorem 4.1 in ^[26] for system (2.4) are written as

Here, we obtained the value of $a$ and $b$ for model system (2.4) as follows:

Here, $b > 0$ and manipulating the expression for $a$, we get a threshold quantity for $\lambda_2$ (i.e., ${\lambda_2}_c$), which determines the direction of transcritical bifurcation. Thus, we summarize the results regarding direction of transcritical bifurcation in the following theorem.

Theorem 2. For ${\lambda_1} = {\lambda_1}_p$, the direction of transcritical bifurcation of model system (2.4) is forward if $\lambda_2 > {\lambda_2}_c$ and backward if $\lambda_2 < {\lambda_2}_c$, where

with

Remark 2. If the growth rate of fish due to consumption of algae (i.e., $\lambda_2$) exceeds its threshold quantity (i.e., ${\lambda_2}_c$), the bifurcation is forward, indicating a persistent presence of algae in the water ecosystem when $\lambda_1$ beyond ${\lambda_1}_p$. In this scenario, the algal population establishes itself and maintains a continuous presence in the aquatic environment. Conversely, when $\lambda_1$ falls below ${\lambda_1}_p$, we observe a mitigation of algal density, suggesting that under these conditions, the ecosystem can naturally reduce algal concentrations. However, the system dynamics become more intricate when $\lambda_2$ is less than ${\lambda_2}_c$. In this case, the bifurcation is backward, leading to a persistent algal presence for values of $\lambda_1$ greater than ${\lambda_1}_p$. Interestingly, this scenario also presents the possibility of algal persistence even when $\lambda_1$ is below ${\lambda_1}_p$, a phenomenon that is contingent upon the initial algae density in the water body.

5.   Numerical simulations

In this section, we conduct numerical simulations for the proposed model system (2.4). Despite the crucial importance of clean water on Earth, there remains a paucity of quantitative data on water pollution. Consequently, validating the model and its outcomes with real-world field data presents challenges. Nonetheless, to visualize and analytically validate the obtained results and derive insights from them, we utilize the parameter values specified in Table 1 until otherwise mentioned. {The parameter values mentioned in Table 1 for $q$, $\alpha_0$, $k_1$, $k_{11}$, $k_{12}$, $\pi_2$, $\lambda_1$, $\pi_1$, and $\alpha_1$ are comparable to those used in previous research, particularly the study of An et al. ^[27] and Shukla et al. ^[23]. The remaining parameters are unique to our model. Notably, the simulation results we obtained in this sections align well with the findings reported in ^[23,27], lending credence to our parameter choices and model formulation.

From the model analysis, we have demonstrated that the proposed model system (2.4) has either one, two, or no coexisting equilibrium points, contingent on the value of $\lambda_2$. To visualize the specific region, where exactly one coexisting equilibrium point exists or there may be two, one, or no coexisting equilibrium point, we generate a curve, where $\lambda_2 = \lambda_2^*$ (depicted as a black curve) in the $r_1-\lambda_2$ parametric plane (Figure 3). In this figure, the red region represents the area where exactly one coexisting equilibrium point exists, while the blue region illustrates the region where model system (2.4) may have either one, two, or no coexisting equilibrium point.

Further, we choose some important parameters of model system (2.4), to see their impact on the equilibrium level of algae density and fish population in the considered aquatic ecosystem. To achieve this objective, we select $\lambda_1$, $\lambda_2$, and $k_1$ to analyze the behavior of the system's variables by varying two parameters simultaneously and integrating up to $t = 100$ days. By varying $k_1$ and $\lambda_1$, we observe that the maximum uptake rate of nutrients by algae and the growth rate of algae through nutrients uptake have a positive impact on algae density and fish population (Figure 4). This implies that for large values of $k_1$ and $\lambda_1$, the equilibrium levels of both algae density and fish population are high, while for small values of $k_1$ and $\lambda_1$, the equilibrium levels of algae density and fish population are low. Moreover, by varying $\lambda_1$ and $\lambda_2$, we generate surface plots illustrating the equilibrium levels of algae density and fish population, shown in Figure 5. From Figure 5(a), we observe that for high values of $\lambda_1$ and $\lambda_2$, the fish population is high, while the density of algae is high for a low value of $\lambda_2$ and a high value of $\lambda_1$ (Figure 5(b)).

By applying the center manifold theorem, we have established that model system (2.4) undergoes a transcritical bifurcation at $\lambda_1 = {\lambda_1}_p$. Theorem 2 further establishes a connection between the direction of the transcritical bifurcation and the parameter $\lambda_2$, which represents the growth rate of fish due to algae uptake. For the given parameter values in Table 1, the value of ${\lambda_1}_p$ is determined to be 1.8443, leading to ${\lambda_2}_c\approx 0.9998$. Consequently, for $\lambda_2 > 0.9998$, the transcritical bifurcation occurs in the forward direction, whereas if $\lambda_2 < 0.9998$, its direction is backward at ${\lambda_1}_p = 1.8443$. Additionally, we generate a plot for ${\lambda_2}_c$ by varying the value of $r_2$, as shown in Figure 6. In this plot, the black curve represents the values of $\lambda_2$ where $\lambda_2 = {\lambda_2}_c$, which divide the whole $r_2-\lambda_2$ parametric plane into two regions. The green color represents the region for forward transcritical bifurcation where ${\lambda_2} > {\lambda_2}_c$, and the maroon color represents the region for backward transcritical bifurcation, where $\lambda_2 < {\lambda_2}_c$.

Furthermore, the equilibrium curve is generated in the $\lambda_1-A$ plane for a fixed value of $\lambda_2 = 1.9 > {\lambda_2}_c = 0.9998$ (Figure 7(a)). Analysis of this figure reveals that at $\lambda_1 = 1.8443$ (denoted as bifurcation point ("BP")), model system (2.4) undergoes a transcritical bifurcation in the forward direction. This observation signifies that for $\lambda_1 > 1.8443$, the proposed system exhibits a coexisting equilibrium, whereas there is no coexisting equilibrium for $\lambda_1 < 1.8443$. Additionally, at the BP, the algae-free equilibrium point transfers its stability to the coexisting equilibrium point as the bifurcation parameter $\lambda_1$ crosses its critical threshold ${\lambda_1}_p = 1.8443$. This phenomenon of forward transcritical bifurcation indicates that the density of algae in the considered aquatic ecosystem vanishes if $\lambda_1 < 1.8443$, while it consistently exists for $\lambda_1 > 1.8443$. To show this stability transfer between the algae-free and coexisting equilibrium point, we generate the variation plots in Figures 7(b), (c) for algae density in the aquatic ecosystem with respect to time "$t$". Figure 7(b) is generated for $\lambda_1 = 1 < 1.8443$, which depicts that all the solution trajectories approach the algae-free equilibrium point. Therefore, whatever be the density of algae in the considered aquatic ecosystem initially, as time $t\to \infty$, it will approach to $0$. Biologically, this implies that if $\lambda_1 < {\lambda}_p = 1.8443$, the algae completely mitigates from the aquatic ecosystem. Moreover, for ${\lambda_1} = 2.8 > 1.8443$, the variation plot for algae density with respect to time "$t$" is shown in Figure 7(c). This figure demonstrates that as time $t\to\infty$, all solution trajectories converge toward the coexisting equilibrium point. From a biological standpoint, this phenomenon illustrates that algae persists consistently in the aquatic ecosystem and reaches to its coexisting equilibrium point level for $\lambda_1 > {\lambda_1}_p$.

For $\lambda_2 = 0.3 < \lambda_{2c} = 0.9998$, the proposed system exhibits the phenomenon of transcritical bifurcation in the backward direction as shown in Figure 8(a). This figure illustrates that the equilibrium curve bends at the "SN" point (where $\lambda_1 = {\lambda_1}_b\approx 1.482$), indicating the occurrence of a saddle-node bifurcation for model system (2.4). At this "SN" point, two coexisting equilibria collide and annihilate each other with one being stable and the other unstable in nature. This scenario delineates that model system (2.4) has one coexisting equilibrium point for ${\lambda_1} > {\lambda_1}_p = 1.8443$, no coexisting equilibrium point for ${\lambda_1} < {\lambda_1}_b$ and two coexisting equilibria when $\lambda_1\in({\lambda_1}_b = 1.482, {\lambda_1}_p = 1.8443)$. Further, we generate the variation plots for algae density for $\lambda_1 = 1.519\in({\lambda_1}_b, {\lambda_1}_p)$ (blue curves), $\lambda_1 = 1.85 > {\lambda_1}_p$ (green curves) and $\lambda_1 = 1.4 < {\lambda_1}_b$ (red curves) as shown in Figures 8(b)–(d). From Figure 8(b), we can see that all the solution trajectories either gravitate toward the stable coexisting equilibrium point or stable algae-free equilibrium point. This implies that algae in the considered aquatic ecosystem mitigates or persists wholly depending on the initial density of the algae. Figure 8(c) illustrates that the algae density in the water reservoir always persists and attains its equilibrium level as $t\to \infty$. Moreover, Figure 8(d) demonstrates that the algae is completely mitigated from the considered water reservoir as $t\to \infty$.

5.1.   Sensitivity analysis

Sensitivity analysis is a valuable technique for examining the relationship between parameters and the solution of a system. The semi-relative sensitivity of a variable $V$ with respect to a parameter $p$ is given by the expression $p\frac{\partial{V}}{\partial{p}}$. This analysis provides insights into how the solution changes when the parameter is doubled. The choice to double parameters in sensitivity analysis offers a balance between simplicity and insight. This approach provides clear, comparable results across parameters and serves as an efficient initial screening method. It's computationally manageable and easily interpretable, revealing the direction and magnitude of the system's response to parameter changes. We perform the basic sensitivity analysis for the parameters $q$, $k_1$ and $\lambda_2$. These three parameters are important; therefore, we performed both bifurcation and sensitivity analysis of these parameters. The parameters $q$, $k_1$, and $\lambda_2$ are critical for sensitivity analysis in this aquatic ecosystem model due to their significant roles in system dynamics. The parameter $q$ represents the rate of nutrient input, directly affecting the nutrient concentration that drives the entire system. Parameter $k_1$, the maximum nutrient uptake rate by algae is a key parameter, which determines both nutrient absorption and algal growth rates. The parameter $\lambda_2$ is a the proportionality constant for fish growth due to algae consumption and directly influences the fish population dynamics. Changes in these parameters can have far-reaching effects on the ecosystem's balance, from nutrient levels to algal growth and fish populations. Analyzing the model's sensitivity to these parameters provides crucial insights into how the ecosystem responds to variations in nutrient loading, algal uptake efficiency, and energy transfer in the food chain, which is valuable for understanding and managing aquatic ecosystems under different conditions. The differential systems of sensitivity corresponding to system (2.4) with respect to $q$, $k_1$ and $\lambda_2$ are written as follows:

To perform a semi-relative sensitivity analysis, we have chosen the parameters $q = 4$, $k_1 = 0.4$, and $\lambda_2 = 1.8$. We have calculated the semi-relative sensitivity solutions to determine the impact of these parameters on the dynamic variables. As shown in Figure 9, doubling the parameter $q$ results in increases in the values of the dynamic variables $N$, $A$, $S$, and $F$ by 22.43, 3.86, 0.3784, and 0.4474, respectively. Similarly, doubling the parameter $k_1$ leads to a decrease in the value of the dynamic variable $N$ by 63.31, while the values of $A$, $S$ and $F$ increase by 24.39, 2.386, and 2.785, respectively. Moreover, doubling the parameter $\lambda_2$ increases the values of $N$ and $F$ by 9.391 and 0.6307, respectively, whereas the values of $A$ and $S$ decrease by 5.687 and 0.534, respectively.

Moreover, other parameters such as $\delta_1$, $r_1$, and $r_0$ have shown considerable influence on the model's behavior, with notable changes in solutions observed when their values were doubled. Thus, we have also conducted sensitivity analysis for proposed model (2.4) with respect to $\delta_1$, $r_1$, and $r_0$. The remaining parameters in the proposed model are either proportionality constants or demonstrated a less substantial impact when their values were doubled. The differential systems of sensitivity corresponding to system (2.4) with respect to $\delta_1$, $r_1$ and $r_0$ are written as follows:

To perform a semi-relative sensitivity analysis, we have chosen the parameters $\delta = 0.5$, $r_1 = 0.8$, and $r_0 = 0.4$. We have calculated the semi-relative sensitivity solutions to determine the impact of these parameters on the dynamic variables. As shown in Figure 10, doubling the parameter $\delta$ results in increase in the values of the dynamic variables $N$ and $F$ by 2.216 and 0.01154, respectively, whereas a decrease in the dynamic variables $A$ and $S$ by 0.8 and 0.2, respectively. Similarly, doubling the parameter $r_1$ leads to an increase in the value of the dynamic variable $N$ and $F$ by 0.27 and 1.02, respectively while the values of $A$ and $S$ decrease by 2.316 and 0.085, respectively. Moreover, doubling the parameter $r_0$ decreases the values of $N$ and $F$ by 8.279 and 1.5, respectively, whereas the values of $A$ and $S$ increase by 3.956 and 0.15, respectively.

6.   Discussion

In this study, the analysis using the center manifold theorem indicates that when the growth rate of the fish population attributed to the uptake of algae, exceeding a threshold value ${\lambda_2}_c$, the bifurcation proceeds in a forward direction at $\lambda_1 = {\lambda_1}_p$ (threshold value of growth rate of algae due to uptake of nutrients). Conversely, if $\lambda_2 < {\lambda_2}_c$, the bifurcation takes a backward direction at $\lambda_1 = {\lambda_1}_p$. The transcritical bifurcation in forward direction elucidates that as the growth rate of algae, propelled by nutrients uptake, surpasses a critical threshold ${\lambda_1}_p$, the equilibrium point undergoes a significant shift. At low growth rate of algae due to nutrients uptake, the algae biomass encounters challenges in flourishing within the aquatic ecosystem. However, with an increase in the growth rate of algae due to nutrients uptake, the algae biomass experiences heightened growth, rendering the coexisting equilibrium point feasible, while destabilizing the previously stable algae-free equilibrium point. Biologically, this phenomenon describes that if the growth rate of fish population propelled by the algae uptake is grater than its threshold value ${\lambda_2}_c$, the algae biomass can be totally removed from the water reservoir if $\lambda_1 < {\lambda_1}_p$ and always persist in the water reservoir when $\lambda_1 > {\lambda_1}_p$. Further, the noteworthy phenomenon of backward transcritical bifurcation manifests when the growth rate of fish due to algae uptake is lower than the threshold quantity ${\lambda_2}_c$. In this scenario, two coexisting equilibria emerge even though the growth rate of the algae biomass due to nutrients uptake is below the threshold quantity ${\lambda_1}_p$. From the emerged two equilibria, the equilibrium point with high algae density is stable, while equilibrium point with low algae density is unstable. These two equilibria collide and disappear by experiencing a saddle-node bifurcation at $\lambda_1 = {\lambda_1}_b$. From a biological perspective, this phenomenon elucidates that despite the low growth rate of algae fueled by nutrients uptake, the presence of algae biomass in the water reservoir can persist for $\lambda_1\in({\lambda_1}_b, {\lambda_1}_p)$. This persistence occurs because of low growth of the fish population driven by algae uptake and initial density of algae in the water body. The occurrence of backward bifurcation highlights that if the growth rate of the fish population, attributed to algae consumption, falls below a critical threshold value ${\lambda_2}_c$, then complete removal of algae biomass from the water reservoir is only achievable if the growth rate of algae, driven by nutrients uptake, is below its threshold quantity ${\lambda_1}_b$. This ecological insight underscores the complex interplay between nutrient-driven algae dynamics and fish population growth in shaping the equilibrium points of the aquatic ecosystem.

This study distinguishes itself in the field of algal bloom modeling by emphasizing the crucial role of fish populations in aquatic ecosystem dynamics. Unlike earlier research that focused primarily on nutrient-phytoplankton interactions, such as O'Brien's 1974 study ^[28] or Huppert et al.'s 2005 work ^[29], this study explicitly incorporates fish population growth rates as a key factor influencing algal bloom dynamics.

This research identifies critical thresholds for both algae and fish growth rates, demonstrating how fish consumption of algae can potentially control or exacerbate algal blooms under various conditions. This approach provides a more nuanced understanding compared to studies, like Chen et al.'s model ^[30] or Shukla et al.'s ^[23] work, which focused more on environmental factors and nutrient inflow respectively.

Moreover, a notable finding is the phenomenon of backward bifurcation when fish population growth falls below a critical threshold. This insight is not explored in studies, like Zhao et al.'s ^[31] on stochastic factors but reveals how low fish growth rates can contribute to algae persistence even under low nutrient conditions. While sharing some conceptual ground with recent works, like An et al.'s resource-based models ^[27], this study's focus on immediate fish-algae interactions offers a unique perspective. By explicitly modeling the fish population's role, it provides valuable insights for ecosystem management, particularly in using fish populations as a potential tool for mitigating algal bloom risks.

Although, our model provides valuable insights into algal bloom dynamics, several limitations warrant consideration. The use of hypothetical parameters, while allowing theoretical exploration, limits direct applicability to specific ecosystems. More critically, the omission of light and temperature effects, which are significant factors in algal growth impacts the model's ecological realism. Incorporating these factors could introduce seasonal variations in bloom patterns, potentially altering equilibrium stability throughout the year and leading to additional bifurcations. Fish-algae interactions would likely vary with temperature, affecting grazing rates and bloom control. Furthermore, light-dependent growth could introduce depth-based variations, necessitating a spatially explicit approach. Including seasonal effects could inform optimal timing for management interventions, while temperature considerations could elucidate system resilience to climate change. Despite these limitations, our model provides a foundational framework for understanding nutrient-algae-fish dynamics, offering valuable theoretical insights and a basis for more comprehensive future studies.

Moreover, future iterations of the model incorporating factors, like temperature and light dependent fish-algae growth would enhance its predictive power for real-world scenarios. Exploring temperature-dependent fish-algae interactions may provide a more nuanced understanding of grazing rates and bloom control mechanisms. Also, developing a spatially explicit model to account for depth-based variations in light-dependent algal growth would offer a more comprehensive view of bloom dynamics in stratified water bodies. These enhancements could inform optimal timing for management interventions and elucidate the system's resilience to climate change. While these additions would increase model complexity, they would also greatly enhance its relevance for real-world scenarios and ecosystem management strategies.

7.   Conclusions

The complex interplay among nutrients, algae, and fish within the aquatic ecosystem is crucial for shaping the dynamics of the aquatic ecosystem. Nutrients, specifically nitrogen and phosphorus, act as vital components for the sustainable growth of algae. When there is an excess of nutrient input from agricultural run-off or household discharges, it can result in an overgrowth of algae, leading to algal blooms. This relationship between nutrients and algae has direct implications on fish populations, as their sustenance relies on the availability of dissolved oxygen in the aquatic ecosystem. Fish play a significant role in the ecosystem by regulating algae levels through grazing, thereby contributing to maintaining a balance. Conversely, the abundance of algae directly influences the growth of fish populations as it impacts their primary food source. This intricate web of interactions emphasizes the fragile equilibrium necessary for a healthy aquatic ecosystem. Disruptions, whether caused by nutrient imbalances or other environmental factors, can have cascading effects, affecting the overall stability and biodiversity of the aquatic ecosystem. Therefore, it is imperative to comprehend and manage the nutrient-algae-fish relationship for sustainable aquatic ecosystem management and conservation efforts. In this article, we have formulated a four-dimensional mathematical model designed to govern the dispersion of algae in water bodies with dynamical variables including nutrients, algae, detritus, and fish. First, we have identified the feasible equilibria of the proposed mathematical model. Subsequently, utilizing the stability theory of differential equations, we have assessed the stability of obtained equilibria. Furthermore, our investigation has revealed that the proposed mathematical model undergoes transcritical bifurcation and saddle-node bifurcation. To illuminate the direction of the transcritical bifurcation, we have employed the center manifold theorem. Additionally, to show the existence of saddle-node bifurcation, Sotomayor's theorem has been applied.

Use of AI tools declaration

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

Acknowledgement

The authors gratefully acknowledge Qassim University, represented by the Deanship of Scientific Research, on the financial support for this research under the number (2023-SDG-1-BSRC36971) during the academic year 1445 AH/2023 AD.

Conflict of interest

The authors declare there is no conflicts of interest.