Annual Report 2022

1.   Introduction

Marine ecosystems are composed of large and complex networks of energy suppliers and consumers ^[1]. Plankton models are receiving increased attention due to the need to protect resources and the environment ^[2]; their populations are affected by a variety of elements such as illumination and temperature. Conversely, the influence of plankton on the diversity of marine systems is also not negligible. Plankton-fish is the most prominent food chain in the marine system, and their interactions lie at the bottom of the system. These are specific types of interactions in food chain models ^[3]. A number of mathematical models have been constructed on the effects of plankton-fish interactions ^[4,5,6,7]. Dubey et al. ^[8] analyzed the rich dynamical behavior of a plankton-fish model with fear and carry-over effects, including bifurcation behavior after adding a delay. Meanwhile, to combine the fear effect and refuge, Kaur et al. ^[9] analyzed the kinetic behavior of a three-dimensional zooplankton-fish model. Toxic substances may directly or indirectly affect the evolution of species; Raw et al. ^[10] demonstrated the kinetic behavior of a three-dimensional ocean model under the effect of toxicity and nonlinear harvesting. Additional food is also often incorporated into food chain models. For example, Sajan et al. ^[11] analyzed the effect of additional food on the plankton-fish model. However, theoretical research on plankton-fish modeling, especially dynamics and control, is still essential for the broader development of marine ecological modeling.

Interacting species tend to evolve on different time scales. Incorporating slow parameters into a mathematical model can greatly change the dynamics of the model and allow to capture some realistic scenarios. While initially attracting attention owing to the large differences in reaction rates of chemical reactions, many scholars have recently introduced fast and slow systems into ecological models. The slow parameter is interpreted as the ratio between predator mortality and prey growth rate ^[12]. Examples include the use of constant slow parameters and time-varying slow-variable parameters. Sahoo et al. ^[13] considered the oscillatory and instantaneous dynamics of a fast-slow model in which carry-over and fear effects were present. In addition, complete bifurcation analysis of small time-scale parameters and the fast-slow nature of the attractor were investigated when the predator had a scale-dependent functional response function ^[14]. Branches and oscillations of a 1 fast-3 slow dynamical system were investigated when adaptive change occurred in the diet of a predator population ^[15]. Nevertheless, few works have introduced fast and slow parameters to model the dynamics of marine planktonic ecosystems. When the system dimensionality is high, Cortez et al. ^[16] showed a method for reducing the dimension of a model by using the fast and slow theory, and investigated system behavior when predators or prey exhibit fast evolution. Moreover, some literature has suggested that some parameters vary slowly with time, compared to changes in population ^[17,18,19]. However, the slow parameter is provided as a fixed value in most of the literature, and replacing it with a variant form can make the model more realistic.

Ecosystem models have very complex dynamical behaviors that sometimes produce chaos. Most theoretical research has shown that the presence of chaotic ecosystems may affect many critical characteristics such as predictability, species persistence, and biodiversity ^[20,21]. Therefore, it is significant to identify chaotic behavior in ecosystems. There are a wide range of methods for identifying chaos in dynamical systems, including the Lyapunov exponent, fractional form dimension, power spectrum method, Poincaré map, and $0-1$ test ^{[22,23,24,25]}. Hastings et al. ^[26] provided detailed information about chaos theory and its components, such as the definition of chaos, singular attractors, and the Lyapunov exponent. When the spatial distribution of species remains regular, Morozov et al. demonstrated that population oscillations can become chaotic. Moreover, this chaos is the product of a range of periodic doubling bifurcations, where small variations of parameters cause alternation between chaotic and regular dynamic behavior ^[27]. Based on topological horseshoes, Herrera et al. ^[28] proved the existence of sufficient conditions for chaos in a model with impulses. Gakkhar and Sunita investigated the existence of chaos in a food-chain model with nonlinear functional responses ^[29]. Chaotic dynamical behavior of predator-prey systems with phase structure, harvesting gains, and self-cannibalization was analyzed in ^[30]. Kaur et al. ^[31] studied the chaos of the system when additional food was available to both zooplankton and fish. Chaotic behavior of zooplankton-fish in the presence of additional food, fear effects, and modified Holling-IV functional response functions were discussed in detail in ^[11]. In addition, competition and refuge are inevitable in real plankton-fish systems.

Ecological chaos can have many significant effects on population dynamics. Chaos is highly disorganized and unpredictable, which can lead to the extinction of certain species, plankton blooms, or anything unpredictable in nature. Therefore, identifying factors that may enhance or suppress chaos is a challenging problem. Many researchers have proposed and applied chaos control methods, such as sliding mode control (SMC) ^[32], delayed feedback control ^[33], backstepping control ^[34], and impulse control ^[35]. For the chaos control of Chua's circuit, Kocamaz et al. ^[36] discussed the performance of classical sliding-mode, cubic sliding mode, and partial cubic sliding mode, and pointed out that the partial cubic sliding mode controller is more effective. Nonetheless, there is some literature on the application of sliding mode control to suppress chaotic behavior in zooplankton-fish systems. For the discrete time model, Akhtar et al. applied a hybrid control device to control the chaos generated by the Neimark-Sacker bifurcation ^[37]. Bounded feedback control to achieve calming of instability in uncontrolled chaotic systems is discussed in ^[38]. In addition, external signals can also control chaos, such as periodic sinusoidal forces or noise signals ^[31,39]. Kaur et al. ^[31] analyzed the effect of the presence of periodic forces and predation delays in zooplankton to weaken the chaos. However, chaos control of a system in the presence of weak periodic forces in phytoplankton, and gestation delays in intermediate predators is rarely seen.

