Dynamics of an eco-epidemiological model with toxicity, treatment, time-varying incubation

1.   Introduction

Since the pioneering work studied by Lotka ^[1] and Volterra ^[2], predator-prey interaction models have become one of the most important and fascinating research areas in ecology. In recent years, researchers have expanded the classic predator-prey model into the eco-epidemiology. They focused on studying predator-prey systems with infectious diseases, further analyzing the dynamic relationship between predator and prey ^{[3,4,5,6,7,8,9,10]}. This interdisciplinary research approach helps us to better understand the impact of disease transmission on the structure and function of food webs. These research findings not only enrich ecological theory but also provide crucial insights for formulating more effective ecological management strategies. Beltrami and Carroll ^[3] established an eco-epidemiological model based on predator-prey, in which the prey population appears to be infected by a virus, forming an infected population. They discovered that such a system has been disrupted by a small amount of infectious factors, otherwise a stable structure would be observed. Meng et al. ^[5] showed that in the predator-prey system, the predator population can survive and spread the disease among the prey species. Recently, Sk and Pal ^[6] studied a predator-prey model with infected prey, and found that high levels of fear and low levels of refuge behavior can eliminate the disease from the system in deterministic and stochastic environments.

Plankton, especially phytoplankton occupying the primary trophic level, forms the foundation of all aquatic food chains. These microscopic algae thrive in marine and freshwater environments, typically existing in low concentrations but capable of proliferating extensively on the water surface, resulting in dense cell aggregations known as algal blooms. While the majority of algae in aquatic environments are beneficial, a small number possess the ability to produce toxins. Viral infections within the plankton community can influence the changes of algal blooms, leading to alterations in the behavior and other aspects of aquatic and marine ecosystems. Based on the concept of plankton diseases, Kermack and McKendrick ^[11] proposed the classic SIR model to explore the impact of infections on populations. In recent years, many scholars have conducted extensive ecological epidemiological models to study the different ecological scenarios where populations are affected by external toxicity, disease transmission, prey refuge, Allee effect, fear effect, and other factors ^{[12,13,14,15]}. Among these studies, phytoplankton-zooplankton systems play a critical role in marine and freshwater ecosystems. Scholars have found that these toxic substances have a significant impact on the growth of zooplankton, and further affect the structure and function of the entire plankton community ^{[15,16,17,18,19]}. Chattopadhyay et al. ^[16] first proposed and studied the toxin-producing phytoplankton-zooplankton system. They suggested that toxin-producing phytoplankton (TPP) exert a biocontrol effect by reducing the grazing pressure exerted by zooplankton, thereby inhibiting zooplankton reproduction. Similarly, Gakkhar and Singh ^[18] proposed an eco-epidemiological time-delay model involving viral infection, TPP, and zooplankton system. Their results demonstrated that viral infections play a crucial role in controlling complex dynamic behaviors, including chaos. Additionally, variations in the toxin release rate of TPP significantly influence the emergence of chaos in the time-delay model. Huang et al. ^[19] studied aquatic toxicity through a simple mathematical model and concluded that high toxin levels can lead to population extinction. Srivastava and Thakur ^[15] proposed an eco-epidemiological model of an aquatic dynamical system:

where $S$, $I$, and $P$ represent the densities of susceptible prey, infected prey, and predator at time $t$, respectively. The parameter $\theta$ represents the toxin release rate of susceptible and infected prey. The significance of other parameters in system (1.1) can be found in reference ^[15]. They analyzed the stability of both spatial and nonspatial models, and found that the density of predator is affected by virus-induced prey and toxicity.

In the coexistence of disease and treatment, controlling the spread of disease is one of the primary objectives in the eco-epidemiological models ^{[5,14,20,21,22,23]}. Ghosh et al. ^[14] investigated an eco-epidemiological model incorporating the effects of fear, treatment, and cooperative hunting:

where $\frac{\rho I}{\sigma+I}$ represents the treatment function, $\rho > 0$ represents the maximum medical resources available for treatment, and $\frac{1}{\sigma} > 0$ represents the saturation factor measuring the impact of delayed treatment on infected individuals. The meanings of other parameters can be found in reference ^[14]. The authors found that the fear of prey for predator, cooperative hunting, and treating infected prey with memory effect play a crucial role in preserving biodiversity.

In an ecosystem with infectious disease, time delay is an important factor to describe the dynamic characteristics of the system. In recent years, many researchers have focused on delay-induced infection models to explore the combined effects of the infection process and time delay on population dynamic ^{[18,24,25,26,27,28,29,30]}, where the delays are typically constant. However, in a real system, time delays are often not fixed and dynamically change over time, making time-varying delays more reflective of the complexity of a natural environment. Consequently, time-varying delays have been receiving increasing attention in predator-prey models, infectious disease transmission models, and broader ecological dynamics research. In recent years, researchers have explored how time-varying delays affect the stability and periodic behavior of systems through theoretical analysis and numerical simulations ^{[31,32,33,34,35,36,37]}. Research on this phenomenon not only helps to understand the fundamental mechanisms of ecosystem and epidemic dynamics, but also provides important theoretical support for devising effective management strategies. Wu et al. ^[35] developed a brucellosis transmission model incorporating seasonal alternation, density-dependent growth, stage structure, maturation delays, and time-varying latency periods. The numerical results showed that if the effects of time-varying latency or maturation delay are ignored, then the transmission of brucellosis would be overestimated (or underestimated). Liu and Meng ^[36] established an eco-epidemiological predator-prey system with time-varying delay, and obtained the permanence and global asymptotic stability of the time-varying delay system.

Motivated by the above work, we will consider an eco-epidemiological model with toxicity, treatment and time-varying delay to analyze the interaction between toxicant prey and predators. This is a novel attempt to consider the effect of time-varying delay on prey and predators in a plankton system. Under the viral infection with and without time delay in a plankton system, investigating such complex dynamics are significant implications for both research on plankton population and ecological modeling.

The organization of this paper is as follows. In Section 2, we perform the formulation of this model. In Section 3, we discuss the positivity and boundedness of the solution of the non-delayed system. Furthermore, we derive the conditions for the existence and local stability of five equilibriums of the system without delay. In Section 4, we focus on investigating the existence of Hopf bifurcation. Additionally, we apply the method of multiple scales (MMS) to the delay differential system. Building on this, we propose a control strategy that involves introducing time-varying perturbation to the delay, which aims to suppress oscillation. Moreover, the global stability of the positive equilibrium for model (2.1) is analyzed through the construction of a suitable Lyapunov function. In Section 5, we present some numerical simulations to validate the theoretical results obtained in the previous sections. We summarize our paper with some discussions in the last section.

