Journal of Applied Mathematics

Volume 2013 (2013), Article ID 608073, 10 pages

http://dx.doi.org/10.1155/2013/608073

## Dynamical Complexity of a Spatial Phytoplankton-Zooplankton Model with an Alternative Prey and Refuge Effect

^{1}School of Mathematics and Information Science, Wenzhou University, Wenzhou, Zhejiang 325035, China^{2}School of Life and Environmental Science, Wenzhou University, Wenzhou, Zhejiang 325035, China

Received 23 January 2013; Accepted 29 March 2013

Academic Editor: Wan-Tong Li

Copyright © 2013 Weiwei Zhang and Min Zhao. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

The spatiotemporal dynamics of a phytoplankton-zooplankton model with an alternative prey and refuge effect is investigated mathematically and numerically. The stability of the equilibrium point and the traveling wave solution of the phytoplankton-zooplankton model are described based on theoretical mathematical work, which provides the basis of the numerical simulation. The numerical analysis shows that refuges have a strong effect on the spatiotemporal dynamics of the model according to the pattern formation. These results may help us to understand prey-predator interactions in water ecosystems. They are also relevant to research into phytoplankton-zooplankton ecosystems.

#### 1. Introduction

In recent years, the degradation of water ecosystems has exacerbated certain algal species in aquatic systems [1]. Increasing attention is being paid to preventing this exacerbation from reaching a critical condition because of adverse impacts on fisheries and aquaculture [1]. In aquatic systems, there is an ecological threshold below which species remains essentially undetected, whereas they seize the opportunity for unchecked growth if the ecological threshold is exceeded [2]. Predator-prey models have been studied widely and it is well known that these model can directly reflect changes in the sizes of populations, which can be used to determine the ecological threshold for certain algal species in aquatic systems.

Since the pioneering work of Lotka and Volterra, the predator-prey model has developed significantly in mathematical ecology [3]. A variety of models have been studied including the Wangersky-Cunningham model, the classical prey-dependent predator-prey model, and the ratio-dependent function response model [4–13]. A previous model studied the global dynamics of a predator-prey model with a nonconstant death rate and diffusion [14]. This model incorporated Holling type II and Leslie-Gower functional responses [15, 16]. Most authors have focused their attention on a single prey item known as the focal prey in these studies, although some studies have show that the presence of alternative food sources can affect biological control via a variety of mechanisms [17]. The population of a prey species may be affected negatively by another prey species because they may lead to higher predation rates for both prey items if the prey shares the predator [17]. However, alternative prey can also offset predation on the focal prey if the predator has a predilection for the alternative prey [17, 18]. In this case, the alternative prey may have a positive effect on the population of the focal prey [17].

In a predator-prey model, there are solution states in some cases. However, the transformation between these states is important. Thus, many researchers have studied the presence of traveling wave solutions in predator-prey models [19–25]. In the previous works, we mainly discussed the dynamics of prey-predator system with impulsive control strategy [3, 26–28].

In this paper, we discuss a two-species predator-prey model and the effect of alternative prey. This paper is organized as follows. In Section 2, we introduce the predator-prey model and the reaction-diffusion term. Section 3 presents the mathematical analysis, including the existence of an equilibrium point and boundedness, a stability analysis of the equilibrium point, and traveling wave solutions. In Section 4, we describe a numerical simulation of our analytical results, which helps us to understand the feasibility of the theorem. The final section provides our conclusions.

#### 2. The Model

First, we introduce the following model: where and are the densities of phytoplankton and zooplankton at time and position , respectively. The parameters of the model (1) can be expressed as follows: represents the inherent growth rate of the phytoplankton without any environment limitations; is the carrying capacity of the phytoplankton in the presence of predators and harvesting; represents the consumption rate; is the half-saturation constant that determines how quickly the maximum is attained when the phytoplankton density increases; is the conversion factor denoting the number of new herbivores produced for each capture; is the growth rate of micro-grazers due to the alternative prey; denotes the food-independent zooplankton mortality rate where ; is the refuge protect the prey, where is constant [29–32].

