Global dynamics of a diffusive phytoplankton-zooplankton model with toxic substances effect and delay

1.   Introduction

In aquatic systems such as rivers, lakes and oceans, harmful algal blooms (HABs) have attracted considerable scientific attention in recent years. The effect of toxin-producing phytoplankton (TPP) on zooplankton is one reason for such attention. Researchers have made great progress in mathematical modeling in this field ^{[1,2,3,4,5,6]}. To fully understand the mechanism of planktonic blooms and to formulate reasonable control measures, Chattopadhyay et al.^[3] build a ODE phytoplankton-zooplankton model with toxin effect, using data gathered from the coastal region of West Bengal and part of Orissa, India. These researchers investigated plankton models with different predational response functions and toxin liberation processes to obtain rich dynamics. They also considered the diffusivity of plankton influenced by ocean currents and tides, and the time delay effect of the phytoplankton toxin on zooplankton. Wang et al. ^[7] consider a ODE plankton system with Holling II response function and linear toxin processes. These authors concluded that the system undergoes Hopf bifurcation at positive equilibrium and Hopf-transcritical bifurcation when the parameters satisfy a particular condition. In spatial plankton models, for example, Chaudhuri et al. ^[4] consider that the system includes both non-toxic and toxic phytoplankton and the system shows toxic phytoplankton-induced spatiotemporal patterns when one phytoplankton releases toxin. Further, to describe the reduction of zooplankton due to toxin-producing phytoplankton, plankton systems with discrete delay are presented in References ^[2,8] to show the effect of delay. There is o large body of literatures describing the dynamics of aquatic models ^[6,9,10,11].

In fact, the problem we analyze is based on Chen et al.^[1], which considers a model as follows:

They conclude that the plankton model (1.1) occurs as a bistable phenomenon.

In real-world conditions, plankton affected by tides and turbulence may move and diffuse across lakes and seas. Thus, such diffusion should be considered in studying the dynamics of plankton models. Meanwhile, the effect of phytoplankton toxin on zooplankton is in the form of a time-delay. We assume that no plankton species enter or leave at the boundaries of their environments. Based on these considerations, we present the plankton model as follows:

for $t > 0, x\in\Omega = [0, l\pi]$. $u$ represents the TPP population, and $v$ represents the zooplankton population. Here, parameters $\beta, r, \theta$ and $c$ are positive. $\beta$ denotes the conversion ratio, $r$ denotes the mortality rate of zooplankton due to natural death and higher predation, $\theta$ denotes the rate of toxin liberation by TPP population, and $c$ denotes the half-saturation constant, here, we consider the case $c > 1$. $d_{1}$ and $d_{2}$ represent the diffusion coefficients of phytoplankton and zooplankton, respectively.

Here, the initial conditions of system (1.2) are considered as

where $\varphi = (\varphi_{1}, \varphi_{2})$ is uniformly continuous, and the homogeneous Neumann boundary conditions are imposed as

where $\dfrac{\partial}{\partial \nu}$ denotes the outward normal derivative on $\partial\Omega$.

This paper proceeds as follows. Section 2 gives the uniform boundedness of solutions. Section 3, by constructing the upper and lower solutions, establishes the global attractivity of positive equilibrium of the system without time-delay. Further, taking time delay as a branch parameter, we give the sufficient condition for which the system undergoes Hopf bifurcation at the interior equilibrium. In Section 4, under the assumption condition, we analyze the existence of Hopf bifurcation. Section 5 considers the global existence of these bifurcating periodic solutions. Finally, in Section 6 presents several numerical simulations.

2.   Basic properties of solutions

This section shows the boundedness of the solutions of system (1.2).

Lemma 2.1. Under the initial boundary conditions $(1.3)$–$(1.4)$, all nontrivial solutions of system $(1.2)$ are positive and uniformly bounded.

Proof. For $v(x, t)$ note that, on $0 \leq t \leq \tau$, $v(x, t-\tau) = \varphi_{2}(x, t-\tau)$ and $u(x, t)$ is bounded on $x\in\bar{\Omega}$. For further proof, we consider the following auxiliary system