Dubey et al. ^[8] developed a plankton-fish model that takes into account both the fear effect and the carry-over effect from predators and also analyzed the impact of the delay in the carry-over effect. Motivated by the above, this paper proposes a fast-slow plankton-fish model with slow variable parameters. The chaotic behavior of the model is discussed and several control methods are applied to suppress the chaos of the system. There are several differences compared to ^[8]. First, fish evolve at a slower rate compared to zooplankton, so we considered two time scales. Specifically, fish evolve with a slow time-scale parameter with a correction parameter related to fish density. Second, the soil, nutrients, and space used for phytoplankton growth are limited and competition for phytoplankton mortality is considered. In addition, prey in nature hide in places predators cannot easily find to protect themselves from harm. Incorporating zooplankton refuge into the model can make it more realistic. Finally, with respect to the chaotic control aspect of the system, we use a gestation delay, a seasonal force to calm the chaotic system to periodic solutions, and apply a partial cubic SMC to stabilize the chaos at the internal equilibrium point.

This manuscript is organized as follows. The established plankton-fish model and positivity and boundedness of the solution are in Sections 2 and 3, respectively. In Section 4, we explore the existence of biologically feasible equilibrium, stability conditions, and Hopf bifurcations. In Section 5, several methods are employed to support the existence of chaos in the system. In Section 6, we discuss the control of chaos in the model. Specifically, we employ competition, refuge, gestation delay, weak seasonal force, and sliding mode control to suppress or eliminate chaotic behavior in the system. In Section 7, we discuss the differences between the model established and the existing models. The mechanisms and benefits of different ways to control chaos are also compared. Finally, we arrive at our conclusion in Section 8.

2.   Mathematical model

The model in ^[8] is depicted as

The densities of phytoplankton, zooplankton, and fish are denoted as $P$, $Z$, and $F$, respectively. $K$ is the environmental carrying capacity, $r$ is the maximum intrinsic growth rate, $\beta_2$ is the net gain rate of zooplankton with respect to phytoplankton, and $\beta_{4}$ is that of fish with respect to zooplankton. $\frac{1+cZ}{1+cZ+\varpi F}$ defines the fear and carry-over effects of zooplankton-fish interactions. The other parameters have their own meaning.

Based on the different forms of the slow parameter in fast-slow systems ^[17,18,19], the slow parameter in our system is extended to the form of the variable parameter associated with the variation in population density. We consider a plankton-fish model with variable slow parameters. In this model, the slow-variable parameter $\gamma e^{-\theta F}$ is used as a time-varying parameter for fish, where $\gamma$ is a positive constant less than 1 and $\theta$ denotes the correction parameter. $\theta$ shows that the change rate of fish population density decreases with increasing fish density per unit time. In addition, space for phytoplankton growth, resources, and nutrients are limited, so mortality due to competition are considered. Refuge is a crucial factor in the predator-prey model ^[1,9,48]. Prey in nature protect themselves from harm by hiding in places that are not easily detected by predators, which are called refuge. Considering the refuge of zooplankton from fish predation in the model may be a more realistic scenario. Based on the idea of modeling inter-species interactions ^[13,14], this paper explores how modeling with slow variable parameter scenarios affects the dynamics of the system. Considering the above assumptions, our proposed phytoplankton-fish model is given by the following differential equation,

where $\gamma \in \left(0, 1\right]$, and $\theta$ is a positive constant. The biological meaning and source of all the parameters are listed in Table 1.

3.   Positivity and boundedness

The boundedness of the solution shows that the ecological model is biologically benign. Meanwhile, it also implies that no single species in the system has experienced unintended growth for a extended period of time, and the densities of all species are constrained by limited resources. In this section, positivity and boundedness of solutions in (2) are derived.

Theorem 3.1. Under positive initial value conditions, all solutions of (2) are non-negative. Further, the set $\Gamma = \{(P(t), Z(t), F(t))\in R_+^3:0$ < $P(t)+Z(t)+F(t) \leq \frac{\Lambda}{k}+\varepsilon \ \mathit{\text{for any}} \ \varepsilon > 0 \}$ is the invariant set for all solutions inside the positive octant.

Proof. From (2), we get

$P(t) = P(0)e^{\int_{0}^{t} (r(1-\frac{P(s)}{K})-\frac{\beta_1 Z(s)}{P(s)+a}-qP(s))ds},$

$Z(t) = Z(0)e^{\int_{0}^{t}(\frac{\beta_2 P(s)}{P(s)+a} \frac{1+cZ(s)}{1+cZ(s)+\varpi F(s)}-\frac{\beta_3(1-r_1)F(s)}{(1-r_1)Z(s)+b}-\delta_1)ds },$

$F(t) = F(0) e^{ \int_{0}^{t}\gamma e^{-\theta F(s)}(\frac{\beta_4 (1-r_1)Z(s)}{(1-r_1)Z(s)+b}-\delta_2)ds }.$

