The dynamics of an aquatic ecological model with aggregation, Fear and Harvesting Effects

1.   Introduction

Ecology is the scientific study of how living organisms interact with one another and with their surroundings. Interactions between species and between organisms and their physical and chemical environments involve all individuals of the same type. Aquatic ecology is the study of these interactions in all aquatic environments, such as rivers, lakes, streams, ponds, wetlands, and oceans. The physical properties of aquatic organisms have an impact on all types of species found there. The organisms in a specific ecosystem are directly impacted by environmental elements like nutrient content, water flow, temperature, and shelter. Only those species that can adapt to and thrive in a specific habitat will flourish because of the resources at hand. Predation and competition for resources (such as food and habitat) have an impact on species abundance and variety in aquatic ecosystems, as allow interactions between living organisms. In turn, living creatures in a given environment can have an impact on certain characteristics of that environment ^[1,2,3,4].

We are all aware of this, algal blooms are triggered by a combination of external and internal forces, algal population aggregation and migration are examples of the latter ^[5]. In the last 10 years, numerous scholars have examined how the growth of algae affects the reproduction of aquatic organisms in environmental studies and aquatic biology ^{[6,7,8,9,10,11,12]}. The tight gel layer of microcystis aggregation in particular, filter feeding fish have a hard time digesting it. Even if the filter-feeding fish eat them, a substantial percentage of undamaged active algal cells will remain in the fish feces. Xie and Liu ^[13] have demonstrated that aggregations of microcystis can generate more algal toxins than solitary cells, reducing the danger of predation. On other hand, Sommer et.al. ^[14] have shown the food web in the aquatic environment that there are fish that feed on algae and predatory fish that feed on small fish (prey).

However, on the mentioned area, there has been little scientific investigation into the population dynamics and mathematical modeling of the dynamic impacts of algal aggregation. Furthermore, population dynamics and the ecological mathematical model a complete research technique combines coordination, organization, and purpose to create three foundations for studying biological systems. Analysis of dietary intake, the view from the model, and dynamic presentation are the three kinds (see ^[15,16]). The basic goal of population dynamics modeling is to gain a better knowledge of how populations interact and how they are affected by internal and external factors ^[17,18]. However, a thorough investigation of ecological mathematical models and population dynamics has been made in recent decade. In addition, a slew of additional articles has yielded key research findings in population dynamics and ecological mathematical models (we refer ^{[19,20,21,22,23]}).

Here we remark that authors ^[24] represented the following system

which has created a fresh aquatic ecological model to explain the microcytic aeruginosa aggregation impact and filter-feeding fish. Further, the model describes the connection between microcytic aeruginosa and filter feeding fish.

By memorizing the foregoing, and depending on the investigation carried out by Sommer et.al. ^[14], we divided the aquatic environment into the following:

According to recent studies, predators have an impact on prey populations not only by killing them directly, but also through instilling fear in prey animals. This significantly affects their pace of reproduction (see ^[25,26]). Predator fear may affect the behavioral patterns and physiology of some prey species. This is more terrible than killing someone directly ^[27,28]. Wang et al. ^[26], which has suggested that raising the cost of fear can help to stabilize the system by incorporating the fear effect into a predator-prey model. Pandey et al. ^[29] investigated a tri-topic food chain model with Holing type-Ⅱ functional response by incorporating fear impact. Mentioned authors have discovered that a low amount of fear causes chaotic behavior, whereas a high level of fear allows the chaotic oscillations to be controlled. Predator fear may affect the behavioral patterns and physiology of some prey species; this is more terrible than killing someone directly.

In this study, we have expanded the first system from the second dimension to the third dimension, where we added a very important component in the aquatic ecological model. The said target is the large predatory fish that has an effect on each of the first two components (the microcytic aeruginosa aggregation and filter feeding fish). In addition, we studied the effect of fear of predatory fish (big fish) on the filter feeding fish because the main source of growth for filters feeding from microcytic aeruginosa fish, when they were fed on microcytic aeruginosa. On other hand, the harvesting effect on the predatory fish is considered.

Keeping the above discussion in mind, we add the function of fear for food filter feeding fish from predatory fish. Therefore, the mathematical formulation of aquatic environment models takes the form

where

● $r$ is the rate of intrinsic growth of microcystis aeruginosa.