For simplicity, using the scaling , , and , we have where , , , and . In model (2), we state that the phytoplankton and zooplankton generally experience the same homogeneous environment. Based on a large number of experimental observations in actual environments, it is known that the distribution of individual organisms often depends on interactions with the physical environment and other organisms. A growing body of research indicates that space can change the dynamics of populations and the structure of communities [33]. Thus, we assume that the phytoplankton and zooplankton move randomly, which is described using random Brownian motion [34, 35]. The system can be described using the following set of differential equations: where is the usual Laplacian operator in two-dimensional space. It is assumed that the environment is uniform and all parameters of the model (3) are independent of space or time, that is, these parameters are constant. Parameters and are positive diffusion coefficients for the phytoplankton and zooplankton, respectively. From the viewpoint of biology, we are only interested in the dynamics of models (2) and (3) in the first quadrant . Thus, model (3) is analyzed using the following initial conditions: and the homogeneous Neumann boundary condition: which indicates that there are no population fluxes through the boundary and that no external input is imposed from the outside [19]. is the outward unit normal vector for the boundary , which is usually assumed to be smooth [35].

#### 3. Mathematical Analysis

First, we will prove that all solutions of model (3) are bounded.

Theorem 1. *Let be any solution of model (3). Thus, we can find
*

*Proof. *Based on the first equation in model (3), we derive

Based on the comparison principle, let be the solution of the equation

Thus, we can obtain . For any , there is a , such that for any . Therefore, is defined for all and .

Next, we will discuss the boundedness of based on the second equation of model (3). We let and we find that

Assuming that there is a such that
we can determine , for all , , where .

Based on the formulae above and the first equation in model (3), we can obtain
so we find that
where is the solution of the equation and .

Then, as , . Thus, there is a such that

Substituting inequality (13) into the second equation of model (3), we can get

Based on the comparison principle, we can find that as , , although this is clearly contradictory.

Let , , . If there is a , such that , for all , , then it is considered to be bounded and obviously we can find that for for all , .

If for a sequence of tending to as , we can say that it is unbounded, that is, the -axis is divided by the zeroes into a sequence of intervals , . Given this assumption, , for all , . For the intervals , arguing the inequality (10) is not possible, so we find that if , for all , , .

Based on the comparison principle, for , where is the solution of the equation given suitable initial conditions. Thus,

Therefore, we can find .

Among the conditions that satisfy Theorem 1, we mainly address the steady states and their stabilities for (3), as follows.

(i), which corresponds to the extinction of the phytoplankton and zooplankton, which is apparently an unstable saddle point. (ii), which corresponds to the extinction of the zooplankton.(iii)Their coexistence will be observed in different states and , where
and . Obviously, if , and are negative. In this study, we assume that is satisfied. and have similar properties, thus we focus only on the positive equilibrium point .

In an actual environment with changing parameters, we will observe different states. In this study, the main aim is to address changes in the state of the equilibrium points . Next, we determine the stability of the equilibrium points for model (3). Substituting
into model (3), where , and selecting all terms that are linear about ,
where
where
Substituting (19) into (18), we obtain:

##### 3.1. Stability Analysis

Theorem 2. * (i) is unstable;** (ii) if , then for model (3) is asymptotically stable;** (iii) is locally asymptotically stable when and *(a)*if , then ,*(b)*if , then .*

*Proof. *The stability of the equilibrium state is determined based on the nature of the eigenvalues in the Jacobian matrix . Substituting into (19), we can determine the eigenvalues of the Jacobian matrix , and . Apparently , so is unstable. The stability of can be established in a similar manner as , so we only discuss the stability of .

We can use to determine the eigenvalues of the operator on , which satisfy the homogeneous Neumann boundary condition and . is the eigenvector relative to , so