Therefore, under the condition that the initial value is positive, all solutions are always positive. That is, all solutions of the plankton-fish system (2) are non-negative. Furthermore, we prove that the above solutions are confined to the octant $\Gamma$.

For any $F\in[0, F(0)]$, $e^{-\theta F}$ is a continuous function in $F\in[0, F(0)]$, according to Weierstrass theorem, $e^{-\theta F}$ is bounded in $[0, F(0)]$, and if the infimum is $m$, then $m \leq e^{-\theta F} \leq 1$. Let $U = P+Z+F$, then the derivative of the solution along Eq (2) is

For any constant $k$, we get

where $k = \min\{\delta_1, m \gamma \delta_2\}$. Based on the definition of environmental carrying capacity, it is obtained that $P\leq K$, then

Noting that $\Lambda = (r+k)K$, after simple deformation, the above formula is obtained

Using the theory of differential inequalities for $U$, there is

where $U(0) = (P(0), Z(0), F(0))$. In addition, $0\leq U\leq \frac{\Lambda}{k}$ when $t\rightarrow \infty$. Thus, all solutions of (2) are constrained to the region

□

4.   Stability analysis

Now, several biologically feasible equilibria of (2) are computed. Then, the stability of these equilibrium is discussed. Finally, we prove the existence of Hopf bifurcations with parameter $r_1$.

4.1.   Existence of equilibrium states

The feasible equilibrium of the ecosystem is what we are interested in. The system (2) has four categories of non-negative equilibrium states as follows:

$(i)$ The zero equilibrium $U_0(0, 0, 0)$ and the predators free-equilibrium state $U_1(K, 0, 0)$ are always present.

$(ii)$ The fish free-equilibrium point $U_2(P_2, Z_2, 0) = (\frac{\delta_1 a}{\beta_2 -\delta_1}, \frac{\beta_2 a (Kr(\beta_2-\delta_1)-(r+Kq)\delta_1 a)}{K\beta_1(\beta_2 -\delta_1)^2}, 0)$ exists when $\beta_2 > \delta_1$ and $Kr(\beta_2-\delta_1) > (r+Kq)\delta_1 a$.

$(iii)$ The internal equilibrium $U_*(P_*, Z_*, F_*)$ exists when conditions $0 < r_1 < 1-\frac{\beta_1 \delta_2 b}{(\beta_4-\delta_2)ar}$, $\delta_2 < \beta_4$, and $\delta_1 a < (\delta_2-\delta_1)P_*$ hold.

According to Eq (2), it is easy to obtain $Z_* = \frac{\delta_2 b}{(\beta_4 -\delta_2)(1-r_1)}$, as $r_1 < 1$ and $\delta_2 < \beta_4$. Substituting $Z_*$ into the first equation of system (2), $P_*$ is the positive root of the below function.

which has a positive solution, when $f(0) < 0$, i.e., $0 < r_1 < 1-\frac{\beta_1 \delta_2 b}{(\beta_4-\delta_2)ar}$.

Furthermore, $F_*$ is the positive root of function (3),

where

$A_2 = \varpi \beta_3(1-r_1),$

$A_1 = \beta_3(1-r_1)(1+cZ_*)+\delta_1 \varpi ((1-r_1)Z_*+b),$

$A_0 = (1+cZ_*)((1-r_1)Z_*+b)(\delta_1-\frac{\delta_2 P_*}{P_*+a}).$

Thus, under the condition $g(0) < 0$, i.e., $\delta_1 < \frac{\delta_2 P_*}{P_*+a}$, (3) has a positive solution.

Remark 1. This implies that fish survival becomes possible when zooplankton mortality is lower than their population response to phytoplankton.

Biologically, in the absence of these two predators, phytoplankton biomass survival is higher than in their presence (see Table 2).

Remark 2. As can be seen from Table 3, in the presence of fish, zooplankton biomass is maintained at a low level and phytoplankton biomass is high. This coincides with realistic food chain models.

4.2.   Local stability

In the previous subsection, we explored the feasible equilibrium points of system (2). Now, the local stability of all these feasible equilibrium states is examined. We first compute the Jacobian matrix for each equilibrium and then establish their stability based on the matrix eigenvalues.

Theorem 4.1. The local stability at the four types of equilibrium points in system (2) is as follows:

$(i)$ The zero equilibrium $U_0$ is unstable.

$(ii)$ The predators' free equilibrium $U_1$ is locally asymptotically stable (LAS) if $K < \frac{\delta_1 a}{\beta_2 -\delta_1}$.

$(iii)$ The equilibrium $U_2$ is LAS if $\frac{\beta_4(1-r_1)Z_2}{(1-r_1)Z_2+b} < \delta_2$ and $r < \frac{2r}{K}P_2 -2qP_2-\frac{\beta_1 a Z_2}{(P_2+a)^2}$.

$(iv)$ The equilibrium $U_*$ is LAS if $\eta_1 \eta_2 > \eta_0$, $\eta_0 > 0$ and $\eta_1 > 0$, where $\eta_0$, $\eta_1$ and $\eta_2$ are provided in the proof.

Proof. $(i)$ Since $r$ is the eigenvalues of $\left.M\right|_{U_0}$ (the Jacobian matrix at $U_0$), the zero equilibrium is unstable.