● k is the maximum environmental capacity.

● ${\alpha }_{1}$ is the grazing coefficient.

● ${m}_{1}$ is the microcystis aeruginosa aggregation parameter.

● ${c}_{1}$ is the half-saturation constant.

● ${\beta }_{1}$ is the absorption coefficient.

● $c$ is the harvesting's capacity coefficient.

● $E$ is the subject of big fish population to a combined harvesting effort.

● $\delta$ is the level of fear which divides anti- predatory fish behaviour of the filter feeding fish.

● ${\beta }_{2}$ is efficiency of conversion of consumed filter-feeding fish into big fish.

In the process of dynamic systems analysis, we try as much as possible to facilitate the calculations. Therefore, since model (2) has 11 parameters, hence, we are going to do the Non-dimensionalization for all variables and parameters as:

Hence, the dimensional system (2) corresponding the following non-dimensional system:

All parameters in model (3) are expected to be positive. The following are the initial conditions for system (3)

2.   Some fundamental outcomes

In this part, we first demonstrate the positivity and boundedness of solutions for the system (3). These are quite significant in terms of the model's validity.

Lemma 1. For every $T > 0$, all solutions $\left(x\left(T\right), y\left(T\right), z\left(T\right)\right)$ of the system (3) with initial values $\left(x\left(0\right), y\left(0\right), z\left(0\right)\right)\in {R}_{+}^{3}$, are positive.

Proof. It is simple to show that the solutions of the system (3) are nonnegative under the initial conditions: $x\left(0\right)\ge 0, y\left(0\right)\ge 0,$ and $z\left(0\right)\ge 0$.

In the first, second and third instances, we have gotten respectively as

So

Therefore, all solutions of the model (3) are positive.

The system (3) exhibits continuous interaction functions and continuous partial derivatives, indicating that the solution exists and is unique. Furthermore, to guarantee convergence of the solution to an attractor, the solution of the system (3) is proven uniformly bounded, as shown in the following theorem.

Lemma 2. The feasible region of system (3), where all solutions can be found is given by

Proof. We have the following result from the first equation: $\frac{dx}{dT}\le x\left(1-\frac{x}{{k}_{1}}\right)$, therefore, when $T\to \infty$, then $x\le {k}_{1}$. Now let us consider ${L}_{1} = x+\frac{1}{{k}_{3}}y+\frac{{k}_{4}}{{{k}_{3}k}_{5}}z$, the derivative of ${L}_{1}$ with regard to time can thus be expressed as

where W $= \mathop {{\rm{min}}}\limits_{T \to \infty } \left\{ {1, {k_6}} \right\}, {\xi _1} = su{p_{T \to \infty }}y\left( T \right)$. Then, using direct computation, one can see that for $T\to \infty$, we obtain

where ${\xi }_{2} = 2{k}_{1}+\frac{\xi }{{k}_{3}}$.

Before starting with the existence of equilibrium points of (3) and its discussion, we present the Jacobian matrix of (3) as:

3.   Existence and stability of equilibriums

The equilibrium of the system (3) is as follows:

● The trivial equilibrium point ${u}_{0}\left(\mathrm{0, 0}, 0\right)$, which always exists. To study its local stability, the characteristic equation is computed from the given Jacobian matrix

as

Clearly, (8) has three distinct roots ${\lambda }_{1} = 1, {\lambda }_{2} = -\frac{{k}_{2}{k}_{3}}{1-{k}_{2}}$ and ${\lambda }_{3} = -{k}_{6}$. The Routh–Hurwitz criterion states that the trivial equilibrium point is locally asymptotically stable, when its all eigenvalue has negative real parts. Notice that only ${\lambda }_{3} < 0$, and ${\lambda }_{1}$ is purely positive, so the trivial equilibrium point ${u}_{0}\left(\mathrm{0, 0}, 0\right)$ is saddle node point.

● The x-axis equilibrium point ${u}_{1}\left({k}_{1}, \mathrm{0, 0}\right)$, which is always exists. The characteristic equation of the following Jacobin matric around ${u}_{1}$

is given by