where $M_{1} = \max_{x\in\bar{\Omega}, t\in[-\tau, 0]}u(x, t)$. By comparison principle, we know $v(x, t)\geq m(x, t)\geq0$ for $x\in\Omega, \; t > 0$. To proceed in this way, for $x\in\Omega, t\in[\tau, 2\tau]$, we have $v(x, t)\geq m(x, t)\geq0$. Then, by mathematical induction, we obtain the positivity of $v(x, t)$ in $x\in\Omega, \; t > 0$. Similarly, we can prove the positivity of $u(x, t)$.

Let

Then, in view of the Hopf boundary lemma and the boundary condition, it follows that

This implies

By adding the two equations in system (1.2) with the form $\beta u+v$, it follows that

By the maximum principle in Reference ^[12], this implies that

This completes the proof.

Clearly, system (1.2) has two boundary equilibria, $E_{0} = (0, 0)$ and $E_{1} = (1, 0)$.

3.   Stability of the positive equilibrium

Let us now analyze the stability of plankton system $(1.2)$ at the interior equilibrium. By referring to [1,Lemma 2.3], we know system $(1.2)$ has a unique positive equilibrium $E^{\ast}(u_{-}, v_{-})$ when the parameters satisfy the following condition

When $\tau = 0$, system (1.2) becomes

Applying the upper and lower solutions method, we show that the interior equilibrium $E^{\ast}(u_{-}, v_{-})$ for system $(3.1)$ is globally attractive under assumption (H1). We denote $E^{\ast}(u_{-}, v_{-}) = (\delta, v_{\delta})$, where

In this section, note that $(u_{1}, v_{1}) > (u_{2}, v_{2})$ means $u_{1} > u_{2}$ and $v_{1} > v_{2}$.

Theorem 3.1. Under assumption (H1), the interior equilibrium $E^{\ast}$ for the plankton system $(3.1)$ is globallyattractive.

Proof. Firstly, we can know from system $(3.1)$

by the convergence of the logistic equation, the comparison principle of parabolic equation and Lemma 2.1, then for a sufficiently small positive number $\varepsilon$, there must be a $t_{1} > 0$ such that, for $t\geq t_{1}$,

Moreover, for any initial value $(u_{0}, v_{0}) > (0, 0)$, there exists a $t_{2} > t_{1}$ such that, for $t\geq t_{2}$, $(u, v) > (\frac{\varepsilon}{2}, \frac{\varepsilon}{2}) > (0, 0)$. Then the solutions $(u, v)$ of system $(3.1)$ satisfy

Denote $(\tilde{c}_{1}, \tilde{c}_{2}) = (1+\frac{\varepsilon}{2}, \frac{\beta(1+4r)}{4r}-\frac{\varepsilon}{2})$ and $(\hat{c}_{1}, \hat{c}_{2}) = (\frac{\varepsilon}{2}, \frac{\varepsilon}{2})$, then $(0, 0) < (\hat{c}_{1}, \hat{c}_{2}) < (\tilde{c}_{1}, \tilde{c}_{2})$ and

and it is evident that there exists a $K > 0$ such that

In applying the upper and lower method in References ^[13,14], we define the two sequences $(\tilde{c}_{1}^{m}, \tilde{c}_{2}^{m})$ and $(\hat{c}_{1}^{m}, \hat{c}_{2}^{m})$ as follows:

where $(\tilde{c}_{1}^{0}, \tilde{c}_{2}^{0}) = (\tilde{c}_{1}, \tilde{c}_{2}), \; (\hat{c}_{1}^{0}, \hat{c}_{2}^{0}) = (\hat{c}_{1}, \hat{c}_{2})$.

Therefore, it follows that

then there exist $(\bar{c}_{1}, \bar{c}_{2}) > (0, 0)$ and $(\check{c}_{1}, \check{c}_{2}) > (0, 0)$ satisfying

and