$(ii)$ The eigenvalues of $\left. M\right|_{U_1}$ are $-(r+2qK)$, $-\delta_2 \gamma$ and $\frac{K\beta_2}{K+a}-\delta_1$, so $U_1$ is LAS when $\frac{K\beta_2}{K+a} < \delta_1$, otherwise a saddle point.

$(iii)$ The Jacobian matrix $M$ at $U_2$ is

where

$a_{11} = r-\frac{2r}{K}P_2-2qP_2-\frac{\beta_1 a Z_2}{(P_2 +a)^2}, \ a_{12} = -\frac{\delta_1 \beta_1}{\beta_2} < 0,$

$a_{23} = -\frac{\beta_2 P_2 Z_2}{P_2 +a} \frac{\varpi}{1+cZ_2}-\frac{\beta_3 (1-r_1)Z_2}{(1-r_1)Z_2+b},$

$a_{21} = \frac{\beta_2 a Z_2}{(P_2+a)^2} > 0, \ a_{33} = \frac{\beta_4 \gamma (1-r_1)Z_2}{(1-r_1)Z_2+b}-\delta_2\gamma.$

One eigenvalue of $\left. M\right|_{U_2}$ is $a_{33}$, and the other two eigenvalues are determined by Eq (4),

using Routh-Hurwitz criterion, if $a_{11} < 0$, each root of the Eq (4) has a negative real part, which implies $U_2$ is locally asymptotically stable when conditions $r < \frac{2r}{K}P_2-2qP_2-\frac{\beta_1 a Z_2}{(P_2 +a)^2}$ and $\frac{\beta_4 (1-r_1)Z_2}{(1-r_1)Z_2+b} < \delta_2$ hold.

$(iv)$ $M\left|_{U_*}\right.$ is calculated as

where

$b_{11} = -\frac{\beta_1 a Z_*}{(P_* +a)^2}-2qP_*-\frac{2r}{K}P_*+r\ , \ b_{12} = -\frac{\beta_{1}\ P_*}{P_*+a},$

$b_{21} = \frac{\beta_2 a Z_*}{(P_*+a)^2} \frac{1+cZ_*}{1+cZ_*+\varpi F_*}\ , \ b_{23} = -\frac{\beta_2 P_* Z_*}{P_*+a}\frac{\varpi(1+cZ_*)}{(1+cZ_*+\varpi F_*)^2}-\frac{\beta_3(1-r_1)Z_*}{(1-r_1)Z_*+b}\ ,$

$b_{22} = \frac{\beta_2 P_*}{(P_*+a)(1+cZ_*+\varpi F_*)}(1+cZ_*+\frac{c\varpi Z_* F_*}{1+cZ_* +\varpi F_*})-\delta_1-\frac{\beta_3 (1-r_1)bF_*}{((1-r_1)Z_*+b)^2} \ ,$

$b_{32} = \gamma e^{-\theta F_*}\frac{\beta_4 b(1-r_1)F_*}{((1-r_1)Z_*+b)^2}.$

After a simple calculation, the characteristic equation of $M\left.\right|_{U_*}$ is

where $\eta_1 = -(b_{11}+b_{22}), \ \eta_0 = b_{11}b_{32}b_{23}, \ $ $\eta_2 = -b_{12} b_{21}-b_{32}b_{23}+ b_{11} b_{22}.\ $

According to Routh-Hurwitz criterion, the roots of (5) have negative real parts when $\eta_1\eta_2 > \eta_0$, $\eta_0 > 0$, and $\eta_1 > 0$. □

Therefore, if these conditions in Theorem 4.1 hold, the corresponding feasible equilibrium point is locally asymptotically stable.

4.3.   Global stability

Global asymptotic stability (GAS) in dynamic systems is often proved by constructing Lyapunov functions. In this subsection, on the basis of the principle of Lyapunov function, the global stability condition of $U_*$ is presented.

Theorem 4.2. The internal equilibrium $U_*(P_*, Z_*, F_*)$ is GAS when conditions

hold, where $M = \frac{\Lambda}{k}$, $k = \min \{\delta_1, m \gamma \delta_2 \}$. $P_m$ and $Z_m$ are provided in the proof.

Proof. According to Theorem 3.1, there exist $T$ and $M$ such that $P(t), Z(t), F(t) < M$ for $t > T$. Based on system (2), we obtain

which implies $P\geq \frac{K}{r+Kq}(r-\frac{\beta_1 M}{a}) = :P_m$. Likewise, we get

which implies $Z \geq \frac{\beta_2 P_m}{(M+a)(1+(c+\varpi) M)}-\frac{\beta_3(1-r_1)M}{b}-\delta_1 = :Z_m$.

Then, the Lyapunov function $V(P, Z, F)$ with respect to $U_*$ is defined by

Thus,

i.e.,

where $K_1 = \frac{r}{q_1}$, and $q_1 = \frac{r}{K}+q$. Because of the size relationship $\beta_3 \geq \beta_4$, $\beta_1\geq \beta_2$ and $\gamma < 1$, therefore

It is noted that if

then $\frac{dV}{dt} < 0$. Following the Lyapunov-Lasalle's invariance principle, the equilibrium $U_*$ can be said to be the GAS. □

Biologically, when the system starts at any initial population level, the solution trajectory of the system eventually reaches the basin where that internal equilibrium is located. That is, if the internal equilibrium is GAS, then no population can become extinct.