Clearly (10) has three distinct roots ${\lambda }_{1} = -1, {\lambda }_{2} = \frac{{k}_{3}\left({k}_{1}-{k}_{2}\right)}{{k}_{1}-{k}_{2}+1}$ and ${\lambda }_{3} = -{k}_{6}$. The Routh–Hurwitz criterion establishes that the x-axis equilibrium point is locally asymptotically stable when its all eigenvalue has negative real parts. Notice that ${\lambda }_{1}, {\lambda }_{3} < 0$, and ${\lambda }_{2}$ is negative if ${k}_{1} < {k}_{2} < {k}_{1}+1$. Hence, under this condition, the x-axis equilibrium point ${u}_{1}\left({k}_{1}, \mathrm{0, 0}\right)$ is locally asymptotically stable.

Theorem 1. The trivial equilibrium point ${u}_{0}\left(\mathrm{0, 0}, 0\right)$ and the x-axis equilibrium point ${u}_{1}\left({k}_{1}, \mathrm{0, 0}\right)$ exist always and exhibit the following local behaviors:

1. the trivial equilibrium point ${u}_{0}\left(\mathrm{0, 0}, 0\right)$ is saddle point.

2. the x-axis equilibrium point ${u}_{1}\left({k}_{1}, \mathrm{0, 0}\right)$ is locally asymptotically stable provided that:

In the next part, we have to discussed the existence of internal equilibrium point. If the internal equilibrium does not exist, the model (3) will not be durable, implying that the filter-feeding fish and predatory fish would become extinct. But this is not good for ecological balance. As a result, it is critical that we investigate the model's complicated dynamic having properties of the Model (3). We know the existence condition of equilibrium from the foregoing analysis. Therefore, the next step is to investigate stability. (Note that we only discuss equilibrium stability when it exists.)

● The interior equilibrium point is denoted by ${u}^{*}\left({x}^{*}, {y}^{*}, {z}^{*}\right),$ where

Hence, ${z}^{*}$ is positive provided that:

and ${x}^{*}$ is a root of the following third degree polynomial

where,

Therefore, according the Descartes role, the following Lemma summarizes the existence of the interior equilibrium point ${u}^{*}\left({x}^{*}, {y}^{*}, {z}^{*}\right)$.

Lemma 3: Under the condition (11), the model (2) has:

1) a unique interior equilibrium point ${u}^{*}\left({x}^{*}, {y}^{*}, {z}^{*}\right)$ where one of the following is satisfied

ⅰ- $1 > {k}_{1}+{k}_{2}, {k}_{6} > {k}_{5}\left(1-{k}_{2}\right);$

ⅱ- $1 < {k}_{1}+{k}_{2}, {k}_{6} < {k}_{5}\left(1-{k}_{2}\right);$

ⅲ- $1 > {k}_{1}+{k}_{2}, {k}_{6} < {k}_{5}\left(1-{k}_{2}\right)$.

2) A triple interior equilibrium points ${u}_{1}^{*}\left({x}_{1}^{*}, {y}^{*}, {z}^{*}\right), {u}_{2}^{*}\left({x}_{2}^{*}, {y}^{*}, {z}^{*}\right)$ and ${u}_{3}^{*}\left({x}_{3}^{*}, {y}^{*}, {z}^{*}\right)$ if the following is satisfied

ⅰ- $1 < {k}_{1}+{k}_{2}, {k}_{6} > {k}_{5}\left(1-{k}_{2}\right)$.

To simplify the calculations, let's defined the following two quantities

The Jacobin matrix (6) around ${u}^{*}\left({x}^{*}, {y}^{*}, {z}^{*}\right)$ is given by

The characteristic equations of (13) is,

where,

Given the coefficient ${A}^{*}$, which according to Routh-Hwartze Criteria must be both it and coefficient ${C}^{*}$ with a positive sign, then ${A}^{*}$ is positive using the condition

It is quite easy to verify that the same condition (15), leads to ${C}^{*}$ being positive as well. According to Ruth-Hwartz, the characteristic equation (14) represents eigenvalues with negative real parts in order for us to achieve the stability property of the positive point, if ${A}^{*}, {C}^{*}$ and ${A}^{*}{B}^{*}-{C}^{*} > 0$. Therefore, by forward calculations its easy to verify that ${A}^{*}{B}^{*}-{C}^{*} > 0$ provided the conditions (11) and (15). Hence, the positive equilibrium point ${u}^{*}\left({x}^{*}, {y}^{*}, {z}^{*}\right)$ is locally asymptotically stable.

