Stability and Bifurcation for a Delayed Diffusive Two-Zooplankton One-Phytoplankton Model with Two Different Functions

In this paper, a diffusion two-phytoplankton one-zooplankton model with time delay, Beddington–DeAnglis functional response, and Holling II functional response is proposed. First, the existence and local stability of all equilibria of such model are studied. -en, the existence of Hopf bifurcation of the corresponding model without diffusion is given by taking time delay as the bifurcation parameter. Next, the direction of Hopf bifurcation and the stability of the bifurcating periodic solutions are investigated by using the normal form theory and center manifold theorem. Furthermore, due to the local bifurcation theory of partial functional differential equations, Hopf bifurcation of the model is investigated by considering time delay as the bifurcation parameter. -e explicit formulas to determine the properties of Hopf bifurcation are given by the method of the normal form theory and center manifold theorem for partial functional differential equations. Finally, some numerical simulations are performed to check out our theoretical results.


Introduction
Generally speaking, plankton is composed of phytoplankton and zooplankton. Phytoplankton is a main producer in the marine ecosystem and the basic link of the food chain. Zooplankton is multiplying, large, and widely distributed. e zooplankton controls the number of phytoplankton by predation and plays an important regulatory role in marine ecosystems.
Plankton systems play an important role in the marine ecosystem, which have always attracted attention of ecologists and mathematicians [1][2][3]. Many researchers have established a lot of mathematical models to study the dynamic behaviors of plankton systems [4][5][6][7][8][9][10]. Phytoplankton can be divided into nontoxic phytoplankton and toxic phytoplankton (TPP). Chattopadhyay et al. [11,12] showed that toxic phytoplankton and the toxins it released affect the growth rate of zooplankton populations and the interaction between phytoplankton and zooplankton. Chattopadhyay et al. [13] proposed a delayed differential equation for toxins and pointed out that toxic phytoplankton or toxic substances may be the key reason for the termination of plankton blooms. Similar results can be found in references [14][15][16][17][18][19]. In their work [15], Saha and Bandyopadhyay considered the following TPP-zooplankton model: where P(t) and Z(t) represent the density of TPP and zooplankton, respectively, r and K denote the intrinsic growth rate and environmental capacity of TPP, μ is the maximum predation rate of zooplankton, β 1 is the ratio of biomass consumed per zooplankton for the production of new zooplankton (0 < β 1 < 1), δ is the natural death rate of zooplankton, α denotes the half saturation constant, ρ is the ratio of toxic substances produced per unit biomass of TPP, and τ represents the discrete time required for maturation of TPP to release toxic substances. ey [15] found that system (1) is locally asymptotically stable at the positive equilibrium when τ is less than its threshold; otherwise, system (1) undergoes Hopf bifurcation.
ere are many types of zooplankton in marine ecology and freshwater biological systems. Each zooplankton has its own living habit. ere is an inseparable relationship among them. erefore, the impact of different species on the ecosystem is the concern for most scholars. Some researchers considered two types of zooplankton [20][21][22][23][24]. Lv et al. [22] investigated the complex dynamics of the two zooplankton-phytoplankton model with delay: where P(t), Z 1 (t), and Z 2 (t) represent the density of TPP, zooplankton 1, and zooplankton 2 at time t, respectively. e parameter μ i is the maximum uptake rate of zooplankton i(i � 1 and 2). β i denotes the ratio of biomass conversion of zooplankton i(i � 1 and 2). ρ i is the ratio of toxic substances produced by per unit biomass of TPP for the zooplankton i(i � 1 and 2). α i represents the half saturation constant of zooplankton i(i � 1 and 2). g i is the interspecific competition coefficient of zooplankton i(i � 1 and 2). d i denotes the natural mortality of corresponding zooplankton i(i � 1 and 2). e biological significance of other parameters is the same as system (1). All parameters are positive constants. System (2) shows that the intrinsic characteristics of toxins, such as the release rate of toxicity and the delay of toxin production, does not irreversibly change the stability of this system. It is finally concluded that toxin-producing phytoplankton may be used as a biocontrol agent for the harmful algal bloom problems.
In nature, species interact with each other to produce different biological effects. Each effect can vividly describe the relationship between species. Beddington [25] and DeAngelis et al. [26] proposed Beddington-DeAnglis functional response (B-D function), which is similar to the functional response of the famous Holling type II. e B-D function can better reflect the prey-dependent and the mutual interference between predators [27][28][29][30][31][32][33]. Cantrell and Cosner [27] investigated the predator-prey model with the B-D function and showed that the degree of mutual interference between predators affects the position and stability of the equilibria. Meng and Wang [31] established a delayed diffusive model with the B-D function and discussed the existence of Hopf bifurcation and its properties.
From a biological point of view, individual organisms are distributed in space, usually interacting with the physical environment and other organisms in their spatial neighborhood. Generally speaking, self-diffusion refers to the transmission from an area of higher concentration to a low concentration area through random motion. e diffusion of individuals may be related to other things, such as finding food and escaping from high risk of infection. Hardy [34] pointed out that the spatial distribution of marine plankton is not uniform. e cooperativity of plankton is weak, but plankton can move freely under the influence of ocean currents and monsoons.
is spatial diffusion is subject to Fick's law. erefore, the influence of spatial diffusion on the phytoplankton-plankton model has been paid more attention by many scholars [35][36][37][38][39][40][41][42][43][44][45]. In the work by Jia et al. [46], a three-component plankton model with spatial diffusion and time delay is proposed, which describes the relationship between a zooplankton and two phytoplankton. Rao [47] studied the complex dynamics of a spatial toxicphytoplankton-zooplankton model with the Holling II function and showed that the interaction between toxic phytoplankton and zooplankton in the marine environment may be partially driven by diffusivity or environmental carrying capacity.
Encouraged by the above work, we will consider the Beddington-DeAngelis functional response and the effect of self-diffusion for phytoplankton and zooplankton into system (2) in this paper.
us, a diffusion toxic phytoplankton-zooplankton model with time delay and two kinds of functional responses is given as follows: 2 Complexity where P(t, x), Z 1 (t, x), and Z 2 (t, x) represent the densities of the toxic-phytoplankton and two kinds of zooplankton at location x and time t, respectively. Δ is the Laplace operator, and the homogeneous Neumann boundary condition means that no plankton species enters or leaves this area. d 1 , d 2 , and d 3 are the diffusion rates of the three species, g i is expressed as natural mortality of zooplankton i(i � 1 and 2), and s is the measure of the degree of mutual interference among zooplankton Z 1 . e remaining parameters have the same biological meanings as model (2). All parameters are positive constants. e remain parts of this paper are organized as follows. e existence and stability of equilibria are discussed in Section 2. e existence and properties of Hopf bifurcation of model (2) with time delay by the normal form theory and the center manifold theorem are given in Section 3. In Section 4, we analyze the Hopf bifurcation of model (3) at the coexistence equilibrium and give the explicit formulas to determine the properties of Hopf bifurcation by using the normal form theory and center manifold theorem for partial functional differential equations (PFDEs). In Section 5, numerical simulations are given to verify the theoretical results. A discussion is given to conclude this work in Section 6.

