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. 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.
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.  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 , 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  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.  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  and DeAngelis et al.  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  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  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  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. , a three-component plankton model with spatial diffusion and time delay is proposed, which describes the relationship between a zooplankton and two phytoplankton. Rao  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
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.
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.
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 , 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 , 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 , 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 .
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. . 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 , 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
In the remainder of this section, through using the method in , we can get the following relevant parameters which can determine the direction of Hopf bifurcation and the stability of the bifurcation periodic solutions:where