Theorem 3. Suppose that the conditions of existence for the positive equilibrium point ${u}^{*}\left({x}^{*}, {y}^{*}, {z}^{*}\right)$, the point ${u}^{*}$ is locally asymptotically stable if the condition (15) is met.

Remark 1. The situation of extinction of all species is continually in flux. Other equilibrium states, on the other hand, are conditionally stable. Even if there is no microcystis aeruginosa in the environment, this might be taken as implying that both filter-feeding and predatory fish cannot go extinct at the same time. Because of maximum environmental capacity and microcystis aeruginosa aggregation, this occurrence may be explained by the model (3).

4.   Local bifurcation analysis

In this part, the system (3) is rebuilt in a vector form to explore the local bifurcation that may occur around each equilibrium point.

where ${g}_{i}, i = \mathrm{1, 2}, 3$ are functions which given in the right hand side of dimensionless the system (3). The general second derivative of the vector $g$ can now be represented as follows using simple computations:

here $U = {({h}_{1}, {h}_{2}, {h}_{3})}^{T}$ be any eigenvector and ${k}_{1}$ be a bifurcation parameter.

Theorem 4: Assume that the condition (15) is not satisfied. Then, system (3) near the positive equilibrium point ${u}^{*}\left({x}^{*}, {y}^{*}, {z}^{*}\right)$ has a saddle-node bifurcation.

Proof. According to (3) and condition (15), a simple calculations reveals that the Jacobian matrix of system (1), stated by (7), at ${E}_{0}$ with ${{k}_{1}\equiv {k}_{1}}^{*} = \frac{2{x}^{*}}{1-{k}_{13}^{*2}{y}^{*}}$ has zero eigenvalue ${\lambda }_{02}\left({{k}_{1}}^{*}\right) = 0$ and can be written as follows:

Assume that $V = {({v}_{1}, {v}_{2}, {v}_{3})}^{T}$ is the eigenvector to ${\lambda }_{*} = 0$, then ${{J}_{3}}^{*}V = {\bf{0}}$ yields:

here ${v}_{3}$ be any real, non-zero number.

Similarly, if $W = {({w}_{1}, {w}_{2}, {w}_{3})}^{T}$ is the eigenvector to, ${\lambda }_{*} = 0$ of the matrix ${\left({{J}_{3}}^{*}\right)}^{T}$, then ${\left({{J}_{3}}^{*}\right)}^{T}W = {\bf{0}}$ results in

where ${w}_{3}$ be any real number, not zero. Moreover, it is easy to verify that

Then

System (3) near u* features saddle – node bifurcation, according to Sotomayer's theorem.

Now by substitute ${u}^{*}$ and ${{k}_{1}}^{*}$ in (16) with $V$ instead of $H$ it's observed that

Accordingly, we obtain that:

Hence due to Sotomayer's theorem a saddle–node bifurcation is occur.

5.   Numerical simulation

Some numerical simulations are offered in this part to back up our analytical and mathematical findings. These simulations also reveal the fascinating complicated of the system's behavior. In the following section, we will show several examples to back up our earlier findings. For a hypothetical collection of data, numerical simulations are run using the matlab 7.0.1 software suite.

Suppose that the following hypothetical collection of data;

Looking at the scalar hypothetical set (22), we notice that, it's given the coefficients of the polynomial ${R}_{1} = 0.15 > 0$ and ${R}_{2} = -0.0068 < 0$. By solving model (3) we have, ${u}^{*} = \left(0.4402, 0.6364, 0.0833\right)$. Since, ${k}_{2} < \mathrm{min}\left\{0.4402, 1\right\}$, so the positive condition ${k}_{2} < \mathrm{min}\left\{{x}^{*}, 1\right\}$ for ${z}^{*}$ is satisfied. Therefore, according the condition (11) and point (1-ⅰ) of lemma 3 all the conditions of existence for the positive point are fulfilled. On other hand, the local stability condition (15) is satisfied. Hence, Figure 1(a) and Figure 1(b) show us the positive equilibrium point ${u}^{*}\left({x}^{*}, {y}^{*}, {z}^{*}\right)$, which appears to be locally asymptotically stable.

Suppose that the following hypothetical collection of data;