By substituting into (18) and (19), and expanding the solution of (18),

Substituting (23) into (18) and balancing the coefficients about each , we can obtain where the matrix satisfies .

Thus, is stable if and only if each decays to zero, that is, there are two eigenvalues with negative real parts for each matrix . Based on simple algebraic computations, the characteristic polynomial of is determined as follow,

If and is satisfied, the above conditions are met. Therefore, with a random , and hold if and .

Thus, we can obtain: (a)if , then ,(b)if , then . This completes the proof.

Theorem 3. *If , and hold, the steady state of the model (3) is globally asymptotically stable, where .*

*Proof. *From La Salle’s theorem [36], the Lyapunov function of model (3) needs to be analyzed to prove the theorem. Thus, the following function is structured,
assuming and are satisfied.
is positive for all in the positive quadrant, except for where . If is proven, then is the Lyapunov function. Calculating the rate of change of for the solution of model (3), we find that,

Based on simple algebraic computations, becomes
where and .

Obviously, increases strictly monotonically, so is negative for all in the positive quadrant if and only if decreases strictly monotonically, and we find that and holds.

In order to determine , we can use the Green formula. If it satisfies the homogeneous Neumann boundary conditions, we find
where . Based on Theorem 1, we know . If is found, . If , and hold, so is negative. Therefore, is globally asymptotically stable.

##### 3.2. Traveling Wave Solutions

In the previous section, we discussed the stability of the equilibrium points for model (3). It can be got that is unstable, is unstable for model (3), and the steady state of the model (3) is globally asymptotically stable in suitable conditions. However, this is an essential condition that leads to traveling wave solutions. In this section, the main aim is to find a heteroclinic orbit between and or and .

To find the heteroclinic orbit between and , we need to show that model (3) has a solution with the special form , , where the wave speed parameter is positive. After substituting , , and into (3), we get the following wave equations

If and are nonnegative and they satisfy the following boundary conditions which shows there is a traveling wave solution between and .

Substituting and into (30),

After computing the Jacobian matrix for model (32) at , we get the following matrix:

Therefore, the characteristic equation of is given by , and are its eigenvalues, as follows

If , then and are a pair of complex conjugate eigenvalues with positive real parts. Based on the theorem given in [37], we can find a two-dimensional unstable manifold, which has the base at . If , the trajectory approaching must have for some . However, this is contradictory if .

This discussion shows that because , traveling wave solutions cannot be found for model (32), whereas if , model (32) has non-negative solutions that satisfy the boundary conditions. In some suitable conditions, therefore, there is a traveling wave solution from to . Next, we discuss the Jacobian matrix of model (32) at , where , , and . Thus, the characteristic equation of is given by and is its eigenvalues, as follows: .

Based on the Routh-Hurwitz stability criterion, when and , the characteristic equation of always has two eigenvalues with positive real parts and two eigenvalues with negative real parts, that is, to say, there is a stable two-dimensional manifold at and a trajectory from to where . Thus, if and hold, there is a traveling wave solution.

#### 4. Numerical Results

In this section, the spatiotemporal dynamics of the phytoplankton-zooplankton model will be simulated, to better understand the theoretical findings. All our numerical simulations use the homogeneous Neumann boundary conditions where . Based on numerical simulations, we found that the spatial pattern distributions of the phytoplankton and zooplankton were always of the same type. This is because the Holling II function reaction term shows that the number of predatory zooplankton is proportional to the number of phytoplankton, while the number of zooplankton is equal to the number of predators. Therefore, we only show the spatial patterns of the phytoplankton and the snapshot with blue parts represents the low density value of phytoplankton , whereas the snapshot with red parts represents the high density value of phytoplankton . To investigate the interplay among several critical parameters, we determine the pattern formation using model (3) where other parameters were fixed as , , , , and .