4.4.   Bifurcation analysis

In this subsection, we investigate how parameter variations impact the dynamics of (2), i.e., Hopf bifurcation. Through the bifurcation, the stability of the system can be switched between different states. Choosing the refuge rate $r_1$ as the bifurcation parameter, while keeping the other parameters constant, we demonstrate that the bifurcation occurs at $U_*$.

Theorem 4.3. For the occurrence of Hopf bifurcation around $U_*$ at $r_1 = r_1^*$, the bifurcation conditions are

$(1)\; \eta_i\left(r_1^ *\right) \geq 0, \ i = 1, 0, \ \eta_1\left(r_1^ *\right) \eta_2\left(r_1^ *\right)- \eta_0\left(r_1^ *\right) = 0,$

$(2)\; \left(\eta_1(r_1^*)\eta_2(r_1^*)\right)' \ne \left(\eta_0(r_1^*)\right)'.$

Proof. Holding other parameters constant, we consider the Hopf bifurcation of the system (2) for $r_1$. First, at the threshold $r_1 = r_1^*$, one has

Thus, the characteristic equation

must have the following form at $r_1 = r_1^*$,

the Eq (7) clearly has roots $-\eta_1(r_1^ *)$, and $\pm i\sqrt{\eta_2(r_1^*)}$. To compute the transversality condition, we assume that the roots $\lambda_i \, (i = 1, 2, 3)$ are functions of $r_1$ such that $\lambda_i(r_1) = u(r_1)\pm iv(r_1)$, $(i = 1, 2)$ and $\lambda_3 = -\eta_1(r_1)$. Substituting the value of $\lambda_i \, (i = 1, 2)$ into Eq (6), yields

Then, deriving the derivative with respect to the parameter $r_1$, we get

where

$L_1 = 3u^2+2\eta_1u-3v^2+\eta_2,$

$L_2 = 6uv+2\eta_1v,$

$L_3 = \eta_1 ' \left[u^2-v^2 \right]+\eta_2 'u+\eta_0 ',$

$L_4 = 2\eta_1' uv+\eta_2 ' v.$

Taking $u(r_1^*) = 0$, and $v(r_1^*) = \sqrt{\eta_2(r_1^*)}$, one obtains

Solving (8) for $u'(r_1^*)$, under the assumptions, yields

□

Thus, the transversality condition holds, leading to a Hopf bifurcation at $r_1 = r_1^*$.

Due to the complexity of system (2), the theoretical proof of the stability of the Hopf bifurcation still needs to be further investigated. However, a numerical approach using MatCont software is described in detail in ^[40]. Inspired by this, we plot the limit cycle generated by the Hopf bifurcation, and the solution curves of $P$, $Z$, and $F$ with respect to the bifurcation parameter $r_1$, as shown in Figure 1. The results show that the system (2) varies from unstable limit cycle to stable equilibrium with the parameter $r_1$, and thus a subcritical Hopf bifurcation occurs around $U_*(15.8659, 17.7767, 18.9283)$ at $r_1^*\approx0.45648$.

Remark 3. Similarly, the system (2) may undergo Hopf bifurcation for the parameter $q$.

5.   Chaotic dynamics of model (2)

Through the previous subsection, we know the conditions for the occurrence of Hopf bifurcation in system (2). Bifurcation and chaos are the most important and fundamental features in nonlinear dynamical systems. In this section, the chaotic dynamics of (2) are discussed. We apply several methods to support the existence of chaos in the system (2).

5.1.   Dissipativity condition ^[22]

Dissipation is a general condition for system (2) to manifest chaotic behavior. Thus, let $W$ be the volume element of the system trajectory flow. To verify the dissipation condition, we consider the function

A simple algebraic calculation at $U_*$, yields

when condition

holds. As time tends to infinity, the volume of the trajectory that contains the system (2) will shrink to zero with an exponential rate. The asymptotic motion of the attractor in the three-dimensional phase space is shown in Figure 2.

5.2.   Sensitivity of initial conditions

Numerous studies have defined chaos as a system that is sensitive to initial conditions, which is one of the most obvious signs of chaotic behavior. In Figure 3, the phase diagram displays the chaotic behavior of the state variables. In addition, time series plots with two sets of initial values $[(10, 15, 10)$ and $(10.1, 15, 10)]$ are shown in Figure 3. From Figure 3, we notice that the time series plots of the two sets of initial values follow almost the same trajectory before $200$ iterations but, as time proceeds, the two sets of trajectories deviate from each other, exhibiting chaotic behavior of the system.

5.3.   Analysis of Lyapunov exponents

The Lyaponuv exponent expresses the average exponential divergence of neighboring trajectories in phase space and is one of the most effective methods for detecting chaos in dynamical systems. For system (2), three Lyaponuv exponents are estimated. The obtained Lyapunov spectrum of (2) is shown in Figure 4(a). The exponent $\lambda_3 = -0.8141$ indicates an increase in the contraction of the attractor in the phase diagram, where the maximum exponent $\lambda_1 = 0.0192$ shows an increase in the unfolding of the attractor. $\lambda_2 = -0.0013$ characterizes the critical nature of the attractor. The complexity of this singular attractor is expressed in Kaplan-Yorke dimension as

Consequently, it can be said that the system (2) exhibits chaos characterized by fractional order dimensions.