Since $(\delta, v_{\delta})$ is the unique equilibrium of system (3.1), we observe $(\bar{c}_{1}, \bar{c}_{2}) = (\check{c}_{1}, \check{c}_{2})$. Thus, by employing in Reference [14,Theorem 2.2], the solutions $(u, v)$ of system (3.1) converge uniformly to $E^{\ast}$ as $t\rightarrow \infty$.

Further, for the dynamics of system (1.2), we fix the parameters $d_{1}, d_{2}, r, \theta$ and $c$, and pick delay $\tau$ as the bifurcation parameter.

Firstly, in the phase space $\mathcal{C} = C([-\tau, 0], X)$, we linearize system $(1.2)$ to analyze the stability of the positive equilibrium $(\delta, v_{\delta})$, and have

where $D = \begin{pmatrix}d_{1} & 0\\ 0&d_{2}\end{pmatrix}$ and define $L:\mathcal{C}\rightarrow X$ as

and the corresponding $L_{1}$ and $L_{2}$ are

where $A_{1}(\delta) = \frac{\delta}{1+c\delta}$ and $B_{1}(\delta) = \frac{1-\delta}{1+c\delta}$. Then by a simple derivation, the characteristic equation of the system $(1.2)$ at $(u_{-}, v_{-})$ is

where

When $\tau = 0$, Eq (3.3) becomes

and then Eq (3.5) has two roots given by

Here, we establish the definition:

Therefore, if (H1) holds, then $\beta A_{1}B_{1}-\theta\delta(1-\delta) > 0$, and adding the condition (H2), and we obtain that $D_{n}-\theta\delta(1-\delta) > 0$. Thus, we obtain the following results.

Lemma 3.2. Under the assumptions (H1) and (H2), the positive equilibrium $E^{\ast}$ of system $(1.2)$ is locally asymptotically stable when $\tau = 0$.

Remark 3.3. When the parameters of system $(1.2)$ satisfy the condition (H1) and (H2), $\lambda = 0$ is not the characteristic root of Eq $(3.5)$.

In the following, with the help from Ruan and Wei ^[15], we analyze the existence of purely imaginary eigenvalues $\lambda = \pm \text{i}\omega(\omega > 0)$ to study the stability of the positive equilibrium $E^{\ast}$. Plugging $\lambda = \text{i}\omega$ into Eq (3.3), we obtain

then it follows from separating the real and imaginary parts that

we rewrite the above Eq (3.6) by $z = \omega^{2}$, and get

and Eq (3.7) has a pair of roots given by

When assumptions (H1) and (H2) hold, we find that

Define

By the above calculation, we have

Lemma 3.4. When the parameters of system $(1.2)$ satisfy conditions (H1) and (H2), the following conclusions hold.

(i) If $2D_{n}-T_{n}^{2} < 0$ for all $n\in\mathbb{N}_{0}$, then all the roots of Eq $(3.3)$ have negative real parts.

(ii) If $2D_{n}-T_{n}^{2} > 0$ and $T_{n}^{4}-4T_{n}^{2}D_{n}+4\theta^{2}\delta^{2}(1-\delta)^{2}\geq0$, thenthe $(n+1)$th of Eq $(3.3)$ has a pair of simple pure imaginary roots $\pm \mathit{\text{i}}\omega_{n}^{+}\; (\pm \mathit{\text{i}}\omega_{n}^{-})$ at $\tau = \tau_{n, k}^{+}\; (\tau = \tau_{n, k}^{-}), \; n\in\mathfrak{N}, k\in\mathbb{N}_{0}$.

Lemma 3.5. When the parameters of system $(1.2)$ satisfy conditions (H1) and (H2), if

then

Proof. Taking the differential of both sides of Eq (3.3) with respect to $\tau$, we obtain

then we simply calculate to obtain

Completing the proof.

Thus, it is straightforward that

and we know that $\tilde{\tau} = \tau_{n, 0}^{+}$ is the smallest point changing the stability of the linearized system $(1.2)$ when the other parameters are fixed.

Theorem 3.6. Assume that (H1) and (H2) hold.

(i) If $2D_{n}-T_{n}^{2} < 0$ for all $n\in\mathbb{N}_{0}$, then the positive equilibrium $E^{\ast}$ of system $(1.2)$ is locally asymptotically stable.