We only changed the value of to test its effects on the spatiotemporal dynamics of the model (3). In Figure 1, we show some snapshots of the numerical simulations where the value of ranged from to . Figures 1(a) and 1(b) show that the value of is , the spatial distribution of prey is a spiral stripy wave, and it is unstable. The instability is uncertain, so we expect that instability can be avoided by increasing the value of . In Figures 1(c) and 1(d), the value of is and we can observe a periodic spiral wave, because the phytoplankton gather in the form of a continuous spiral stripe to avoid the predator. Finally, in Figures 1(e) and 1(f) the value of increases to and the spatial distribution of prey is a periodic spiral wave, where the width of the stripe increases. Thus, the living space of the prey grows with increasing values of . To better understand the effects of prey refuges, we present the spatiotemporal series and phase diagram in Figure 2, using the same parameter values shown in Figure 1. This series of snapshots also shows that we can increase the number of prey refuges to avoid exacerbating algal growth.

Figure 3 shows the spatial patterns of the phytoplankton when the value of increases to 1.2. Based on the above theoretical analysis, we can see that phytoplankton and zooplankton will coexist. With increasing iterations, we can obtain the stable wave pattern shown in the interior of Figure 3(c), which is limited to a rectangular area. Thus, the phytoplankton is limited to this area and its value is almost constant in this area. Figure 3(d) also further that the state of the phytoplankton is changed from its original state to a final stable state.

Based on the analysis above, we can see that model (3) has traveling wave solutions, which supports all of the results above. If the value of is 1.2, the traveling wave solution of model (3) and the corresponding spatial trend are as shown in Figure 4. It is important to note that the phytoplankton density in Figure 4 changes from to with spatial variation and time variation , before finally reaching a steady state. We also used different time traveling waves , , , at times , , , , respectively. Figure 4 shows that the trajectories of the phytoplankton population are different at different times, although all trajectories started at a steady state , before finally moving to a steady state . Furthermore, it is important to note that model (3) has a traveling wave solution connecting the steady state and the coexisting equilibrium points . These results further show that the coexisting equilibrium points becomes more stable as the value of increase.

#### 5. Discussion

In this paper, we considered a phytoplankton-zooplankton model given by two coupled reaction-diffusion equations with a Holling’s type II functional response and alternative prey and refuge effects with homogeneous Neumann boundary conditions. Mathematical analyses were used to investigate the stability of the equilibrium points and the critical conditions for traveling wave solution. The numerical analysis indicated that parameter has an important effect on the spatiotemporal dynamics of model (3). Earlier in the article, . When the carrying capacity and half-saturation constant are fixed, the value of is only determined by the refuge protecting the prey , which is changed readily by different biological and chemical factors. These result shows that the prey refuge has an important role in the spatiotemporal dynamics of model (1), where it can make the system more stable thereby stabilizing the density of phytoplankton, which helps to avoid the exacerbation of algal growth. These results may help us to understand the effects of the undoubted susceptibility to diffusion in real ecosystems.

#### Acknowledgments

This work was supported by the National Key Basic Research Program of China (973 Program, Grant no. 2012CB426510), by the National Natural Science Foundation of China (Grant no. 31170338 and 30970305), by the Key Program of Zhejiang Provincial Natural Science Foundation of China (Grant no. LZ12C03001) and by the Program for Wenzhou Science & Technology Innovative Research Team of China (Grant no. C20120007).

#### References