Looking at the scalar hypothetical values given in (23), we notice that at the coefficients of the polynomial ${R}_{1} = -7 < 0$ and ${R}_{2} = 10.4 > 0$. By solving model (3), we have, ${u}^{*} = \left(5.4396, 1.4857, \;0.1654\right)$. Since, ${k}_{2}\le \mathrm{min}\left\{5.4396, 1\right\}$, so the positive condition ${k}_{2} < \mathrm{min}\left\{{x}^{*}, 1\right\},$ for ${z}^{*}$ is satisfied. Therefore, according to the condition (11) and point (1-ⅱ) of Lemma 3, all the conditions of existence for the positive point are fulfilled. On other hand, the local stability condition (15), at $\frac{{k}_{2}}{2} = 3.5 < 5.4396$ is satisfied. Hence, Figure 2 (a) and (b) shows us the positive equilibrium point ${u}^{*}\left({x}^{*}, {y}^{*}, {z}^{*}\right)$, which appears to be locally asymptotically stable.

Suppose that the following hypothetical collection of data

The reader can notice that the imposed set (24) that contains ${k}_{1} < {k}_{2} < {k}_{1}+1$ makes the second root of the polynomial (10) with a negative sign. Therefore, the point ${\mathrm{u}}_{1},$ where the solution went is a stable point (see Figure 3).

Suppose that the following hypothetical collection of data;

When we chose the hypothesis data set (25), one can notice that ${k}_{1} < {k}_{2}$ but ${k}_{2}-{k}_{1} > 1$ which makes the sign of the second eigenvalue of the characteristic equation (10) negative. So, according to Roth-Hwartize criteria, the point ${u}_{1}$ will be completely unstable equilibrium point (see Figure 4 (a). On other hand, This data set gives values for the parameters ${R}_{1} = -31.5$ and ${R}_{2} = 17.0484$ which are represent the coefficients of a polynomial (12) whose number of roots determines the existence of the positive equilibrium point ${u}^{*}$. As ${R}_{1} = -31.5 < 0$ and ${R}_{2} = 17.0484 > 0$, so, according to Lemma 3, the model (3) admit a unique positive equilibrium point. The local stability of positive point is investigate since the condition (15) which gives a sufficient condition for routh howartize theorem is satisfied as 14.1606 = ${x}^{*} > \frac{{k}_{1}}{2} = 7.5$ (see Figure 4 (b)).

Suppose that the following hypothetical collection of data;

This data set gives values for the parameters ${R}_{1} = -31.5$ and ${R}_{2} = 17.0484$ which are represent the coefficients of a polynomial (12) whose number of roots determines the existence of the positive equilibrium point ${u}^{*}$. As ${R}_{1} = 0.1000 > 0$ and ${R}_{2} = -0.0516 < 0$, so, according to Lemma 3-ii of the model (3) admit a unique positive equilibrium point. The local stability of positive point is investigate since the condition (15) which gives a sufficient condition for routh howartize theorem is satisfied as $0.4640 = {\mathrm{x}}^{*} > \frac{{k}_{1}}{2} = 0.25$ (see Figure 5).

6.   Conclusions

In this article, we investigated the fear effect on filter feeding fish from predatory fish and harvesting effect of predatory fish of an aquatic ecological model of kind microcytic aeruginosa-filter feeding fish-predatory fish model with aggregation effect of microcytic aeruginosa. In the begging, to reduce the number of parameters of the model (2), we use non dimensionlization. The positivity and feasible region for the solution of the model (3) have been derived. In addition, we have proved some qualitative analysis by using Routh-Hurwitz stability criterion. In the numerical section, we have presented the approximate solution and their graphical presentation in various Figures including Figure 6 (a), Figure 6 (b), Figure 7 (a), Figure 7 (b) and Figure 8 to Figure 13 respectively. Some graphs of time series and phase portrait are given in this section also. As, we explained the outcomes one at a time by figuring out the solution and displaying the figures, this section demonstrated the validity of the theoretical conclusions obtained by this study. In future, the proposed study will be updated by using the concept of fractional calculus.

Conflict of interest

There does not exist any conflict of interest.

Acknowledgment

Author Kamal Shah and Thabet Abdeljawad would like to thank Prince Sultan University for paying the APC and support through TAS research lab.