2.   Model formulation

Motivated by the above papers, we take Noctiluca scintillan as the prey population and Paracalanu as the predator population. Results have shown that Noctiluca scintillan, a phytoplankton, produces toxins during its breeding season and can have adverse effect on zooplankton ^[17]. We give the following assumptions:

(ⅰ) Assuming that only susceptible prey has the ability to reproduce, the disease only spreads within the prey population, and the disease is not inherited. When susceptible phytoplankton come into contact with infected phytoplankton, they also become infected. Therefore, the virus transmission is supposed to be Susceptible-Infected-Susceptible ($SIS$) type. However, the infected phytoplankton population still continues to grow at the same carrying capacity $K$.

(ⅱ) The zooplankton consumes susceptible and infected phytoplankton in regard to Holling Ⅱ functional response. In addition, we hypothesize that the additional mortality of zooplankton is due to the phytoplankton liberating toxin instantaneously with rate $\theta$. TPP populations do not release toxic chemicals unless there is a dense zooplankton population nearby ^[38].

(ⅲ) In most eco-epidemiological models, researchers assume that infected prey recovers at a linear rate, and they do not consider the effect of treatment. In the real world, the process of treatment recovery can be complex. Thus, we consider a nonlinear treatment function, which may be more realistic.

(ⅳ) Due to factors such as climate, environment, and individual differences, the incubation period of the virus may vary over time. Therefore, considering the time-varying incubation delay may be more in line with the actual situation.

Thus, a time-varying delay eco-epidemiological model with toxicity and treatment is given as follows:

where $S(t)$, $I(t)$, and $P(t)$ represent the densities of the susceptible phytoplankton, infected phytoplankton, and the zooplankton at time $t$, respectively. $\tau(t)$ represents incubation period of disease and it is a continuously differentiable function. All parameters appeared in system (2.1) are positive. The biological significance of the remaining parameters of system (2.1) are given in Table 1. In addition, system (2.1) satisfies the nonnegative conditions $S(m) = \varphi_1(m) > 0, I(m) = \varphi_2(m) > 0, P(m) = \varphi_3(m) > 0, m\in[-\tau, 0]$, $0\leq\tau(t)\leq\tau, \varphi_i(m)\in\mathcal{C}\big([-\tau, 0), R_{+0}^{3}\big), i = 1, 2, 3$, where $R_{+0}^{3} = \{(S(t), I(t), P(t):S(t) > 0, I(t) > 0, P(t) > 0)\}$.

3.   Dynamics of non-delayed system

When $\tau(t) = 0$, the non-delayed system is shown as follows.

In this section, we will be devoted to discussing some biological properties of non-delayed system (3.1), such as positivity and boundedness of solutions. We also will analyze the local stability and global stability of system (3.1).

3.1.   Positivity and boundedness of solutions

Before analyzing the boundedness of solutions of system (3.1), we discuss the positivity of solutions.

Lemma 3.1. The solution $(S(t), I(t), P(t))$ of system (3.1) with positive initial conditions $(S(0), I(0), P(0))\in \mathbb{R}_{+}^{3}$ remains positive throughout the region for all $t > 0$.

Proof. From system (3.1), we get

Hence, all the solutions of system (3.1) are positive for all $t > 0$. □

Next, we will investigate that all positive solutions are uniformly bounded.

Lemma 3.2. Assume that condition $\mu_{2}\mu_{3}\geq\mu_{1}\mu_{4}$, then the set $\Xi = \big\{(S, I, P)\in \mathbb{R}_{+}^{3}: S+I+ \frac{\mu_1}{\mu_3}P\leq\frac{K}{4r\delta^{*}}(r+\delta^{*})^2\big\}$ is a region of attraction for all solutions initiating in the positive quadrant, where $\delta^{*} = min\{\rho_1, \rho_2\}$.

Proof. To examine the boundedness of system (3.1), we define a function

Differentiating with respect to $t$, we have

From system (3.1), we get

Introducing the positive number $\delta^*$ and rewriting, we obtain

which implies

Choosing $\delta^{*} = min\{\rho_1, \rho_2\}$ and $\mu_{2}\mu_{3}\geq\mu_{1}\mu_{4}$, we have

Using the method of variation of constant for Eq (3.2), we get

where $W_{0} = W(S(0), I(0), P(0))$.

Therefore, it follows that

Then, all the solutions initiating in $\mathbb{R}_{+}^{3}$ of system (3.1) are defined in the region

Thus, all solutions of system (3.1) are uniformly bounded for all $t\geq0$ and $\Xi$ is a positive invariant set. □

3.2.   Existence of equilibriums

The eco-epidemiological system (3.1) has possible nonnegative feasible equilibria as follows: the trivial equilibrium $E_{0}(0, 0, 0)$ always exists. The boundary equilibrium $E_{1}(\tilde{S}, 0, 0)$ exists, where $\tilde{S} = K$.

Assume that $E_{2}(\bar{S}, \bar{I}, 0)$ is the positive equilibrium of system (3.1), where $\bar{S}$ and $\bar{I}$ satisfy the following equations:

According to the second equation of (3.3), we can obtain

Substituting (3.4) into the first equation of (3.3), we get

where

If $(H_1)$: $\beta < r$ holds, then $b_4 < 0$. By Descartes' rule of signs, Eq (3.5) has only one positive root $\bar{S}$, if and only if, one of the following conditions is met:

We assume the assumption $(H_2)$: one of conditions (1)–(4) is true. Therefore, if $(H_1)$ and $(H_2)$ are satisfied, then system (3.1) has only one predator-free equilibrium $E_2(\bar{S}, \bar{I}, 0)$, which is determined by Eqs (3.4) and (3.5).

From the first equation of (3.1) under the case $I = 0$, we have

Substituting (3.6) into the third equation of (3.1), we obtain

If the assumption $(H_3): K > S'$, $\mu_3 > \theta\alpha+\rho_2$ holds, then the equilibrium $E_3(S', 0, P')$ exists.

The positive equilibrium $E^{*}(S^*, I^*, P^*)$ exists if $S^*$, $I^*$, and $P^*$ are the positive solutions of the following equations:

By the second equation of (3.7), we get

The value of $P^*$ is positive if the assumption $(H_4): \rho_1+\frac{\delta}{\sigma+I^*} < \beta S^*$ holds.

Substituting (3.8) into the first and third equations of (3.7), then $S^*$ and $I^*$ are positive solution of following isoclines, such that