Remark 4. Compared to ^[8], Figure 4(a) has a larger maximum Lyapunov exponent, indicating a higher complexity of system (2) ^[51].

5.4.   Alternative characterization of chaotic dynamics

There are several ways to characterize chaos, including Poincaré map, power spectral density (PSD), and $0-1$ test. In this subsection, we will give images of these characterization methods to further underpin the view of the chaos behavior of system (2), with parameter values using Table 1.

To begin with, the $0-1$ test algorithm is employed to check the chaotic behavior in system (2). The method is applicable to both discrete and continuous systems, and it does not concern the reconstruction of phase space. If a time series $\psi(k), (k = 1, 2, 3, ...)$ denotes the data set gained in the plankton-fish model (2), then the translation variables $p$, $q$ are defined as ^[24,25]

where

for $n = 1, 2, 3, ...$. In general, Brownian (unbounded) trajectories in the $(p, q)$ plane imply that the dynamics are chaotic, whereas bounded trajectories imply regular dynamics.

Further, a random selection of $c$ from $(0, \pi)$ allows one to compute the mean-square displacement $M_c$, defined by

Adopting $D_c$, which has better convergence than $M_c$. $D_c$ is obtained from the following equation

where $E\psi = \lim\limits_{N\rightarrow \infty} \frac{1}{N} \sum_{j = 1}^{N} \psi(j)$, The asymptotic growth rate $\tilde{D}_c$ of $D_c$ is then calculated using the correlation method. Eventually, the values of $K$ corresponding to the different parameter values are obtained from the median value of $\tilde{D}_c$. Values of $K$ closer to $1$ indicate chaotic dynamics and closer to $0$ indicate no chaos.

The $p-q$ graph is shown in Figure 5(c), and the diffusion behavior of $p$ and $q$ shows that $M_c$ is linear with respect to time, indicating that the system is in chaos.

The Poincaré map is defined by the motion of the trajectory line in phase space. The main approach is to choose the proper section for the trajectory, simplify the evolution of the trajectory in space to the motion of the intersection of the trajectory line with this section, and succinctly determine the state of the motion based on the trajectory passing over the section. First, we compute a series of discrete points of the trajectory line, then choose a suitable plane, compute the intersection of the rail line with the plane, and finally plot the intersection to obtain the Poincaré cross section. Specifically, the plane we choose is the $P-F$ plane with $Z = 19$. The map drawn using the parameters of Table 1 is shown in Figure 5(a), where dense points with fractal structure appear in the section, indicating the chaotic state of system (2). PSD is also an effective approach to detect the chaotic behavior of a system. The PSD curves in Figure 5(b) have no obvious peaks and are in the form of thick bands, signaling that (2) is chaotic.

6.   Chaos control

Despite the fact that it is helpful to detect chaos in a system, it is necessary to control chaotic systems in order to improve their effectiveness. Various active control methods and external signals to control chaos have been proposed. In this section, how to effectively transform chaos into stable oscillations or steady states is explored. Adjusting the system parameter values within the feasible interval is one of the most commonly used methods. We first exploit competition and refuge as chaos control parameters. Furthermore, the following techniques are applied to maintain the stability of the system: (ⅰ) Gestation delay control, (ⅱ) weak periodic force control, and (ⅲ) sliding mode control (SMC).

6.1.   Chaotic control through refuge and competition of two plankton

Chaos has become an important issue, with numerous studies on chaos and chaos control techniques. Various active control methods and external signals to control chaos have been proposed. It is common to control chaos by adjusting the values of system parameters. In this subsection, we exploit the mortality due to competitive of phytoplankton ($q$) and the refuge of zooplankton ($r_1$) as chaos control parameters.

Remark 5. In Figure 6, the bifurcation maps of the three state variables expose the rich dynamical behavior, including chaos, in the system (2). The maximum Lyapunov exponent (Figure 4(b)) is also computed to verify the switching between chaos and steady state in the bifurcation diagram.

Remark 6. The bifurcation diagram in Figure 7 reveals that the refuge parameter $r_1$ influences the dynamical behavior of the system. The maximum Lyapunov exponential spectrum (Figure 7(a)) and the $0 - 1$ test plot (Figure 7(b)) imply the accuracy of the bifurcation diagram.

Figure 8 illustrates the chaotic system shifting to a stable system as the intra-species mortality rate $(q)$ of phytoplankton due to competition increases, with other parameters from Table 1. Specifically, Figure 8(a) shows that the system remains in a chaotic state at $q = 0.015$. The dynamical behavior changes from chaotic to bi-periodic orbitals (Figure 8(b)) at $q = 0.03$ and to a stable limiting cycle at $q = 0.06$ (Figure 8(c)). With $q$ still increasing, the system switches to the stable focus at $q = 0.08$ (Figure 8(d)). In ecology, when phytoplankton overpopulation occurs, the chaos of the system can be effectively weakened by reducing the amount of nutrients or living space available for phytoplankton growth, which is consistent with reality.

