#### Abstract

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. Then, 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. The 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.

#### 1. 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. The 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–3]. Many researchers have established a lot of mathematical models to study the dynamic behaviors of plankton systems [4–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–19]. In their work [15], Saha and Bandyopadhyay considered the following TPP-zooplankton model:where and represent the density of TPP and zooplankton, respectively, and denote the intrinsic growth rate and environmental capacity of TPP, is the maximum predation rate of zooplankton, is the ratio of biomass consumed per zooplankton for the production of new zooplankton , 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. They [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.

There are many types of zooplankton in marine ecology and freshwater biological systems. Each zooplankton has its own living habit. There is an inseparable relationship among them. Therefore, the impact of different species on the ecosystem is the concern for most scholars. Some researchers considered two types of zooplankton [20–24]. Lv et al. [22] investigated the complex dynamics of the two zooplankton-phytoplankton model with delay:where , , and represent the density of TPP, zooplankton 1, and zooplankton 2 at time , respectively. The parameter is the maximum uptake rate of zooplankton . denotes the ratio of biomass conversion of zooplankton . is the ratio of toxic substances produced by per unit biomass of TPP for the zooplankton . represents the half saturation constant of zooplankton . is the interspecific competition coefficient of zooplankton . denotes the natural mortality of corresponding zooplankton . The 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. The B-D function can better reflect the prey-dependent and the mutual interference between predators [27–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. The 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. The cooperativity of plankton is weak, but plankton can move freely under the influence of ocean currents and monsoons. This spatial diffusion is subject to Fick’s law. Therefore, the influence of spatial diffusion on the phytoplankton-plankton model has been paid more attention by many scholars [35–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 toxic-phytoplankton-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. Thus, a diffusion toxic phytoplankton-zooplankton model with time delay and two kinds of functional responses is given as follows:where , and represent the densities of the toxic-phytoplankton and two kinds of zooplankton at location and time , respectively. is the Laplace operator, and the homogeneous Neumann boundary condition means that no plankton species enters or leaves this area. , and are the diffusion rates of the three species, is expressed as natural mortality of zooplankton , and is the measure of the degree of mutual interference among zooplankton . The remaining parameters have the same biological meanings as model (2). All parameters are positive constants.

The remain parts of this paper are organized as follows. The existence and stability of equilibria are discussed in Section 2. The 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.

#### 2. Existence and Stability of Equilibria of ODE

In this section, we investigate the dynamic behavior of system (3) without delay and diffusion. The corresponding ordinary differential model of system (3) is as follows:

The Jacobian matrix of system (4) is

##### 2.1. The Trivial Equilibriums

It is obvious that system (4) has the trivial equilibrium and the semitrivial equilibrium .

The characteristic equation of system (4) at equilibrium is

Therefore, we can obtain three eigenvalues:

Thus, the trivial equilibrium is unstable.

The local stability of system (4) at the semitrivial equilibrium is as follows.

We can obtain that the characteristic equation of system (4) at the semitrivial equilibrium is as follows:

From equation (8), we can get three eigenvalues:

So and if and only if . Thus, the semitrivial equilibrium is locally asymptotically stable; otherwise, is unstable.

Theorem 1. *System (4) is unstable at the trivial equilibrium . System (4) is locally asymptotically stable at the semitrivial equilibrium if , but system (4) is unstable if .*

##### 2.2. The Boundary Equilibriums

System (4) has two boundary equilibriums. From the second equation of system (4), we can obtain that . Then, by substituting it into the first equation of system (4), we can obtain that the following quadratic equationwhere .

If the assumption , then and ; thus, equation (10) has at least one positive root, defined as . From the biological meaning, is nonnegative when . If the assumption holds, then system (4) has the boundary equilibrium .

Now, under the assumption , we prove the stability of the boundary equilibrium . Linearizing system (4) at , we get the characteristic equation:where

Further , and are the roots of the quadratic equation:where . Let ; we assume that : , , and . Under the assumption , , and . Thus, equation (13) has two negative roots. Then, we have the following conclusion.

Theorem 2. *The boundary equilibrium exists if holds. Furthermore, of system (4) is locally asymptotically stable if assumptions and hold.*

Similarly, the other boundary equilibrium exists when the assumption : holds, where and .

Next, we investigate stability of the boundary equilibrium . The method is the same to that of the boundary equilibrium .

The characteristic equation of system (4) at the boundary equilibrium iswhere

If , then and are the roots of the following equation:

Let ; we assume that : and . Under the assumption , , , and ; thus, equation (16) has two real negative roots. Then, we have the following result.

Theorem 3. *If assumption is true, the boundary equilibrium of system (4) exists. Further, if assumptions and are true, then the boundary equilibrium of system (4) is locally asymptotically stable.*

##### 2.3. The 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 , where , , and are the positive solution of the following equations:

From the third equation of (17), we can get that if holds. By substituting into the second equation of (17), we can get that if , where . By substituting and into the first equation of (17), we can get that if holds, where . Thus, system (4) has one positive equilibrium when the condition : , , and is true.

From a biological point of view, the stable coexistence of multiple species is very meaningful. Now we discuss the stability of the positive equilibrium . The characteristic equation of system (4) at iswhere

By the Routh–Hurwitz criterion [48], is locally asymptotically stable if , , and . Hence, we assume that : . Then, we have the following conclusion.

Theorem 4. *Under the assumption , the coexistence equilibrium of system (4) is locally asymptotically stable if the assumption holds.*

Now, we consider the parameter as a bifurcation parameter to study the existence of Hopf bifurcation. According to Liu [49], we can get that the following result.

Theorem 5. *If the characteristic equation of the coexistence equilibrium is given bywhere , and are smooth function of in an open interval about such that*(1)*, , and *(2)*; then, a simple Hopf bifurcation occurs at **Proof. In fact, , , and . It is easy to see that there exists such that , , , , , and . So, we can get . Thus, the Hopf bifurcation occurs at .□*

#### 3. Hopf Bifurcation of System (3) without Diffusion

##### 3.1. The Existence of Bifurcation

In this section, under the assumption , we regard time delay as the bifurcation parameter to study its influence on the stability of coexistence equilibrium . System (4) with time delay is as follows:

Let . Substituting into system (21), we can get the following linearized system:wherewhere and

Therefore, the characteristic equation of system (22) iswhere

When , the distribution of roots of equation (25) has been discussed above. Next, we will discuss the distribution of characteristic roots of equation (25) when . First, we suppose that is a root of equation (25). Then, we obtain that

Separating the real and imaginary parts of equation (27), we have

From equation (28), we get that

Let , , , and . Then, equation (29) can be rewritten as

To obtain the existence of Hopf bifurcation, equation (30) has at least one positive root. Due to , , and . Hence, we assume that : ; : , , and .

From Lemmas 2.2 and 4.2 in [50], we obtain the following lemma.

Lemma 1. *For polynomial equation (30), the following conclusions are true.*(1)*If the assumption holds, equation (30) has no positive root*(2)*If the assumption holds, equation (30) has at least one positive root*

Without losing of generality, we assume that equation (30) has three positive roots, denoted by , , and . Hence, equation (29) has three positive roots . By equation (28), we can get that the critical values of time delay arefor . Therefore, is a pair of purely imaginary roots of equation (25) when . Denote and .

Let be the root of equation (25) near with and . Then, we have the following transversality condition.

Lemma 2. *Suppose that ; then, the following transversality condition is satisfied: .*

*Proof. *Differentiating the two sides of equation (25) with respect to , we have

Then, we obtain

By some calculation, we getwhere .

Hence,

By the theory of Hopf bifurcation, we have the following result.

Theorem 6. *Assume that and hold, then the following statements are true.*(1)*If the hypothesis holds, then the coexistence equilibrium of model (21) is locally asymptotically stable for .*(2)*If the hypothesis and hold, then the coexistence equilibrium of model (21) is locally asymptotically stable for and unstable for . Moreover, model (21) undergoes Hopf bifurcation at the coexistence equilibrium when .*

##### 3.2. The Properties of Hopf Bifurcation

In the previous parts, we have discussed the existence of Hopf bifurcation of model (21) at the coexistence equilibrium . 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]. Let with ; then, is a Hopf bifurcation value of model (21). Let , , and ; then, model (21) can be rewritten in the phase space :where has the following:whereand and are given in (23).

According to the Riesz representation theorem [52], there exists a function whose elements are of bounded variance such that

Thus, we chose

For , we define the operator:

Then, system (36) is equivalent to

For , is defined as the adjoint operator of as

The bilinear product is given bywhere .

By the discussion in Section 3.1, are the eigenvalues of and , which satisfy and . Suppose that is the eigenvector of corresponding to the eigenvalue and is the eigenvector of corresponding to the eigenvalue . By computation, we have

From equation (44) and , we can get

Thus,

In the remainder of this section, through using the method in [51], we can get the following relevant parameters which can determine the direction of Hopf bifurcation and the stability of the bifurcation periodic solutions:where