- K. Das and N. H. Gazi, “Random excitations in modelling of algal blooms in estuarine systems,”
*Ecological Modelling*, vol. 222, no. 14, pp. 2495–2501, 2011. View at Publisher · View at Google Scholar · View at Scopus - E. Beltrami, “Unusual algal blooms as excitable systems: the case of “brown-tides,”
*Environmental Modeling and Assessment*, vol. 1, pp. 19–24, 1996. View at Google Scholar - C. Dai and M. Zhao, “Mathematical and dynamic analysis of a prey-predator model in the presence of alternative prey with impulsive state feedback control,”
*Discrete Dynamics in Nature and Society*, vol. 2012, Article ID 724014, 19 pages, 2012. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - R. May,
*Stability and Complexity in Model Ecosystems with a New Introduction by the Author*, Princeton University Press, Princeton, NJ, USA, 1973. - H. I. Freedman,
*Deterministic Mathematical Models in Population Ecology*, vol. 57, Marcel Dekker, New York, NY, USA, 1980. View at MathSciNet - M. Gatto, “Some remarks of models of piankton densities in lake,”
*The American Naturalist*, vol. 137, no. 3, pp. 264–267, 1991. View at Google Scholar - M. Haque, “A detailed study of the Beddington-DeAngelis predator-prey model,”
*Mathematical Biosciences*, vol. 234, no. 1, pp. 1–16, 2011. View at Publisher · View at Google Scholar · View at MathSciNet - H. R. Akcakaya, R. Arditi, and L. R. Ginzburg, “Ratio-dependent predation: an abstraction that works,”
*Ecology*, vol. 76, no. 3, pp. 995–1004, 1995. View at Google Scholar · View at Scopus - C. Cosner, D. L. Deangelis, J. S. Ault, and D. B. Olson, “Effects of spatial grouping on the functional response of predators,”
*Theoretical Population Biology*, vol. 56, no. 1, pp. 65–75, 1999. View at Publisher · View at Google Scholar · View at Scopus - A. P. Gutierrez, “Physiological basis of ratio-dependent predator-prey theory: the metabolic pool model as a paradigm,”
*Ecology*, vol. 73, no. 5, pp. 1552–1563, 1992. View at Google Scholar · View at Scopus - P. A. Abrams, “The fallacies of
*″*ratio-dependent*″*predation,”*Ecology*, vol. 75, no. 6, pp. 1842–1850, 1994. View at Google Scholar · View at Scopus - S. B. Hsu, T. W. Hwang, and Y. Kuang, “Rich dynamics of a ratio-dependent one-prey two-predators model,”
*Journal of Mathematical Biology*, vol. 43, no. 5, pp. 377–396, 2001. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - Y. Kuang, “Rich dynamics of Gause-type ratio-dependent predator-prey system,”
*The Fields Institute Communications*, vol. 21, pp. 325–337, 1999. View at Google Scholar - C. Duque and M. Lizana, “On the dynamics of a predator-prey model with nonconstant death rate and diffusion,”
*Nonlinear Analysis*, vol. 12, no. 4, pp. 2198–2210, 2011. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - B. I. Camara, “Waves analysis and spatiotemporal pattern formation of an ecosystem model,”
*Nonlinear Analysis*, vol. 12, no. 5, pp. 2511–2528, 2011. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - X. Guan, W. Wang, and Y. Cai, “Spatiotemporal dynamics of a Leslie-Gower predator-prey model incorporating a prey refuge,”
*Nonlinear Analysis*, vol. 12, no. 4, pp. 2385–2395, 2011. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - T. K. Kar and S. K. Chattopadhyay, “A focus on long-run sustainability of a harvested prey predator system in the presence of alternative prey,”
*Comptes Rendus*, vol. 333, no. 11-12, pp. 841–849, 2010. View at Publisher · View at Google Scholar · View at Scopus - P. A. Abrams and H. Matsuda, “Positive indirect effects between prey species that share predators,”
*Ecology*, vol. 77, no. 2, pp. 610–616, 1996. View at Google Scholar · View at Scopus - W. T. Li and S. L. Wu, “Traveling waves in a diffusive predator-prey model with Holling type-III functional response,”
*Chaos, Solitons & Fractals*, vol. 37, no. 2, pp. 476–486, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - S. R. Dunbar, “Traveling waves in diffusive predator-prey equations: periodic orbits and point-to-periodic heteroclinic orbits,”
*SIAM Journal on Applied Mathematics*, vol. 46, no. 6, pp. 1057–1078, 1986. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - R. A. Gardner, “Existence of travelling wave solutions of predator-prey systems via the connection index,”
*SIAM Journal on Applied Mathematics*, vol. 44, no. 1, pp. 56–79, 1984. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - K. Mischaikow and J. F. Reineck, “Travelling waves in predator-prey systems,”
*SIAM Journal on Mathematical Analysis*, vol. 24, no. 5, pp. 1179–1214, 1993. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - M. R. Owen and M. A. Lewis, “How predation can slow, stop or reverse a prey invasion,”
*Bulletin of Mathematical Biology*, vol. 63, no. 4, pp. 655–684, 2001. View at Publisher · View at Google Scholar · View at Scopus - A. I. Volpert and V. A. Volpert,
*Traveling Wave Solutions of Parabolic Systems*, vol. 140 of*Translations of Mathematical Monographs*, American Mathematical Society, Providence, RI, USA, 1994. View at MathSciNet - F. F. Zhang, G. Huo, Q. X. Liu, G. Q. Sun, and Z. Jin, “Existence of travelling waves in nonlinear SI epidemic models,”
*Journal of Biological Systems*, vol. 17, no. 4, pp. 643–657, 2009. View at Publisher · View at Google Scholar · View at MathSciNet - J. Yang and M. Zhao, “A mathematical model for the dynamics of a fish algae consumption model with impulsive control strategy,”
*Journal of Applied Mathematics*, vol. 2012, Article ID 575047, 17 pages, 2012. View at Publisher · View at Google Scholar - H. Yu and M. Zhao, “Seasonally perturbed prey-predator ecological system with the Beddington-DeAngelis functional response,”
*Discrete Dynamics in Nature and Society*, vol. 2012, Article ID 150359, 12 pages, 2012. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - M. Zhao, X. Wang, H. Yu, and J. Zhu, “Dynamics of an ecological model with impulsive control strategy and distributed time delay,”
*Mathematics and Computers in Simulation*, vol. 82, no. 8, pp. 1432–1444, 2012. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - P. P. Liu, “An analysis of a predator-prey model with both diffusion and migration,”
*Mathematical and Computer Modelling*, vol. 51, no. 9-10, pp. 1064–1070, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - N. H. Gazi and K. Das, “Structural stability analysis of an algal bloom mathematical model in tropic interaction,”
*Nonlinear Analysis*, vol. 11, no. 4, pp. 2191–2206, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - C. S. Holling, “The components of predation as revealed by a study of small mammal predation on the European pine sawfly,”
*The Canadian Entomologist*, vol. 91, no. 5, pp. 293–320, 1989. View at Google Scholar - W. Wang, Y. Lin, F. Yang, L. Zhang, and Y. Tan, “Numerical study of pattern formation in an extended Gray-Scott model,”
*Communications in Nonlinear Science and Numerical Simulation*, vol. 16, no. 4, pp. 2016–2026, 2011. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - R. S. Cantrell and C. Cosner,
*Spatial Ecology via Reaction-Diffusion Equations*, John Wiley & Sons, 2003. View at Publisher · View at Google Scholar · View at MathSciNet - W. Wang, Y. Cai, M. Wu, K. Wang, and Z. Li, “Complex dynamics of a reaction-diffusion epidemic model,”
*Nonlinear Analysis*, vol. 13, no. 5, pp. 2240–2258, 2012. View at Publisher · View at Google Scholar · View at MathSciNet - E. E. Holmes, M. A. Lewis, J. E. Banks, and R. R. Veit, “Partial differential equations in ecology: spatial interactions and population dynamics,”
*Ecology*, vol. 75, no. 1, pp. 17–29, 1994. View at Google Scholar · View at Scopus - H. K. Khalil,
*Nonlinear Systems*, Macmillan Publishing Company, New York, NY, USA, 1992. View at MathSciNet - P. Hartman,
*Ordinary Differential Equations*, John Wiley & Sons, New York, NY, USA, 1973. View at MathSciNet