(ii) If $2D_{n}-T_{n}^{2} > 0$ and $T_{n}^{4}-4T_{n}^{2}D_{n}+4\theta^{2}\delta^{2}(1-\delta)^{2} < 0$, thenthe system $(1.2)$ undergoes a Hopf bifurcation at $E^{\ast}$ for $\tau = \tau_{n, k}^{+}\; (\tau = \tau_{n, k}^{-}), \; n\in\mathfrak{N}, k\in\mathbb{N}_{0}$.Further,

4.   Analysis of local Hopf bifurcation

In this section, by employing the bifurcation theory in References ^[16,17], we analyze the property of local Hopf bifurcation near the interior equilibrium $E^{\ast}$. For fixed $k\in\mathbb{N}_{0}$ and $n\in\mathfrak{N}$, we denote $\hat{\tau} = \tau_{n, k}^{\pm}$. Setting $\tau = \hat{\tau}+\mu$, then $\mu = 0$ is the Hopf bifurcation value of system (1.2). We rescale the time by $t\rightarrow \frac{t}{\hat{\tau}}$ to normalize the delay. Meanwhile, let

System (1.2) becomes

furthermore, for simply research, system (4.1) is again transformed abstractly into

where

for $\varphi = (\varphi_{1}, \varphi_{2})^{\text{T}}\in C([-1, 0], X)$.

We know that $\pm \text{i}\omega_{n_{0}}^{\pm}\hat{\tau}$ are two pairs of simple, purely imaginary eigenvalues corresponding to the following linear system at (0, 0)

and the linear functional differential equation

By Riesz representation theorem, there exists a $2\times2$ matrix $\eta(\vartheta, \hat{\tau})\; (\vartheta\in[-1, 0])$, we know that the elements of $\eta(\vartheta, \hat{\tau})$ are of bounded variation functions such that

and we have

where $L_{1}$ and $L_{2}$ are defined in the previous section.

Here, for the sake of derivation, we give the following notation:

and

for $\varphi = (\varphi_{1}, \varphi_{1})^{\text{T}}\in\mathcal{C}$.

Define a bilinear form

where

and

Since

then

with the bilinear form

For $\varphi(\vartheta)\in C^{1}([-1, 0], \mathbb{R}^{2})$, denote

for $\psi$, define $A^{\ast}$ as

Here,

are the eigenvectors of $A^{\ast}$ and $A(0)$ corresponding to the eigenvalue $-\text{i}\omega_{n_{0}}^{+}\hat{\tau}$ and $\text{i}\omega_{n_{0}}^{+}\hat{\tau}$, respectively, where

In view of $\{\pm \text{i}\omega_{n_{0}}^{+}\hat{\tau}\}$, the center subspace of linear system (4.4) is given as

and $\mathcal{C} = P\oplus Q$, where $Q$ is the stable subspace.

In the following, when we let $\mu = 0$ in Eq (4.2), there exists a center manifold

from the theory of ^[16,17], system (4.4) can be rewritten as

then it follows that

with $\langle G, \xi_{n}\rangle = (\langle G_{1}, a_{n}\rangle, \langle G_{2}, a_{n}\rangle)^{\text{T}}$. We can also write Eq (4.9) as

with

According to the Taylor formula,

where

By simple calculation, we have

In the following, we further compute $W_{11}(\vartheta)$ and $W_{20}(\vartheta)$ to obtain $g_{21}$. Here,

where

Therefore, we have

and from Eq (4.11), obtain

and

For Eq (4.13), when $\vartheta = 0$, we have

where $E_{1} = \sum_{n = 0}^{\infty}E_{1}^{n}b_{n}$ and $E_{1} = \sum_{n = 0}^{\infty}E_{1}^{n}b_{n}$. Further, we obtain

where

We then compute $E_{1}^{n}$ and $E_{2}^{n}$ given by

where

and

Thus, by a series of derivations and calculations, we can determine the value of $g_{21}$.

Then, we analyze the bifurcation property according to the following expression:

Theorem 4.1. For system $(1.2)$, the following statements hold.

(i) If $\mu_{2} > 0 (<0)$, then the direction of the Hopf bifurcation is forward(backward), that is, there exist bifurcating periodic solutions for $\tau > \hat{\tau}(\tau < \hat{\tau})$;

(ii) If $\beta_{2} < 0 (>0)$, then the bifurcating periodic solutions on the center manifolds are orbitally stable(unstable);

(iii) If $T_{2} > 0 (<0)$, then the bifurcating periodic solutions are period increases(decreases).

Further, by Theorem 3.6, we obtain the following result.

Corollary 4.2. If $\mathit{\text{Re}}(c_{1}(0)) < 0 (>0)$, then $\mu_{2} > 0 (<0)$ and $\beta_{2} < 0 (>0)$.

5.   Analysis of global Hopf bifurcation

In the above section, the sufficient condition for the occurence of Hopf bifurcation at $E^{\ast}$ is given. In this section, we analyze the global dynamics of system (1.2) near the equilibrium $E^{\ast}$. First, we cite the global Hopf bifurcation result in Reference^[17].

Lemma 5.1. Let $S$ denote the closure of the set $\{(z, \alpha, \beta)\in E\times\mathbb{R}\times(0, \infty)\}; u(t) = z(\beta t)$ is a nontrivial $2\pi/\beta$ periodic solutions. Then for each connected component $\mathcal{C}$, at least one of the following holds:

(1) $\mathcal{C}$ is unbounded, i.e. $\sup\{\max_{t\in\mathbb{R}}|z(t)|+|\alpha|+\beta+\beta^{-1};(z, \alpha, \beta)\in\mathcal{C}\} = \infty$;

(2) $\mathcal{C}\cap(M^{\ast}\times(0, \infty))$ is finite and for all $k\geq1$, one has the equality $\Sigma_{(x_{0}, \alpha_{0}, \beta_{0})\in\mathcal{C}\cap(M^{\ast}\times(0, \infty))}\mu_{k}(x_{0}, \alpha_{0}, \beta_{0}) = 0$.

For convenience, let $z_{t} = (u_{t}, v_{t})$, the system (1.2) is rewriren as follows

where $z_{t}(\vartheta) = z(t+\vartheta)\in C([-\tau, 0], X)$, $X$ is the Banach space. With the help of the bifurcation theory in Reference ^[17], we define

and $l_{(z^{\ast}, \tau_{n_{0}, k}^{+}, \frac{2\pi}{\omega_{n_{0}}^{+}})}$ is the connected component of $(z^{\ast}, \tau_{n_{0}, k}^{+}, \frac{2\pi}{\omega_{n_{0}}^{+}})$ in ${\bf{\Sigma}}$, where $\tau_{n_{0}, k}^{+}$ and $\omega_{n_{0}}^{+}$ are defined in (3.8).

Lemma 5.2. Assume condition (H1) holds. Then system $(1.2)$ has no any nontrivial $\tau$-periodic solution.

Proof. By contradiction, assume that system (1.2) has $\tau$-periodic solution, in other words, system (3.1) then has a periodic solution. We know that system (3.1) also has the same equilibria as system (1.2), that is,

and an interior equilibrium $E^{\ast}$. For system (3.1), it is straightforward that the $u$-axis and $v$-axis are the invariable manifold, and the orbits of the system do not intersect each other, thus, the solutions can not cross the coordinate axes. Hence, it follows that it must be the interior equilibrium $E^{\ast}$ if there exists any periodic solution within the first quadrant. From the above discussion in Lemma 2.1, the positive equilibrium $E^{\ast}$ is globally attractive, then system (3.1) has no nontrivial positive periodic solution, this means that system (1.2) has no nontrivial $\tau$-periodic solution. This leads to contradiction. Completing the proof.

Theorem 5.3. Assume condition (H1) holds and $2D_{n}-T_{n}^{2} > 0$, $T_{n}^{2}-4D_{n}+4\theta^{2}\delta^{2}(1-\delta)^{2} > 0$. When $\lambda = \mathit{\text{i}}\omega_{n_{0}}^{+}$, then for each $\tau > \tau_{n_{0}, k}^{+}(k = 0, 1, \ldots)$ system $(1.2)$ has at least $k+1$ periodic solutions.