In addition, Figure 9 shows that as the zooplankton refuge rate $r_1$ changes, the system switches from chaotic to periodic orbits, which can eventually be converted to a stable focus. We observe through simulations that the chaotic dynamics (2) remain chaotic at $r_1 = 0.18$ (Figure 9(a)), switches to triple period at $r_1 = 0.27$ (Figure 9(b)), and switches to single period at $r_1 = 0.383$ (Figure 9(c)). Periodic oscillations of the Hopf bifurcation occur at $r_1 = 0.455$ (Figure 9(d)). As $r_1$ continues to increase to $0.47$, the system switches to a stable focus (Figure 9(e)). Biologically, when there are fewer zooplankton and too many of their natural enemies, the loss of zooplankton can be effectively slowed down by increasing the number of refuges. This is equally effective in suppressing chaotic behavior.

6.2.   Chaos control through gestation delay

Delay is an important factor in water ecology models, where the growth of zooplankton after ingestion of phytoplankton does not occur instantaneously and there is a discrete time delay representing the gestation period of the zooplankton. A time delay in the adaptation of intermediate predators to variable environments can make the model inherently more realistic. Meanwhile, oscillatory dynamics with a specific period and amplitude can be captured by varying the value of the delay. In this subsection, we explore the control of chaos by the gestation delay in the system. The delay $\tau$ is incorporated into (2), and we obtain

Remark 7. With increasing time delay $\tau$, the dynamics appear chaotic, with a periodic solution and a stabilizing focus ^[43]. In this paper, we apply the gestation delay to stabilize the chaos of system (2).

Remark 8. Figure 10 bifurcation plot of gestation delay shows that increasing gestation delay in some ranges can suppress chaos.

We demonstrate numerically that gestation delay can stabilize chaos in the plankton-fish system. Figure 11 shows that when the time delay is small, the model exhibits chaos at $\tau = 0.3$. When increasing to $0.68$, the chaotic dynamics transforms into a multi-periodic one, which exhibits a second-order periodic solution at $\tau = 1.95$, and continues to increase up to $\tau = 2.15$ exhibiting a stable limit cycle.

Remark 9. It is implied that, all other conditions being constant, the longer it takes to convert to zooplankton growth after swallowing the plankton, the less phytoplankton is reduced due to predation for a period of time, which can make the system more stable.

6.3.   Chaos control via weak periodic forces

It is important to study cyclically externally driven ecosystems since the various populations are almost always in cyclically changing environments. In particular, many of the seasonal mechanisms are inherent in real predator-prey models. In this subsection, we will consider how chaotic dynamics change when there are periodic fluctuations in phytoplankton growth. A weakly periodic force function $u = A_1+A_2\sin (\omega t)$ is employed, where $A_2$ and $\omega$ express the amplitude and frequency of the external weak periodic force, respectively, and $A_1$ is a constant deviation.

The chaotic state can be transformed into a stable state by selecting appropriate parameter values for the weakly periodic force function. We assume that the appropriate parameter sets are $[S_1]:A_1 = 0.6, A_2 = 0.5$, $[S_2]:A_1 = 0.67, A_2 = 2$ and $[S_3]:A_1 = 1.3, A_2 = 1$, with other parameters from Table 1 and $\omega = 5.5$. Numerical methods show that weak periodic forces can control chaos. Figure 12(a) indicates the presence of multi-periodic orbits in the parameter set $[S_1]$. Figure 12(b) exhibits the double-limit cycle behavior at $[S_2]$, and at $[S_3]$ the system is stabilized to the limit cycle, as shown in Figure 12(c).

Remark 10. Systems can become chaotic due to seasonal or external cyclical forces ^[45,46]. Furthermore, Rajinder et al. presented the effect of seasonal forcing on the control of chaos in marine populations. They discussed that a given chaotic dynamical system can be controlled as stable, single-periodic and double-periodic orbits with different seasonal parameters ^[31].

6.4.   Chaos control by sliding mode control

Sliding model control is robust to external perturbations and parameter uncertainties, and has been widely used in a multitude of applications. When phytoplankton reproduce rapidly, algal blooms can be harmful to the marine ecosystem due to the deterioration of seawater quality, hypoxia in the water, and death of fish, even possibly affecting human health. In this subsection, we utilize the population density of zooplankton as a feedback mechanism to control phytoplankton growth rates. For example, when phytoplankton overpopulation occurs, measures that can be taken include harvesting phytoplankton, using phytoplankton-killing drugs, or increasing zooplankton populations. We apply a partial cubic sliding mode controller that performs well in chaotic control. By adding the control function, the plankton-fish model defined by (2) can be written as

where $u_1$ and $u_2$ are the control function, and it will be determined by the theory of SMC.

Remark 11. Sliding model control has been widely used in areas such as power systems, network communications, circuits, and predator-prey systems. However, the sliding model control of chaotic states is rarely seen in marine populations.

The purpose is to control chaos and stabilize the system to the point of internal equilibrium $U_*(P_*, Z_*, F_*)$. The error functions are $e_1 = P-P_*$, $e_2 = Z-Z_*$, $e_3 = F-F_*$, and then the error dynamics of system (11) are stated as:

In the defined error dynamics (12), if $e_1 = e_2 = 0$, then $\dot{e_3} \leq 0$, which means that the error $e_3$ will approach zero ($e_3 \rightarrow 0$) as the time tends to infinity $(t\rightarrow \infty)$. Hence, the following SMC signals are appropriate

where $k_1$ and $k_2$ are chosen positive constants.