We denote the two functions of Eq (3.9) as $F_1(S^*, I^*)$ and $F_2(S^*, I^*)$, respectively. Since $F_1$ and $F_2$ consist of polynomials and rational functions with nonzero denominators, then they are continuous. By selecting $S^*\in[S_{min}, S_{max}]$ and $I^*\in[I_{min}, I_{max}]$, $F_1(S^*, I^*)$ and $F_2(S^*, I^*)$ have opposite signs at the endpoints of the interval, satisfying the conditions of the intermediate value theorem (see Figure 1(a), (b)).

Next, we theoretically prove that $F_1$ and $F_2$ have opposite signs at the endpoints of the interval. Taking the derivative with respect to $S^*$ for $F_1$, we get

Obviously, $\frac{\partial F_1}{\partial S^*} < 0$ when the assumption $(A_1)$ is true, where $(A_1): S^*$ is asymptotically close to K and $\beta S^*(2\alpha + S^*) > \alpha\left(\rho_1 + \frac{\delta}{\sigma + I^*}\right)$.

Taking the derivative with respect to $I^*$ for $F_1$, we have

If the assumption $(A_2)$ holds, where $(A_2): \rho_1 < \beta S^*$ and $\sigma < \alpha$, then $\frac{\partial F_1}{\partial I^*} < 0$. In conclusion, the above analysis shows that $F_1$ is a decreasing function of $S^*$ and $I^*$ under the assumptions $(A_1)$ and $(A_2)$, respectively. Furthermore, we plot the monotonicity curves of $F_1$ with respect to $S^*$ and $I^*$ at fixed values of $I^*$ and $S^*$ (see Figure 2(a), (b)), respectively. Consequently, $F_1$ exhibits opposite signs at the endpoints of the interval.

Next, taking the derivative with respect to $S^*$ for $F_2$, then we get

If the assumption $(A_3)$ is true, where $(A_3): \beta S^* < \frac{\mu_3 \alpha (\rho_1 + \frac{\delta}{\sigma + I^*} - \beta S^*)}{S^* + \alpha}$, $\rho_1 + \frac{\delta}{\sigma + I^*} < \beta (2S^* + I^*)$ and $\mu_4 I^* < \rho_2(I^* + \alpha)$, then $\frac{\partial F_2}{\partial S^*} < 0$. Without loss of generality, we treat $S^*$ as a constant, whereupon $F_2$ becomes a function of $I^*$. Substituting $f^*$ into $F_2$ yields

where $V_3 = \frac{\mu_3 \delta S^*}{\sigma}$, $V_2 = (\theta-\beta) S^*$, $A_1 = \beta \alpha S^*+\mu_4 \rho_1$, $V_0 = -\rho_1 \alpha \cdot \frac{\mu_3 S^*}{S^* + \alpha}+ \rho_1 \rho_2 \alpha + \rho_1 \alpha \theta S^*+ \beta \alpha \cdot \frac{\mu_3 S^{*2}}{S^* + \alpha}- \beta \alpha \rho_2 S^* - \beta \alpha \theta S^{*2}$. It is direct from Eq (3.10) that

It is obvious that $V_3 > 0$ and $V_1 > 0$. If $V_2\geq0$, then $F_2'(I^*) > 0$ for $I^*\in(0, \infty)$. It means that $F_2(I^*)$ is an increasing function of $I^*$ in $\in(0, \infty)$. If $V_2 < 0$ and $V_2^{2} > 3V_{3}V_{1}$, let $F_2'(I^*) = 0$, and we derive

When $I^*\in(0, \frac{-2\sqrt{V_{2}^2-3V_{3}V_{1}}-V_2}{3V_{3}})$, $F_2'(I^*) > 0$. However, when $I^*\in(\frac{-2\sqrt{V_{2}^2-3V_{3}V_{1}}-V_2}{3V_{3}}, \frac{2\sqrt{V_{2}^2-3V_{3}V_{1}}-V_2}{3V_{3}})$$F_2'(I^*) < 0$. That is, $F_2(I^*)$ is an increasing function as $I^* \in (0, \frac{-2\sqrt{V_{2}^2-3V_{3}V_{1}}-V_2}{3V_{3}})$ and a decreasing function as $I^* \in (\frac{-2\sqrt{V_{2}^2-3V_{3}V_{1}}-V_2}{3V_{3}}, \frac{2\sqrt{V_{2}^2-3V_{3}V_{1}}-V_2}{3V_{3}})$. Similarly, we plot the monotonicity curves of $F_2$ with respect to $S^*$ and $I^*$ at fixed values of $I^*$ and $S^*$ (see Figure 2(c), (d)), respectively. Therefore, $F_2$ exhibits opposite signs at the endpoints of the interval.

Additionally, the location of their intersection point is showed in Figure 1(c). From Figure 1(c), the figures of the two functions $F_1$ and $F_2$ intersect as $S^*$ and $I^*$ vary. As shown in the phase plane projection (see Figure 3), the intersection points of $F_1(S^*, I^*)$ and $F_2(S^*, I^*)$ are confined to the positive domain, confirming that the positive equilibrium $E^{*}(S^*, I^*, P^*)$ exists.

Furthermore, in Figure 4, we plot the isoclines of Eq (3.9) by choosing the parametric values mentioned in Table 1. It is clear from Figure 4 that the isoclines of Eq (3.9) intersect in the interior of the positive quadrant. Therefore, system (3.1) has at least a unique positive equilibrium.

3.3.   Stability analysis

In this subsection, the local stability of all equilibria of system (3.1) and the global stability of the positive equilibrium will be stated in the following theorems.

Theorem 3.1. The trivial equilibrium $E_0$ is always unstable.

Proof. The Jacobian of system (3.1) is given by

where

The characteristic equation of system (3.1) at $E_0$ is

Since Eq (3.12) has a unique positive root $\lambda = r$, the trivial equilibrium $E_0$ is always unstable. □

Theorem 3.2. If the assumption $(H_5): \beta K < \rho_1+\frac{\delta}{\sigma}$, and $\frac{\mu_3 K}{K+\alpha} < \rho_2 + \theta K$ is satisfied, then $E_1$ is stable, but it is unstable when one case of the reverse inequality holds.

Proof. From (3.11), the characteristic equation of system (3.1) at $E_1$ is

Clearly, the characteristic values around $E_1$ are $\lambda_1 = -r < 0$, $\lambda_2 = \beta K -\rho_1-\frac{\delta}{\sigma}$, $\lambda_3 = \frac{\mu_3 K}{S+\alpha}-\rho_2-\theta K$. If the assumption $(H_5)$ is true, then all characteristic values of Eq (3.13) have negative real parts, thus $E_1$ is locally asymptotically stable. Conversely, if $\beta K > \rho_1+\frac{\delta}{\sigma}$ or $\frac{\mu_3 K}{K+\alpha} > \rho_2 + \theta K$, then Eq (3.13) has at least one characteristic value with the positive real part and $E_1$ is unstable. □

