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 .

Consequently, (74) and (81)-(83) lead to Since has a pair of purely imaginary characteristic roots, (83) has the unique solution such that For , it follows from (84) that Thus, for , we have and Combining with the definition of and (86), we have Since , we have where From the definition of and (86) and (92), we have Noticing that we have From previous formulas, we can obtain

With the same method, we also have and where At this point, we are able to completely establish the value of and thereby can determine the properties of Hopf bifurcation.

System (3) has the following PoincarƩ normal form: where Hence, By the Hopf bifurcation theory [25], we know that determines the bifurcation direction: if (), then the Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for (); determines the stability of bifurcating periodic solutions: if (), then the periodic solutions are stable (unstable); determines the monotonicity of the period of periodic solutions: if (), then the period increases (decreases).

6. Numerical Simulation

In this section, we conduct the numerical simulations with the help of MATLAB. We first choose the following parameter value: By direct calculation, we have that system (3) has the unique positive equilibrium solution and the corresponding Hopf bifurcation value is .

Figures 1ā€“3 show that the constant equilibrium solution of (3) is asymptotically stable when the time delay is zero or appropriately small. On the other hand, once the time delay is larger than the critical value , the equilibrium solution would no longer be stable and spatially homogeneous periodic solution will bifurcate at the equilibrium solution (see Figure 4).

We reselect , , , and keep other coefficients the same. From Figures 5-6, it is shown that the equilibrium solution is asymptotically stable when and spatially inhomogeneous periodic solution exists when . It confirms that both spatial diffusion and toxin delay have significant effects on the spatiotemporal dynamics of system (3). Besides, if we set , then the spatially inhomogeneous periodic solution vanishes and the equilibrium solution becomes asymptotically stable (see Figure 7). It can be concluded that spatially inhomogeneous pattern is more prone to occurring in small space.

7. Conclusions

In this paper, we have proposed a delayed reaction-diffusion model incorporating three plankton populations. Mathematical analysis indicates that the three populations can coexist when the natural growth rate of nontoxic phytoplankton is large and the death rate of zooplankton by toxin is small. In this case, the balance is finally achieved by interdependence and mutual restraint. Otherwise, the zooplankton will become extinct if the natural growth rate of nontoxic phytoplankton is small and the death rate of zooplankton by toxin is large. In the latter case, the biomasses of nontoxic phytoplankton and toxin-producing phytoplankton reach saturation and algae bloom occurs.

And when considering the toxin delay, the spatiotemporal dynamics of system is almost unaffected when the time delay is sufficiently small. However, when the time delay passes through some critical value, spatially homogeneous or inhomogeneous periodic solution may arise. This means that algal bloom erupts periodically under certain conditions. Therefore, some measures can be adopted to control or defer the occurrence of Hopf bifurcation, such as reducing the toxin delay or adding feedback control. Numerical simulations complementally illustrate that when there is no or only tiny time delay, the population distribution pattern is eventually spatially homogeneous even if the initial distribution is inhomogeneous. When the time delay is sufficiently large, the distribution pattern is time periodic. Moreover, nontoxic phytoplankton has the largest oscillation amplitude, toxin-producing phytoplankton has the smallest oscillation amplitude, and the zooplankton has modest amplitude. These phenomena show that toxin delay has the greatest effects on nontoxic phytoplankton.

Data Availability

No data were used to support this study.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This research was supported by the National Natural Science Foundation of China (61703001), the MOE Project of Humanities and Social Sciences of China (17YJC630175), the Natural Science Research Projects of Universities in Anhui Province (KJ2017ZD35, KJ2017A432), and Anhui Provincial Natural Science Foundation (1508085MA09).