$s\dot{s}\leq 0$ is the reachability condition design for sliding mode. In the initial stage of the system, when the dynamic error is greater than $1$, the cubic SMC function has better control. When the dynamic error is less than $1$, the cubic version loses its performance and can be switched to the classical SMC solution, so the partially cubic SMC exhibits better performance. The partial cubic feedback control function can be taken as ^[36]

where $k_3$ and $k_4$ are constant parameters, and a larger value of $k_3$ shortens the time to reach the sliding surface but causes chattering. Conversely, a smaller value of $k_4$ will reduce chattering but increase the time. This means that it is neccessary to select the appropriate $k_3$ and $k_4$.

Maintaining the system on the sliding surface is the primary goal of the control signal. Here, the control (13) guarantees that the system (11) is on the sliding mode surface $s = 0$.

Noting that $\text{sign}(s^3) = \text{sign(s)}$, so there is

We take into account the Lyapunov function $V = \frac{1}{2} s^2$, and get

Since $s \ast \text{sign}(s) = \left| s \right| \geq 0 \nonumber$, then $\dot{V} \leq 0$ as long as $k_3 > 0$ and $k_4 > 0$. This assures that $s$ is asymptotically stabilized to the zero error point.

Ultimately, the system (11) is controlled by the control function in Eq (13). After repeated adjustments, when we select the gain parameters as $k_1 = 0.6$, $k_2 = 0.5$, $k_3 = 0.05$ and $k_4 = 0.1$, the added controller at $t = 1000$ can control the chaotic orbit to the stabilized internal equilibrium point $U_*$ (Figure 13(a)). The error time series graphs are shown in Figure 13(b). The above results illustrate that the partial cubic SMC can stabilize the chaos to the internal equilibrium point in a relatively short time.

7.   Discussion

According to the 3D water ecology model ^[8], phytoplankton are the prey, zooplankton are the intermediate predators, and planktivorous fish are the top predators. It assumes that the growth of the three groups takes place at the same time scale. In the literature, many researchers have dealt with planktivorous fish models with various ecological factors and have obtained interesting results. Nevertheless, only a few studies have independently explored the effects of phytoplankton competition and zooplankton refuge on the system where fear and carry-over effects are present. Moreover, it was shown that the Hopf bifurcation threshold is related to the value of the slow parameter in the fast-slow system ^[14]. With the above motivation, we explored the dynamic behavior of the fast-slow system when variable slow parameters are present in the system. Meanwhile, our plots of the $0-1$ chaos test, PSD, and Poincaré map support this view of the existence of chaotic behavior of the system (2). In our model, the value of the slow variable parameter does not change the feasible equilibrium point of the system, but it changes the stability of the equilibrium state.

In terms of chaos control in the system, combined with biological significance, we use four control methods to calm chaos. First, by human control, dispersing or aggregating phytoplankton thereby alters the competitive death of phytoplankton. The calming of chaos is achieved by artificially increasing or decreasing the refuge of zooplankton. Second, growth of zooplankton after ingestion of phytoplankton does not occur instantaneously and there is a gestation delay. The gestation delay is generally different for different species of zooplankton, and chaos can be controlled by artificially putting in zooplankton with different gestation periods. In addition, biomes in ecosystems are almost always in cyclically changing environments, and different classes of phytoplankton have different abilities to perceive seasonality. Artificially increasing or harvesting phytoplankton of appropriate seasonality achieves the control of chaos in the system by different seasonal parameters. Finally, when phytoplankton overpopulation or even bloom occurs, the growth of phytoplankton is suppressed by artificially harvesting the phytoplankton or by increasing its natural enemies, so we propose SMC.

The results show that the gestational delay of zooplankton can convert the chaotic behavior of system (2) into multi-cycle, double-cycle, and single-cycle, as shown in Figure 11. From Figure 12, when periodic fluctuations in phytoplankton growth occur, the chaotic behavior manifests itself as multi-periodic, bi-periodic, and single-periodic at $[S_1]$, $[S_2]$, and $[S_3]$, respectively. Meanwhile, chaos can be stabilized to periodic solutions and a stable focus by selecting the competition and refuge in a suitable range. Moreover, a partial cubic SMC is applied to suppress chaos in the plankton-fish system. Figure 13(a) shows that by activating the controller at $t = 1000$, the original system is stabilized to the internal equilibrium point $U_*$. The designed controller stabilizes the chaotic system to the internal equilibrium point of the system, and this result supports that a sliding mode controller can control chaos in marine biological systems.

8.   Conclusions

In this article, we consider a class of ecological models with slow variable parameters to study the chaos of plankton dynamics. We find that chaos can be controlled by different factors/methods, including competition, refuge, gestation delay, weak seasonal forces, and sliding mode control. Analytically, the positivity and boundedness of the solutions are proved, while the existence and stability conditions of each feasible equilibrium point are discussed. We have also supported the chaotic behavior of the proposed model using several methods. In the presence of the slow variable parameter case, we observe that chaos can be stabilized to a periodic solution by suitable parameter (competition and refuge) ranges, weak periodic forces, and gestation delay, and calmed to an internal equilibrium point by sliding mode control. Therefore, it can be concluded that the current biological model exhibits several different intervals of internal parameters and external factors, demonstrating rich dynamical behavior and multiple chaotic controls in the plankton-fish model.

Use of AI tools declaration

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

Conflict of interest

The authors declare no conflict of interest.