#### Abstract

Due to the different roles that nontoxic phytoplankton and toxin-producing phytoplankton play in the whole aquatic system, a delayed reaction-diffusion planktonic model under homogeneous Neumann boundary condition is investigated theoretically and numerically. This model describes the interactions between the zooplankton and two kinds of phytoplanktons. The long-time behavior of the model and existence of positive constant equilibrium solution are first discussed. Then, the stability of constant equilibrium solution and occurrence of Hopf bifurcation are detailed and analyzed by using the bifurcation theory. Moreover, the formulas for determining the bifurcation direction and stability of spatially bifurcating solutions are derived. Finally, some numerical simulations are performed to verify the appearance of the spatially homogeneous and nonhomogeneous periodic solutions.

#### 1. Introduction

Oceans have a major role in the global carbon cycling and so directly impact the pace and extent of climate change [1]. Marine organism can bring great economic and social values. It has the irreplaceable function in the global food processing, tourism, nutrient cycling, gas regulation, and so on. As the basis of the marine food chain system, plankton can supply food and oxygen to a myriad of marine life and can also absorb about half of the climate-warming carbon dioxide [2]. Besides, unlike fish or some intertidal creatures, plankton has rarely been commercially exploited. Moreover, this free-floating plankton can respond quickly to the temperature and change of oceanic system. As a consequence, mathematical modeling is a valuable tool for research fields of the marine ecology.

In general, plankton can be broadly divided into autotrophic phytoplankton and herbivorous zooplankton. The interactive process between phytoplankton and zooplankton is very complicated. It depends not only on the foraging style and feeding behavior, but also on other processes in the system. To describe the interactive relationship between phytoplankton and zooplankton, the following general model can be utilized:Here, and denote the concentration of phytoplankton and zooplankton respectively, and represent the growth rate and death rate of phytoplankton, respectively, is the death rate of zooplankton, and is the predation rate of zooplankton on phytoplankton. Recently, various particular cases of model (1) have been well studied [3–5]. Some interesting results about the stability and Hopf bifurcation type periodic oscillations have been obtained.

It should be noted that some phytoplankton can release toxic substances which will result in poisoning in both fish and shellfish. By accumulating in marine food webs, the toxins may have hazards on animal and human health. On the other hand, due to the rapid growth of plankton, harmful algal bloom may cause massive death of marine animals. However, most previous studies did not directly consider the toxin-producing phytoplankton population and the role of toxic phytoplankton can not be ignored. While in the establishment of planktonic mathematical models, it is necessary to include the nontoxic phytoplankton, toxin-producing phytoplankton, and zooplankton. In 2004, Chattopadhyay et al. [6] proposed the following kinetic model which has three interacting components and is with an additional factor that the release of toxic substance reduces the growth of zooplankton: where denotes the concentration of the nontoxic phytoplankton at time and and denote the concentration of toxic phytoplankton population and zooplankton at time , respectively. It is assumed that the two phytoplankton populations share the same resource. In model (2), all the coefficients are positive constants, and are the growth rates of two phytoplankton populations, respectively, is the carrying capacity, and are the maximum zooplankton ingestion rate and maximum zooplankton conversion rate, respectively, is the death rate of zooplankton, is the rate of toxin liberation by toxic phytoplankton, and is the specific predation rate of zooplankton population on toxic phytoplankton. This model shows that toxic substances released by phytoplankton have negative effects on the grazing pressure of zooplankton. It is finally concluded that toxin-producing phytoplankton may be used as a biocontrol agent for the harmful algal bloom problems. It is also mentioned that the role of time delay and environmental fluctuation in the planktonic dynamics may arouse some interesting results and needs further investigations.

Motivated by (2), some modified models have been proposed recently. For instance, Sarkar et al. [7] established a new model made of two harmful phytoplankton populations and one zooplankton population. Roy et al. [8] investigated the model where the two phytoplankton populations compete with each other. Pal et al. [9] further considered the three-component model with both nonlinear predation functions by zooplankton. Further, the study was also extended from the perspectives of stochastic dynamics and plankton-nutrient interactions, respectively [10, 11].