Existence and Stability of Equilibria of ODE
In this section, we investigate the dynamic behavior of system (3) without delay and diffusion. e corresponding ordinary differential model of system (3) is as follows: e Jacobian matrix of system (4) is

e Trivial Equilibriums.
It is obvious that system (4) has the trivial equilibrium E 0 (0, 0, 0) and the semitrivial equilibrium E 1 (K, 0, 0). e characteristic equation of system (4) at equilibrium E 0 is erefore, we can obtain three eigenvalues: us, the trivial equilibrium E 0 is unstable. e local stability of system (4) at the semitrivial equilibrium E 1 (K, 0, 0) is as follows.
us, equation (13) has two negative roots. en, we have the following conclusion.
Similarly, the other boundary equilibrium E 01 (P, 0, Z 2 ) exists when the assumption (H 3 ): Next, we investigate stability of the boundary equilibrium E 01 . e method is the same to that of the boundary equilibrium E 10 . e characteristic equation of system (4) at the boundary equilibrium E 01 is where If λ 1 � B 22 , then λ 2 and λ 3 are the roots of the following equation:

e Coexistence Equilibrium.
In biology, we are interested in the existence and stability of the coexistence equilibrium. Next, we will analyze the existence of the coexistence equilibrium of system (4). System (4) has a coexistence equilibrium E * (P * , Z * 1 , Z * 2 ), where P * , Z * 1 , and Z * 2 are the positive solution of the following equations: From the third equation of (17), we can get that P * � By substituting P * and Z * 1 into the first equation of (17), we can get that From a biological point of view, the stable coexistence of multiple species is very meaningful. Now we discuss the stability of the positive equilibrium E * . e characteristic equation of system (4) at E * is where Complexity By the Routh-Hurwitz criterion [48], E * is locally as- . en, we have the following conclusion.

Theorem 4.
Under the assumption (H 5 ), the coexistence equilibrium E * of system (4) is locally asymptotically stable if the assumption (H 6 ) holds. Now, we consider the parameter s as a bifurcation parameter to study the existence of Hopf bifurcation. According to Liu [49], we can get that the following result.