Proof. It is straightforward to calculate that the characteristic matrix of system (1.2) at a constant steady state $\tilde{E} = (\tilde{u}, \tilde{v})$ is shown by

Then we find that system (1.2) has no the center of the form $(E_{0}, \tau, p)$ and $(E_{1}, \tau, p)$.

For discussion of $E^{\ast}$, we first know that $\Big(E^{\ast}, \tau_{n_{0}, k}^{+}, \frac{2\pi}{\omega_{n_{0}}^{+}}\Big)$ is a isolated center (see [32,Definition 3.1]) by Lemma 3.4. When $\lambda(\tau_{n_{0}, k}^{+}) = \text{i}\omega_{n_{0}}^{+}$, we show the result that $\text{Re}\lambda'(\tau_{n_{0}, k}^{+}) > 0$ in Lemma 3.5. We define the smooth curve of the solution of the characteristic equation corresponding to the characteristic matrix (5.2) $\lambda:(\tau_{n_{0}, k}^{+}-\zeta, \tau_{n_{0}, k}^{+}+\zeta)\rightarrow \mathbb{C}$ satisfying as follows: there exist $\varepsilon > 0, \; \zeta > 0$ such that, for all $\tau\in[\tau_{n_{0}, k}^{+}-\zeta, \tau_{n_{0}, k}^{+}+\zeta]$

Furthermore, we define

and

Then we obtain that the crossing number of the centers $\Big(E^{\ast}, \tau_{n_{0}, j}^{+}, \frac{2\pi}{\omega_{n_{0}}^{+}}\Big)$ is

Consequently, by [32,Theorem 3.3], the connected component $l_{(z^{\ast}, \tau_{n_{0}, k}^{+}, \frac{2\pi}{\omega_{n_{0}}^{+}})}$ of $\Big(z^{\ast}, \tau_{n_{0}, k}^{+}, \frac{2\pi}{\omega_{n_{0}}^{+}}\Big)$ in ${\bf{\Sigma}}$ is unbounded. Then, from (3.3), we have

and combining with Lemmas 2.1 and 5.2, it follows that for each $\tau > \tau_{n_{0}, k}^{+}(k = 0, 1, \ldots)$ system (1.2) has at least $k+1$ periodic solutions.

6.   Simulations

In the following, let $\Omega = (0, \pi)$, and we choose the set of parameters as

to simulate the dynamics of system (1.2). For the above values and $n = 0$, we can verify that conditions (H1) and (H2) hold, $2D_{0}-T_{0}^{2}\approx0.0014 > 0$ and $T_{0}^{4}-4T_{0}^{2}D_{0}+4\theta^{2}\delta^{2}(1-\delta)^{2}\approx 0.002 > 0$, by calculation, Eq (3.8) is obtained as

then the simulations are shown as (see Figures 1 and 2)

7.   Summary and discussion

This paper's main contribution is that it provides analytic results for the reaction-diffusion TPP-zooplankton model with Holling II response function. Here, for system (3.1), when time delay is considered, it is found that time delay is a factor that causes the dynamic behavior of the system to destabilize, that is, the positive equilibrium of system (1.2) changes from globally asymptotically stable into unstable. Furthermore, from a biological point of view, plankton populations fluctuates periodically over time. Finally, numerical simulation shows that the plankton system (1.2) with time delay discussed in this paper better describes real-world problems.

Acknowledgments

The author was supported by the China Postdoctoral Science Foundation (No. 2020M681521), Postdoctoral Science Foundation of Jiangsu Province (No. 2021K456C), Postdoctoral Science Foundation of Lianyungang City (No. LYG20210012) and Open Project of Jiangsu Key Laboratory of Marine Biological Resources and Environment (No. SH20201209).

Conflict of interest

The author declare there is no conflict of interest.

Declarations

Data sharing not applicable to this article as no datasets were generated or analysed during the current study.