In view of the ocean current and monsoon, the plankton can freely drift and this spatial dispersal is subject to Fickian diffusion. So the effect of spatial diffusion has been investigated by many authors [12–17]. The results indicate that spatial diffusion has a vital role in the spatiotemporal dynamics of the planktonic model and spatial pattern may occur. Besides, the impact of time delay can not be ignored because it usually causes periodic oscillations, even chaotic behaviors, and time delay is ubiquitous in the real ecosystem [18–20].

According to the above factors, we consider the following three-component planktonic model with spatial diffusion and time delay: where , , and denote the densities of nontoxic phytoplankton, toxin-producing phytoplankton, and zooplankton at location and time , respectively, is the usual Laplace operator, denotes the depth of the water column, and the homogeneous Neumann boundary condition means that no plankton species is entering or leaving the column at the top or the bottom.

All the parameters are positive constants, , , are the three species’ diffusion rates, respectively, is the higher mortality of zooplankton, and is the time needed for zooplankton from ingesting toxic phytoplankton to dying. The other coefficients have the same meanings as in model (2). Note that the zooplankton may get eaten by higher predators, whose population is not being explicitly modelled [21–23]. So, we adopt the quadratic closure term to describe the higher mortality of zooplankton.

In this paper, we mainly investigate the spatiotemporal dynamics of delayed and diffusive system (3). The rest of the paper is organized as follows. In Section 2, the permanence and nonpersistence of system (3) are derived. In Section 3, the sufficient conditions for existence of positive constant equilibrium solution are obtained. In Section 4, the stability of equilibrium solution and delay-induced Hopf bifurcation are explored. In Section 5, the detailed formulae for determining the bifurcation properties are given by calculating the normal form on the center manifold. In Section 6, some numerical simulations are conducted to illustrate the theoretical results. Finally, some conclusions are given in Section 7.

#### 2. Long-Time Behavior

In this section, we shall show that any nonnegative solution of system (3) lies in a bounded region as , for all .

Theorem 1. *If and hold, then system (3) is permanent; that is, there exist positive constants and independent of solution such that for any nonnegative solution.*

*Proof. *From the first equation of system (3), we have The standard comparison principle implies and thus for every real number , there exists a such that , for all .

Similarly, from the second equation of (3), we have And thus for every real number , there exists a such that , for all .

The third equation of (3) can be reduced to for all . Then, we have , which indicates as . Therefore, for every real number , there exists a such that , for all .

Again from the first equation of (3), we have for all . This implies thus, Then, for every real number , there exists a such that , for all .

From the second equation of (3), we have which implies From the third equation of (3), we also have for all . Thus, Finally, if we set and then the proof is complete.

*Definition 2. *System (3) is said to be not persistent if for some of its nonnegative solutions.

Next, we discuss the nonpersistence of system (3).

Theorem 3. *If holds, then system (3) is not persistent.*

*Proof. *According to the process of Theorem 1, for an arbitrary positive constant , there exists a such that for all . Further, we have for all . Then, by the arbitrariness of , , and , we have as . The proof is complete.

#### 3. Steady State

In consideration of the biological significance of (3), we focus on the positive constant equilibrium solution. To determine the equilibrium solution of (3), we only need to solve the following algebraic equations: Simplifying the first equation of (23) and substituting it into the third equation, we have and

From (24) and the second equation of (23), we have and thus, which is equivalent to From (25) and (28), we have and whereAccording to Descartes’ rule of signs, cubic equation (31) has at least one positive real root when , that is,

For convenience, we denote any positive root of (31) by . Combining (24) and (29), we can obtain the positive solution of (23) under the conditions and .

From above analyses, we can establish the existence of positive constant equilibrium solution of (3).

Theorem 4. *If the following assumption * (*H1*)*, and ** holds, then system (3) has positive constant equilibrium solution .*

Based on the aim of this study, we always assume that condition (H1) is satisfied in the following sections.