We give some assumptions as follows:

$(H_6):\frac{\mu_3\bar{S}}{\bar{S}+\alpha}+\frac{\mu_4\bar{I}}{\bar{I}+\alpha} < \rho_2+\theta(\bar{S}+\bar{I})$,

$(H_7):r(1-\frac{2\bar{S}+\bar{I}}{K})+\beta(\bar{S}+\bar{I}) < \rho_1+\frac{\sigma\delta}{(\bar{I}+\sigma)^2}$,

$(H_8):r\big(1-\frac{2\bar{S}+\bar{I}}{K}\big)\big(\beta\bar{S}-\rho_1-\frac{\sigma\delta}{(\bar{I}+\sigma)^2}\big) +2\beta^2\bar{S}\bar{I}+\frac{r\beta\bar{S}\bar{I}}{K} > \beta\bar{I}\big(\rho_1+\frac{2\sigma\delta}{(\bar{I}+\sigma)^2}\big)$.

Theorem 3.3. If the assumptions $(H_1)$, $(H_2)$, $(H_6)$–$(H_8)$ are satisfied, then the equilibrium $E_2(\bar{S}, \bar{I}, 0)$ is locally asymptotically stable.

Proof. The characteristics equation of system (3.1) at the equilibrium $E_2$ is given by

One eigenvalue of Eq (3.14) is $\frac{\mu_3\bar{S}}{\bar{S}+\alpha}+\frac{\mu_4\bar{I}}{\bar{I}+\alpha}-\rho_2-\theta(\bar{S}+\bar{I})$, which is negative when the assumption $(H_6)$ holds.

In Eq (3.14), the quadratic term gives two additional eigenvalues which satisfy the following conditions:

Now, $\lambda_1+\lambda_2 < 0$ and $\lambda_1\cdot\lambda_2 > 0$ if the assumptions $(H_7)$ and $(H_8)$ hold.

Therefore, the equilibrium $E_2$ is locally asymptotically stable if the conditions $(H_6)$–$(H_8)$ hold. □