e Existence of Bifurcation.
In this section, under the assumption (H 5 ), we regard time delay τ as the bifurcation parameter to study its influence on the stability of coexistence equilibrium E * (P * , Z * 1 , Z * 2 ). System (4) with time delay is as follows: Let (21), we can get the following linearized system: where L 1 � a 11 a 12 a 13 a 21 a 22 0 where a 1k � C 1k (k � 1, 2, 3) and erefore, the characteristic equation of system (22) is where When τ � 0, the distribution of roots of equation (25) has been discussed above. Next, we will discuss the distribution of characteristic roots of equation (25) when τ > 0. First, we suppose that iω 1 (ω 1 > 0) is a root of equation (25).
en, we obtain that Separating the real and imaginary parts of equation (27), we have From equation (28), we get that 6 Complexity Let (29) can be rewritten as To obtain the existence of Hopf bifurcation, equation (30) has at least one positive root. Due to (1) If the assumption (H 7 ) holds, equation (30) has no positive root (2) If the assumption (H 8 ) holds, equation (30) has at least one positive root Without losing of generality, we assume that equation (30) has three positive roots, denoted by y 1 , y 2 , and y 3 . Hence, equation (29) has three positive roots ω k 1 � �� y k √ (k � 1, 2, and 3). By equation (28), we can get that the critical values of time delay τ are for k � 1, 2, and 3, j ∈ N 0 . erefore, ±iω k 1 is a pair of purely imaginary roots of equation (25) when τ � τ j k (k � 1, 2, and 3; j � 1, 2, . . .). Denote , and 3 and j ∈ N 0 . en, we have the following transversality condition.

e Properties of Hopf Bifurcation.
In the previous parts, we have discussed the existence of Hopf bifurcation of model (21) at the coexistence equilibrium E * . Next, we will study the direction of Hopf bifurcation and the stability of bifurcated periodic solution by using the normal form theory and center manifold theory developed by Hassard et al. [51].
Under the condition (H 6 ), h 1n + c 1n > 0, h 0n + c 0n > 0, and h 2n (h 1n + c 1n ) − h 0n − c 0n > 0, for all n ∈ N 0 . rough the Routh-Hurwitz criterion [48], all the roots of equation (55) have negative real parts. From eorems 5.1.1 and 5.1.3 of Henry [54], system (3) at E * is locally asymptotically stable when τ � 0. To sum up, diffusion does not change the stability of system (3) at equilibrium E * . en, we can get the following theorem. Theorem 8. Assume the conditions (H 5 ) and (H 6 ) hold, then the coexistence equilibrium E * (P, Z 1 , Z 2 ) of system (3) with τ � 0 is locally asymptotically stable. Now, we study the effect of the delay τ ≠ 0 on the stability of the positive equilibrium E * of system (3). First, we seek critical values of τ such that equation (53) has a pair of purely imaginary roots.
From Lemmas 2.2 and 4.2 in [50], we obtain the following lemma. For polynomial equation (60), we can get the following conclusions about the distribution of its positive roots.
To sum up, system (3) without diffusion undergoes a spatially homogeneous Hopf bifurcation at the coexistence equilibrium E * . System (3) undergoes a spatially inhomogeneous Hopf bifurcation at the coexistence equilibrium E * .

Direction and Stability of Bifurcated Periodic Solutions.
We have get some conditions, such that a spatially homogeneous or inhomogeneous periodic solutions bifurcate from the coexistence equilibrium E * of system (3) when τ crosses through the critical values τ n k,j (k � 1, 2, and 3 and j ∈ N 0 ), n � 0 and n 0 . In this part, we will study the direction of Hopf bifurcation and the stability of bifurcated periodic solutions by utilizing the normal form theory and center manifold theorem for partial functional differential equations [55]. For convenience, we fix k � 1, 2, and 3 and j � 0, 1, 2, . . ., and denote τ n k,j (n � 0 and n 0 ) by τ and ω 2 k,n (n � 0 and n 0 ) by ω n .
Let ℘: � C 1 ([− 1, 0], R 3 ). Next, we consider the functional differential equation on ℘: Obviously, Γ(τ) is continuous linear function mapping According to the Riesz representation theorem [52], there exists a 3 × 3 matrix function η 1 (θ, τ, n)(− 1 < θ < 0), whose entries are bounded variation that where R 3 * is the three-dimensional vector space of row vectors. A τ denotes the infinitesimal generator of semigroup induced by the solutions of equation (72), and A * τ is the adjoint matrix of A τ under the bilinear pairing: where ϕ ∈ C 1 , φ ∈ C 1 ([− 1, 0], R 3 * ). us, A τ and A * τ are a pair of adjoint operators, and their eigenvalues are Υ n � ± iω n τ . Let F and F * be the center space, that is, the generalized characteristic space of A τ and A * τ , respectively. en, F * is adjoint space of F, and dim F � dim F * . us, we have the following lemma.

Complexity
Data Availability e data used to support the findings of our study are included within the article. You can find the corresponding values of some parameters in the beginning of the "Numerical Simulation," which are taken from the numerical simulation section in Reference [22]. Of course, every reader can freely access the parameter values and data supporting the conclusions of the study in [22].

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.