#### 4. Hopf Bifurcation Induced by Time Delay

Here, we will regard time delay as the bifurcation parameter to investigate its effect on the stability of coexistence equilibrium solution .

Linearizing system (3) at leads to the corresponding characteristic equation: which can be simplified as the following transcendental equation: where

The special case of (35) with is

If all the roots of (37) have negative real parts for every nonnegative integer , then the positive equilibrium solution without time delay is asymptotically stable. With the help of Routh-Hurwitz criterion, is stable if and only if , , , and .

Assume that (H2), and ;

then, we have the following stability conclusion.

Theorem 5. *If assumption (H2) is satisfied, then the equilibrium solution of (3) is asymptotically stable for .*

Next, we shall discuss the distribution of characteristic roots when . Suppose is a root of (35). Then, for any nonnegative integer , we have andSeparating the real and imaginary parts results in thus, For simplicity, we set ; then, (41) can be rewritten in the form of where

If the assumption (H3) and for any

holds, then we have by combining with (H2). And the cubic equation (42) of has no positive root, so (35) has no purely imaginary root. It can be concluded that the equilibrium solution is always asymptotically stable for any and that system (3) has no spatially nonhomogeneous periodic solution.

On the other hand, when , based on Lemma 2.2 in [24], the following condition is needed to ensure the existence of positive root of (42): (H4) and ,

Denote any positive root of (42) by ; then, is a pair of purely imaginary roots of (35), where .

When , rewrite (35) as and we have Definewhere and . Then, (44) has a pair of purely imaginary roots when .

We claim that if (H5),

then

In fact, differentiating both sides of (35) with respect to , it follows that Thus, we have

From what has been discussed above and the Hopf bifurcation Theorem by Hassard et al. [25], we can draw the conclusion on the existence of spatially homogeneous Hopf bifurcation.

Theorem 6. *Suppose that conditions (H1)-(H5) are satisfied.**(i) The equilibrium solution is locally asymptotically stable for .**(ii) The equilibrium solution is unstable for and is the Hopf bifurcation value.**(ii) System (3) undergoes spatially homogeneous periodic solutions at when .*

#### 5. Stability and Direction of the Bifurcation

In this section, we investigate the properties of spatially homogeneous periodic solutions, including bifurcation direction, stability of periodic solutions, monotonicity of periodic solutions. Here, we mainly apply the normal form theory and center manifold theorem for partial functional differential equations [25, 26].

For fixed , denote bifurcation value by and introduce the new parameter ; then, is the new Hopf bifurcation value. Let and rewrite as ; system (3) can be transformed into where and are given byand

Then, the linearized system of (51) at origin is From the discussion in Section 4, we can find that characteristic equation (35) has a pair of purely imaginary roots when .

Let . Then, we consider the functional differential equation on : It is obvious that is a continuous linear function mapping into . By the Riesz representation theorem, there exists a matrix function () such that here we choose

Let be the infinitesimal generator of the semigroup induced by the solutions of (55) and be the adjoint matrix of under the bilinear pairing: where , . Then, and are a pair of adjoint operators and they both have characteristic roots . Let and be the generalized eigenspaces of and , respectively; then, is the adjoint space of and .

By direct calculation, we have the following lemma.

Lemma 7. *Let then a basis of with is and a basis a of with is *

Let , , where and From (58), we can obtain , , and further that is,

Therefore, we have

Now, we define () and construct a new basis for by . Moreover, we define , where

Let be defined by for , . Then, the center space of linear equation (54) is given by , where and , denotes the complementary subspace of .

If where and then is the infinitesimal generator induced by the solutions of (51) and (54), which can be written as the following operator differential equation: Therefore, the solution of (51) can be written in the form of where , . Specifically, the solution of (51) on the center manifold is Let , and notice that ; then, we can rewrite (73) as where . Furthermore, from [26], also satisfies where Let and

By (74), it is not difficult to compute that where , , . Let . We can get the following expressions: Since and appear in , we should also establish them. According to (77), we have and In addition, by [26], also satisfies where with and .