Theorem 3.4. If the assumptions $(H_3)$, $(H_9):\beta S' < \frac{\mu_2 P'}{\alpha}+\rho_1+\frac{\delta}{\sigma}$, and $(H_{10}):r < \frac{2rS'}{K}+\frac{\mu_1 \alpha P'}{(S'+\alpha)^2}, \frac{\mu_3 \alpha P'}{(S'+\alpha)^2} > \theta P'$ are true, then the equilibrium $E_3(S', 0, P')$ is locally asymptotically stable.

Proof. The Jacobian matrix corresponding to the disease-free equilibrium $E_3$ is:

where

One eigenvalue is $\beta S'-\frac{\mu_2 P'}{\alpha}-\rho_1-\frac{\delta}{\sigma} > 0$ if the assumption $(H_9)$ holds. The others are the eigenvalue of sub-matrix

The eigenvalues of sub-matrix $J'(E_3)$ are negative if $tr\left(J'(E_3)\right) < 0$ and $det\left(J'(E_3)\right) > 0$, where

In fact, $tr\left(J'(E_3)\right) < 0$ and $det\left(J'(E_3)\right) > 0$ if the condition $(H_{10})$ holds. □

Theorem 3.5. If the assumptions $(H_4)$ and $(H_{11}): V_1 > 0, V_2 > 0, V_1 V_2 -V_3 > 0$ hold, then the positive equilibrium $E^*(S^*, I^*, P^*)$ is locally asymptotically stable.

Proof. The Jacobian matrix corresponding to the equilibrium $E^{*}$ is given by

where $e_{ij} = j_{ij}$ at $E^{*}$. Here, $e_{13} < 0$, $e_{23} < 0$ and $e_{33} = 0$. The characteristic equation of $J(E^{*})$ is given by

where

From Routh-Hurwitz criteria, if the assumption $(H_{11})$ holds, then system (3.1) is locally asymptotically stable at the positive equilibrium $E^{*}$. □

Next, we give the following lemma before proving the global stability of the positive equilibrium $E^{*}$ of system (3.1).

Lemma 3.3. If $f$ is defined on $[0, \infty)$ and nonnegative such that $f$ is integrable and uniformly continuous on $[0, \infty)$, then $\lim\limits_{t\to\infty}f(t) = 0$.

Theorem 3.6. If the assumptions $(H_4)$ and $(H_{12}): o_{11} > 0, o_{11}o_{22}-o_{12}^{2} > 0$, $o_{11}o_{22}o_{33}+o_{12}o_{23}o_{31}+o_{13}o_{21}o_{32}-o_{13}^2o_{22}-o_{23}^2o_{11}-o_{12}^2o_{33} > 0$ hold, then the positive equilibrium $E^*$ of system (3.1) is globally asymptotically stable.

Proof. Define the Lyapunov function for system (3.1) as

where $k$ is an applicable positive constant.

Differentiating Eq (3.16) with respect to time $t$, we have

where $\psi_1 = (S+\alpha)(S^*+\alpha)$, $\psi_2 = (I+\sigma)(I^*+\sigma)$, $\psi_3 = (I+\alpha)(I^*+\alpha$).

The above expression can be reduced as

where $Y^{T} = (|S(t)-S^{*}|, |I(t)-I^{*}|, |P(t)-P^{*}|)$, and $O = (o_{ij})_{3\times3}, i, j = 1, 2, 3$. All entries of the matrix $O_{3\times3}$ are defined as

The matrix $O$ is positive if, and only if, all the principal minors of $O$ are positive. If $(H_{12})$ is held, then all the principal minors of $O$ are positive.

Therefore, from (3.17), we have that

where $\Lambda$ is the smallest eigenvalue. From (3.18), we have

By definition of $V$ and (3.19), $S(t)$, $I(t)$, and $P(t)$ are bounded uniformly on $[0, +\infty)$, and $\frac{\mathrm{d}S}{\mathrm{d}t}$, $\frac{\mathrm{d}I}{\mathrm{d}t}$, and $\frac{\rm{d}P}{\mathrm{d}t}$ are also bounded on $[0, +\infty)$. By Lemmas (3.3) and (3.19), we conclude that

Thus, system (3.1) is globally asymptotically stable around the positive equilibrium $E^{*}$. □

4.   Dynamics of delayed system and oscillation suppression

In this section, we will investigate the stability and the existence of Hopf bifurcation of the delayed system.

4.1.   Local stability and Hopf bifurcation

When the incubation period of the virus is constant, system (2.1) becomes as follows:

Suppose $\tilde{S}(t) = S(t)-S^*, \tilde{I}(t) = I(t)-I^*, \tilde{P}(t) = P(t)-P^*$. For simplicity, $\tilde{S}$, $\tilde{I}$, and $\tilde{P}$ still are written as $S$, $I$, and $P$, respectively. Then, the linearized system at the positive equilibrium $E^*(S^*, I^*, P^*)$ is given by

where

Let

The characteristic equation of system (4.2) at $E^*$ is

where

Suppose that Eq (4.3) admits a root $\lambda = i\omega(\omega > 0)$. Then, substituting it in Eq (4.3) and separating real and imaginary parts, we have

Squaring two equations of (4.4), respectively, and then summing them, we obtain

where $Q_1 = P_1^2-2 P_2-P_4^2, \; Q_2 = P_2^2-2 P_1 P_3+2P_4 P_6-P_5^2, \; Q_3 = P_3^2-P_6^2$.

Choosing $\omega^2 = z$, Eq (4.5) can be rewritten as

Differentiating both sides of the equation with respect to $z$ yields:

where the roots of $f^{\prime}(z)$ can be expressed as,

According to the reference ^[39] and the formula Fan in ^[40], we can get the discriminant of Eq (4.6)

where $\tilde{A} = Q_{1}^2-3Q_2, \tilde{B} = Q_1 Q_2 -9Q_3, \tilde{C} = Q_{2}^2-3Q_1 Q_3.$

Lemma 4.1. From Eq (4.6), the following conclusions are true.

$(i)$ If $(H_{13}): Q_3\geq 0, Q_{1}^2-3Q_{2}\leq 0$ holds, then $f(z)$ is monotonically increasing for $z\in[0, +\infty)$. Therefore, $f(z)$ has no positive roots on $[0, +\infty)$, i.e., Eq (4.6) has no positive equilibria.

$(ii)$ If $\Delta > 0$, then Eq (4.6) only has one positive equilibrium.

$(iii)$ If $\Delta = 0$, then Eq (4.6) has two positive equilibria.

$(iv)$ If $\Delta < 0$, then Eq (4.6) has three positive equilibria.

Without loss of generality, suppose that Eq (4.6) has three positive roots, denoted as $z_1$–$z_3$, then Eq (4.5) has three roots: $\omega_1 = \sqrt{z_1}, \omega_2 = \sqrt{z_2}, \omega_3 = \sqrt{z_3}$. According to Eq (4.4), we have

Therefore, when $\tau = \tau_{k}^{j}$, then $\pm i \omega_{k}$ is a pair of pure imaginary roots of Eq (4.3), and the other roots have nonnegative real parts.

Next, we will find the transversality condition for the occurrence of Hopf bifurcation. Let $\lambda(\tau) = \alpha(\tau)+i\beta(\tau)$ be the root of Eq (4.3) when $\tau = \tau_{k}^{j}$, satisfying $\alpha(\tau_{k}^{j}) = 0$ and $\beta(\tau_{k}^{j}) = \omega_k$, $k = 1, 2, 3, j = 0, 1, 2, \cdots$. Denoting $\tau_{0} = \mathop{min}\limits_{k = 1, 2, 3, j = 0, 1, \cdots}\{\tau_{k}^{j}\}$, we have the transversality condition.

Theorem 4.1. If $z_k = \omega_{k}^2$ and $f'(z_k)\neq0$ hold, then $\frac{\mathrm{d}Re(\lambda(\tau))}{\mathrm{d}\tau}\mid_{\lambda = i\omega_0}\neq0$.

Proof. Differentiating Eq (4.3) with respect to $\tau$, we have

Substituting $\lambda = i\omega_{0}$ into Eq (4.7) and getting the real part, we obtain

So, we have

which is the required transversality condition. □

Based on the above analysis, the following theorem on stability of $E^{*}$ and Hopf bifurcation can be obtained.

Theorem 4.2. Consider the system (4.1) with $\tau > 0$.

(i) If assumption $(H_{13})$ holds, then the positive equilibrium $E^{*}$ is locally asymptotically stable for all $\tau > 0$.

(ii) If $\Delta\geq0$ or $\Delta < 0$ holds, then the positive equilibrium $E^{*}$ is locally asymptotically stable when $\tau\in[0, \tau_{0})$, but is unstable when $\tau > \tau_{0}$.

(iii) If all the conditions in $(ii)$ hold and $f'(z_k)\neq 0$, then the system (4.1) undergoes a Hopf bifurcation at the positive equilibrium $E^{*}$ when $\tau = \tau_{0}$.

4.2.   Hopf bifurcation analysis via MMS

The utilization of MMS is frequently employed for the examination of phenomena across varying temporal and spatial dimensions. In this subsection, we will apply this approach to system (2.1). The method is based on references ^{[41,42,43,44,45,46]}.

According to the previous analysis, the system has only one positive equilibrium $X^{*} = E^{*}(S^{*}, I^{*}, P^{*})$. To begin, it is imperative to linearize the system (4.1) at the equilibrium $E^*$, defined as $\hat{S}(t) = S(t)-S^{*}$, $\hat{I}(t) = I(t)-I^{*}$, $\hat{P}(t) = P(t)-P^{*}$. For simplicity, $\hat{S}$, $\hat{I}$, and $\hat{P}$ still are written as $S$, $I$, and $P$, respectively. So we obtain a differential equation about the variable $X$:

where $X = (S, I, P)^{T}$, $X_{\tau} = X(t-\tau)$.

For MMS, we first define the new timescales and the time derivatives, respectively, as:

where $D_{k} = \frac{\partial}{\partial T_{k}}, k\in N\cup\{0\}$. To investigate Hopf bifurcation and use MMS, we need to perturb the bifurcation parameter in the system (4.8) near its critical value, that is, $\tau = \tau_{0}+\epsilon\tau_{\epsilon}$, where $\epsilon$ is a small positive parameter. According to reference ^[41], we assume that the solution for system (4.8) is of the following form:

Furthermore, adopting the different timescales, the delayed solutions of system (4.8) can be asymptotically expanded up and expressed in the following form:

Substituting Eqs (4.10) and (4.11) into system (4.8), using multivariate Taylor expansion, and collecting the terms by different powers of $\epsilon$, we first consider the lowest order of $\epsilon$,

where $F_X = \frac{\partial F}{\partial X}$, $F_{X_\tau} = \frac{\partial F}{\partial X_{\tau}}$. Further, we suppose that the solution of Eq (4.12) can be denoted as:

where

For simplicity, we will write the above equations in the following form:

Substituting Eq (4.13) into Eq (4.12), we get

Then, we derive that $M_{12}$, $M_{13}$, $N_{12}$, and $N_{13}$ can be represented by $M_{11}$ and $N_{11}$, respectively. Furthermore, we get the following system of linear equations:

Similar to Eq (4.9), denote $\frac{\mathrm{d}}{\mathrm{d}T_k}$ by $D_{k}$ for $k = 0, 1, \cdots.$ From Eq (4.14), we have

Substituting Eqs (4.10) and (4.11) into system (4.8), and using Taylor expansion, we get the following equation at the second order of $\epsilon$:

where

In order to ensure that the equation is solvable, the asymptotic oscillation terms including $\sin(\omega T_{0})$ and $\cos(\omega T_{0})$ must be eliminated, that is, suppose that the coefficients of $\sin(\omega T_{0})$ and $\cos(\omega T_{0})$ in Eq (4.15) are zero. From Eq (4.15), we can obtain the following equations:

After eliminating the asymptotic oscillation terms, the Eq (4.15) can be rewritten in the following form:

In order to solve Eq (4.16) and satisfy the solvability conditions, the solution of Eq (4.16) is as follows:

where $X_2 = X_2(T_0, T_1, T_2, \cdots)$, $M_2 = M_2(T_1, T_2, T_3, \cdots)$, $N_2 = N_2(T_1, T_2, T_3, \cdots)$, $O_2 = O_2(T_1, T_2$, $T_3, \cdots)$, and $O_2$ can be denoted by $M_1$ and $N_1$, namely,

Therefore, $O_2$ can be denoted by $M_{11}$ and $N_{11}$.

Next, similar to $\epsilon$ and $\epsilon^2$, we repeat the above procedures and collect the terms of powers of $\epsilon^3$ from Eq (4.8),

where $i = 1, 2$,

where $D_{00} = \frac{\partial^2}{\partial T_{0}^2}$, $D_{01} = \frac{\partial^2}{\partial T_{0} \partial T_{1}}$, and $D_{11} = \frac{\partial^2}{\partial T_{1}^2}$. In general, $g_2$ contains the Taylor expansion of the delay term $X_\tau$, the interaction of the system states $X_1$ and $X_2$, and the higher order time derivative term.

Similar to previous power, we require to eliminate the asymptotic oscillation terms to satisfy the solvability. After calculation, we obtain that the expressions of $D_{2}M_{11}$ and $D_{2}N_{11}$ are related to $M_{11}$ and $N_{11}$ at $\epsilon^3$, which are too complex to be showed here. Ultimately, restoring to the original timescale, we can obtain

So as to conduct the normative analysis, the following polar coordinate transformations on Eq (4.17) are:

As a consequence, when $\tau = \tau_0+\epsilon\tau_\epsilon$ the amplitude frequency equation of system (4.17) can be written as

where $p_i(\varepsilon, \tau_\varepsilon), i = 0, 1, 2, 3$, depends on the parameters $\varepsilon$ and $\tau_\varepsilon$. They describe the dynamic behavior of the system in polar coordinates, including the changes in amplitude $R(t)$ and phase $\varphi(t)$. By analyzing these coefficients, the stability and periodic behavior of the system can be understood.

4.3.   Oscillation suppression based on periodic delay

Bifurcation is a fundamental concept in dynamical systems. It can cause system instability, lead to periodic oscillations, or even chaotic behavior ^[24,47,48]. To enhance system performance, reliability, stability, and to address potential challenges in practical applications, many researchers have focused on bifurcation control. Generally speaking, the adopted methods of controlling bifurcations include: parameter tuning ^[47,49], feedback control ^[50,51], delayed feedback control ^[43,44,52], hybrid control ^[48,53], and dual control ^[54,55,56].

As we know, the delay in system (4.1) can lead to periodic oscillations and Hopf bifurcation. Furthermore, we set $\tau = \tau_0+\epsilon\tau_\epsilon$, where $\epsilon\tau_\epsilon > 0$. In this subsection, we attempt to introduce a perturbation control. That is, applying time-varying perturbation to the delay achieves the goal of oscillation suppression ^[43,57]. Specifically, let

where $L$ and $\Omega$ represent magnitude and frequency of the perturbation, respectively, and $\tau_m = \tau_0 + \epsilon \tau_\epsilon$. Assume that $\Omega$ and $L$ are small, so that $L \sin(\Omega t)$ can be considered as a perturbation to $\tau_m$. Consequently, the method of multiple timescales can be applied to analyze the time-varying delay system. Substituting Eq (4.19) into system (4.8), we have

To derive the amplitude-frequency equation for system (4.17) when $\tau(t) = \tau_0 + \epsilon \tau_\epsilon + L \sin(\Omega t)$, a re-scaling is necessary: $\epsilon R(t) \to R(t)$. Following the procedure in Subsection 4.2, we obtain the amplitude-frequency equation as follows:

where $p_1(t) = q_0+q_{11}\sin(\Omega t)+q_{12}\cos(\Omega t)$, $p_0(t) = r_0+r_{11}\sin(\Omega t)+r_{12}\cos(\Omega t)$, $q_0 = q_0(\epsilon\tau_\epsilon, L, \Omega)$, $q_{11} = q_{11}(\epsilon\tau_\epsilon, L, \Omega)$, $q_{12} = q_{12}(\epsilon\tau_\epsilon, L, \Omega)$, $r_0 = r_0(\epsilon\tau_\epsilon, L, \Omega)$, $r_{11} = r_{11}(\epsilon\tau_\epsilon, L, \Omega)$, $r_{12} = r_{12}(\epsilon\tau_\epsilon, L, \Omega)$. $p_2(t)$ and $p_3(t)$ are functions of $\epsilon$.

Because the first equation of (4.21) is only relative to $R(t)$, the second equation depends on $R(t)$. Therefore, we only pay attention to the first equation of (4.21):

Note that Eq (4.22) is a Bernoulli equation with variable coefficients, and its solution can be derived by using the following lemma.

Lemma 4.2. The solution to Eq (4.22) can be expressed analytically as follows:

where $C$ is the initial value.

It is worth noting that the decay of the solution of Eq (4.22) corresponds to the decay of a Bursting solution of Eq (4.8). In this subsection, a sufficient condition is derived for the solution to Eq (4.22) to decay. So, we have the following theorem.

Theorem 4.3. A sufficient condition for the decay of solutions of system (4.21) is $q_0 < 0$, where $q_0 = q_0(\epsilon \tau_\epsilon, L, \Omega)$.

Proof. $C > 0$ in Eq (4.22). Therefore, Eq (4.22) can be written as

Utilizing the Fourier series expansion and combining with Eq (4.23), we get

where $\bar{C} = C^{-2}$, $Y_1(t)$, and $Y_2(t)$ are periodic functions with a period of $\frac{2 \pi}{\Omega}$. Clearly, $\bar{C} > 0$, and the long-term behavior of $\tilde{v}(t) = \bar{C}e^{-2q_0(t)}Y_1(t) + Y_2(t)$ is governed by the sign of $q_0(t)$. When $q_0(t)\geq0$, $\bar{C}e^{-2q_0(t)}Y_1(t)$ decays exponentially, thus the long-time behavior of $\tilde{v}(t)$ is determined by $Y_2(t)$, leading to periodic oscillation in $R(t)$. On the other hand, when $q_0(t) < 0$, then $\tilde{v}(t)$ grows exponentially, which causes the oscillation of $R(t) = \frac{1}{\sqrt{\tilde{v}(t)}}$ to die out exponentially over time. In summary, if the perturbation parameters of the delay satisfy certain conditions that ensure $q_0(t) < 0$, then the oscillation of system (2.1) can be controlled. Consequently, a critical boundary can be determined by solving $q_0(\epsilon\tau_\epsilon, L, \Omega) = 0$. □

4.4.   Global asymptotic stability of the positive equilibrium

Theorem 4.4. If $(H_4)$, $(H_{14})$, and $(H_{15})$ hold, then the positive equilibrium $E^*$ of system (2.1) is globally asymptotically stable.

Proof. Suppose a small disturbance occurs near the positive equilibrium $E^*$ of the system (2.1); let $w_1(t) = S(t)-S^*$, $w_2(t) = I(t)-I^*$, $w_3(t) = P(t)-P^*$.

By means of neglecting the higher-order nonlinear terms and applying a linear approximation, the following linearized system is obtained:

According to Lyapunov stability theory, the stability of the positive equilibrium of system (2.1) is equivalent to the stability of the zero solution of system (4.26).

Define a Lyapunov-Krasouskii function near the positive equilibrium $E^*$ as

where $r_1, r_2$ are appropriate positive constants.

Taking the derivative of the above $V$ function with respect to $t$, we have

where $\lvert\dot{\tau}(t)\rvert\leq\delta$, and $\delta$ is a positive constant. Substituting (4.26) into (4.28), we can obtain the following formulas

Using the Young inequality to deal with the cross terms in the above expression and perform scaling, we can obtain

where $\varepsilon > 0$, $W^T = (\lvert S-S^* \rvert, \lvert I-I^* \rvert, \lvert P-P^* \rvert)$, and $Q_{3\times3}$ is a symmetric matrix can be defined as:

where

The matrix $Q$ is positive if, and only if, all the principal minors of $Q$ are positive. If $(H_{14}): q_{11}, q_{11}q_{22}, q_{11}q_{22}q_{33}-q_{13}^2q_{22}-q_{23}^2q_{11}$ are positive, then all the principal minors of $Q$ are positive. Further, if $(H_{15}): \frac{r_1}{2}(1-\dot{\tau}_{max}) > \frac{\beta^2}{2\varepsilon}, \frac{r_2}{2}(1-\dot{\tau}_{max}) > 0$ is true, then system (2.1) is globally asymptotically stable near the positive equilibrium $E^*$. □

5.   Numerical simulations

In this section, we have taken the same set of parameter values given in Table 1 to verify the theoretical results, and give some biological explanations. To start, we perform a sensitivity analysis on some parameters of system (3.1).

To identify the parameters that significantly influence the output variables of system (3.1), we conduct the global sensitivity analysis on selected parameters. From system (3.1), we calculate the partial rank correlation coefficients (PRCCs) among parameters $r$, $K$, $\beta$, $\delta$, $\sigma$, and $\theta$. Using the Latin hypercube sampling (LHS) method, we performed $5000$ simulations of parameters in system (3.1). The baseline values of these parameters are presented in Table 2.

Based on the parameter values in Table 2, we analyze the impact of selected parameters on the correlation of the infected phytoplankton. Then, by generating scatter plots and conducting $5000$ parameter samplings, we obtained the sampling results (see Figure 5) and dot plots (see Figure 6). The trend of monotonic increase (decrease) indicates the positive (negative) correlation between parameters and the model outputs. From Figure 6, these parameters show periodic correlations with system outputs, where $r$, $K$, and $\theta$ exhibit positive correlations with model outputs, $\beta$ demonstrates negative correlations with model outputs, and $\delta$ and $\sigma$ show no significant correlation with model outputs.

Using the parameter values in Table 1, we obtain the unique predator-free equilibrium $E_2(56.6756, 92.7445, 0)$ and the positive equilibrium $E^*(118.9001, 44.8221, 13.6165)$. For non-delayed system (3.1), the predator-free equilibrium $E_2$ and the positive equilibrium $E^*$ are locally asymptotically stable according to Theorems 3.3 and 3.5, respectively (see Figures 7 and 8). Thus, the susceptible phytoplankton, the infected phytoplankton, and zooplankton population are stable.

Based on Theorem 3.6, $E^*$ is globally asymptotically stable (see Figure 9). This implies that regardless of variations in initial conditions, both the phytoplankton and zooplankton populations will converge to a stable equilibrium state.

For the time-varying delay system (2.1), the positive equilibrium of such a system remains unchanged by selecting the parameter values in Table 1. First, when we choose incubation period of disease $\tau(t) = \frac{\lvert \sin(t^3) \rvert}{1+t^2}$, we find that system (2.1) is globally asymptotically stable at the positive equilibrium (see Figure 10), which shows that the conclusion of Theorem 4.1 is correct. However, the solution of system (2.1) is unstable when $\tau(t) = \lvert\sin(t^2)\rvert+\frac{5}{2}$ (see Figure 11). At the same time, when $\tau(t) = \lvert\sin(t^2)\rvert+\frac{5}{2}$, the susceptible and infected phytoplankton population show periodicity (see Figure 11(a), (b)), and the zooplankton become extinct. Through numerical simulation, we find that time-varying delay is a key factor affecting the stability of system (2.1). When the value of $\tau(t)$ is small enough, system (2.1) can maintain stability, which provides the possibility for taking some measures to control the widespread outbreak of the disease. On the contrary, when the time-varying delay is relatively large, the system may become unstable, and even chaos may occur. Therefore, early prevention and diagnosis are crucial for effectively preventing the spread of the disease and reducing it.

For delayed system (4.1), the system has the same positive equilibrium $E^*$. We can calculate that $\tau_0 = 0.1093$. Furthermore, we find that when $\tau = \tau_0 = 0.1093$, system (4.1) undergoes a Hopf bifurcation at $E^*(118.9001, 44.8221, 13.6165)$, i.e., periodic solutions bifurcate from $E^*$, as shown in Figure 12. According to the Theorem 4.2, we find that $E^*$ is locally asymptotically stable when $\tau = 0.05 < \tau_0$ (see Figure 13), but $E^*$ is unstable when $\tau = 0.12 > \tau_0$ (see Figure 14). That is, when the incubation period of disease transmission is less than the critical value, phytoplankton and zooplankton can coexist for a long time. When the delay exceeds the threshold, susceptible phytoplankton, infected phytoplankton, and zooplankton can still coexist at positive densities, but their densities oscillate periodically over time.

Utilizing the method of multiple timescales discussed in Subsection 4.2, we can derive the following amplitude-frequency equation:

Multiplying both sides of the first equation of Eq (5.1) by $\epsilon$ simultaneously yields:

Let $\epsilon\frac{\mathrm{d}R}{\mathrm{d}t} = 0$, and we can get a nontrivial steady-state solution for $\epsilon R(t)$, that is,

By substituting (5.3) into the second equation of (5.1) and integrating it, we obtain

where $\phi_0$ is initial phase.

Next, we will primarily focus on the case of time-varying delay by utilizing periodic time-delay perturbation to suppress the oscillations of system. Furthermore, by employing the method of multiple scales to investigate the time-varying delay system, the following amplitude-frequency equation can be derived:

where

and

Assume that $\Omega = 0.02$ and $\epsilon\tau_{\epsilon} = 0.0001$. By solving $q_0 = 0$, we can obtain the critical value of $L$ that induces attenuation of oscillation, which is $L = 0.9622$. When $L_1 = 0.9986$ exceeds the critical value $L = 0.9622$ (see Figure 13), we observe that the oscillation degree is significantly weakened, as shown in Figure 15. When $L_2 = 0.9484 < L = 0.9622$, the oscillation gradually decays and the system eventually becomes stable (see Figure 16).

The above analysis reveals that when time-varying perturbation is applied to the time delay, the oscillation of the system gradually decreases and even disappears. Thus, the oscillation suppression by periodic delay is effective in this paper. Moreover, when the disease latency delay is an approximate periodic disturbance function, the degree of instability of phytoplankton and zooplankton can gradually decrease, ultimately leading to a stable coexistence.

6.   Discussion and conclusions

To the best of our knowledge, only a few articles have incorporated time delays under different fundamental assumptions and analyzed their impact on system dynamics. For example, time delay in the viral infection process may destabilize phytoplankton density while zooplankton density remains unchanged ^[24], time delay in toxin liberation could destabilize an otherwise stable equilibrium ^[25], and a sufficiently large delay in fear-induced prey reproduction may trigger double Hopf bifurcation ^[30]. However, the influence of climate and environmental factors on time delays has not yet been considered. To address this gap, we propose treating delays as dynamic variables, which significantly improves the accuracy of population behavior characterization. To enhance the realism of the model, we considered the time delay as the incubation period of disease transmission, and replaced the constant delay with the time-varying delay. Thus, we developed a time-varying delay eco-epidemiological model incorporating toxicity, treatment, and Holling Ⅱ functional response in this paper.

Then, we explored the dynamics of systems without delay, with constant delay, and with time-varying delay. First, we demonstrated the positivity and boundedness of the solution, and analyzed the local stability of five equilibriums of the non-delayed system (3.1), respectively. When there is no delay, we found that both the predator-free equilibrium and the positive equilibrium of system (3.1) are locally asymptotically stable (see Figures 7 and 8). Additionally, we investigated the global stability of the positive equilibrium by constructing an appropriate Lyapunov function (see Figure 9). Next, we analyzed the impact of time-varying delay on the stability of the system. When the time-varying delay is sufficiently small, the positive equilibrium is global asymptotically stable (see Figure 10). However, as the time-varying delay increases, the system may exhibit complex dynamic behaviors (see Figure 11). For the model with time delay, we studied the local stability of the positive equilibrium and Hopf bifurcation of system (4.1) by taking time delay as bifurcation parameter. We observed that Hopf bifurcation occurs when $\tau = \tau_0$ (see Figure 12). Our analysis revealed that when $\tau = 0.05$ (less than the critical value 0.1093), the system remains locally asymptotically stable at the positive equilibrium (see Figure 13). However, when $\tau = 0.12$ (greater than 0.1093), the system becomes unstable (see Figure 14). Moreover, numerical simulations were conducted using computational software to validate the theoretical findings. Our results demonstrated that time-varying delays generate rich and complex dynamics.

We found that time delay can cause oscillations in system (4.1) through Hopf bifurcation. For our eco-epidemiological model, sustained oscillations between phytoplankton and zooplankton can undermine the stability of the system. Therefore, it is crucial to understand the factors that contribute to these oscillations and to identify strategies that can suppress the oscillation. Applying time-varying perturbation to the delay could serve as an effective control strategy to suppress oscillation, which is a delay-based control strategy. Using the method of multiple timescales, we derived the quantitative relationship between time delay and the periodic oscillation resulting from the Hopf bifurcation, as well as the critical value of the perturbation amplitude necessary to effectively control the oscillatory behavior of system. As shown in Figure 15, when $L_1 = 0.9986 > L = 0.9622$, the oscillation is significantly reduced compared to those in Figure 13. Furthermore, when $L_2 = 0.9484 < L = 0.9622$, the system eventually stabilizes (see Figure 16). The numerical results demonstrate that periodic perturbation of the time delay can successfully suppress the oscillations in systems. Thus, by combining the method of multiple timescales with periodic delay perturbations, we effectively suppress oscillatory behavior induced by Hopf bifurcation, providing a novel approach for stability control in time-delay systems.

Moreover, several interesting topics could be explored in the future. For example, by considering that the release of toxins is not an instantaneous process and may be influenced by factors such as climate and temperature, we will consider another time-varying delay, denoting the delay in the release of phytoplankton toxin. Additionally, we will take into account the fear response of prey to predator and investigate its effect on prey population growth. Thus, we introduce the fear effect, where $k_1$ represents the degree of fear to reduce the growth of prey, and $k_2$ stands for reduced fear of disease transmission. Because fear of predator reduces the foraging activity of prey populations, which reduces disease transmission between prey, it is meaningful to consider the following time-varying delay eco-epidemiological model incorporating fear effect,

The impact of the revisions on the the dynamics of model remains to be explored, and this will be discussed in the future.

Use of AI tools declaration

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

Acknowledgments

This work is supported by the National Natural Science Foundation of China (Grant Nos. 12161054 and 12161011), Funds for Innovative Fundamental Research Group Project of Gansu Province(Grant No. 24JRRA778), and the Doctoral Foundation of Lanzhou University of Technology.

Conflict of interest

All authors declare no conflicts of interest in this paper.