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

Authors: Weiwei Zhang; Min Zhao
Journal: Journal of Applied Mathematics
 (2013)

---

## 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.


---

## Body

## 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 Section2, 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:(1)dXdτ=r1X(1-Xk1)-a1(1-m1)Xb1+(1-m1)XY,dYdτ=a1e1(1-m1)Xb1+(1-m1)XY+s1Y(1-Xk1)-d1Y,
where X and Y are the densities of phytoplankton and zooplankton at time τ and position (x,y), respectively. The parameters of the model (1) can be expressed as follows: r1 represents the inherent growth rate of the phytoplankton without any environment limitations; k1 is the carrying capacity of the phytoplankton in the presence of predators and harvesting; a1 represents the consumption rate; b1 is the half-saturation constant that determines how quickly the maximum is attained when the phytoplankton density increases; e1 is the conversion factor denoting the number of new herbivores produced for each capture; s1 is the growth rate of micro-grazers due to the alternative prey; d1 denotes the food-independent zooplankton mortality rate where d1>s1; m1 is the refuge protect the prey, where m1∈[0,1) is constant [29–32].For simplicity, using the scalingX=k1N, Y=(r1k1/a1)P, and τ=t/r1, we have
(2)dNdt=N(1-N)-Nb+NP,dPdt=eNb+NP+sP(1-N)-dP,
where b=b1/(1-m1)k1,  e=a1e1/r1,  s=s1/r1, and d=d1/r1. 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:
(3)∂N∂t=N(1-N)-Nb+NP+D1∇2N,∂P∂t=eNb+NP+sP(1-N)-dP+D2∇2P,
where ∇2=∂2/∂x2+∂2/∂y2 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 D1 and D2 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 R+2. Thus, model (3) is analyzed using the following initial conditions:
(4)N0=N(x,y·0)≥0,P0=P(x,y·0)≥0,(x,y)∈Ω=[0,L]×[0,L]
and the homogeneous Neumann boundary condition:
(5)∂N∂n=∂P∂n=0,(x,y)∈∂Ω,
which indicates that there are no population fluxes through the boundary and that no external input is imposed from the outside [19]. n 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(N,P) be any solution of model (3). Thus, we can find
(6)limt→∞supmaxx∈Ω-N(x,t)≤1,limt→∞supmaxx∈Ω-P(x,t)≤b+1.Proof.
Based on the first equation in model (3), we derive
(7)∂N∂t-D1∇2N≤N(1-N).
Based on the comparison principle, letz(t) be the solution of the equation
(8)z′(t)=z(t)(1-z(t)),z(0)=maxx∈Ω-N(x,0).
Thus, we can obtainN(x,t)≤z(t). For any δ>0, there is a Tδ>0, such that z(t)<1+δ for any t≥Tδ. Therefore, N(x,t) is defined for all t>0 and limt→∞sup maxx∈Ω-N(x,t)≤1.
Next, we will discuss the boundedness ofP(x,t) based on the second equation of model (3). We let β=d-s>0 and we find that
(9)∂P∂t-D2∇2P≤eNPb+N-βP.
Assuming that there is aT=T(M0,φ)>0 such that
(10)P(x,t)≥M0,∀t≥T,x∈Ω-,whereM0=b+1,we can determine P/(b+N)≥1, for all t≥T*, x∈Ω-, where T*=max{Tδ,T}.
Based on the formulae above and the first equation in model (3), we can obtain
(11)∂N∂t-D1∇2N=-N2+N(1-pb+N)≤-N2,∀t≥T*
so we find that
(12)0≤N(x,t)≤z(t)⟶0,t⟶∞,x∈Ω-,
where  z(t) is the solution of the equation z′(t)=-z2(t) and z(T*)=maxx∈Ω-N(x,T*).
Then,N(x,t)→0 as t→∞,   x∈Ω-. Thus, there is a t0≥T*such that
(13)eNb+N≤β2,∀t≥t0,x∈Ω-.
Substituting inequality (13) into the second equation of model (3), we can get
(14)∂P∂t-D2∇2P≤β2P-βP=-β2P,∀t≥t0,x∈Ω-.
Based on the comparison principle, we can find thatP(x,t)→0 as t→∞, x∈Ω-, although this is clearly contradictory.
LetΨ(x,t)=P(x,t)-M0, t>0, x∈Ω-. If there is a t0>0, such that Ψ(x,t)≠0, for all t>t0, x∈Ω-, then it is considered to be bounded and obviously we can find that 0<P(x,t)<M0 for for all t>t0, x∈Ω-.
IfΨ(x,t)=0 for a sequence of tn tending to +∞ as n→∞, we can say that it is unbounded, that is, the t-axis is divided by the zeroes into a sequence of intervals Ln, n>0. Given this assumption, Ψ(x,t)≥0, for all t∈Ln, n=1,3,5,…. For the intervals L2n-1, arguing the inequality (10) is not possible, so we find that P(x,t)≥M0∃A∈N+ if n>A, for all t∈L2n-1,   x∈Ω-,  ∂P/∂t-D2∇2P≤-(β/2)P.
Based on the comparison principle,P(x,t)≤z(t)  for t∈L2n-1, where z(t) is the solution of the equation z′(t)=-(β/2)z(t) given suitable initial conditions. Thus,
(15)P(x,t)≤z(t)≤z(t2n-1)=M0,n>A,∀t∈L2n-1,x∈Ω-.
Therefore, we can findlimt→∞sup maxx∈Ω-P(x,t)≤M0=b+1.
Among the conditions that satisfy Theorem1, we mainly address the steady states and their stabilities for (3), as follows.
(i)
E
0
=
(
0,0
), which corresponds to the extinction of the phytoplankton and zooplankton, which is apparently an unstable saddle point.
(ii)
E
1
=
(
1,0
), which corresponds to the extinction of the zooplankton.
(iii)
Their coexistence will be observed in different statesE2=(n2,p2) and E3=(n3,p3), where
(16)n2=s+e-d-bs-(s+e-d-bs)2+4bs(s-d)2s,n3=s+e-d-bs+(s+e-d-bs)2+4bs(s-d)2sand pi=(ni+b)(1-ni)  (i=2.3). Obviously, if s+e-d-bs<0,   n2 and n3 are negative. In this study, we assume that s+e-bs>d>s is satisfied. E2=(n2,p2) and E3=(n3,p3) have similar properties, thus we focus only on the positive equilibrium point E2=(n2,p2).
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 pointsEi(i=0,1,2). Next, we determine the stability of the equilibrium points Ei(i=0,1,2) for model (3). Substituting
(17)(N(x,y·t)P(x,y·t))=EiT+(U1(x,y·t)U2(x,y·t))
into model (3), where |Ui|≪1(i=1,2), and selecting all terms that are linear about U,
(18)∂U∂t=D∇2U+J(Ei)U,
where
(19)D=(D100D2),J(Ei)=(J11J12J21J22)|Ei,
where
(20)J11=1-2N-bP(b+N)2,J12=-Nb+N,J21=beP(b+N)2-sP,J22=eNb+N+s-d-sN.
Substituting (19) into (18), we obtain:
(21)∂U1∂t=D1∇2U1+J11U1+J12U2,∂U2∂t=D2∇2U2+J21U1+J22U2.


### 3.1. Stability Analysis

Theorem 2.
(i)E0=(0,0) is unstable;
(ii)  ifd>e/(b+1), then  E1=(1,0) for model (3) is asymptotically stable;
(iii)E2=(n2,p2)  is locally asymptotically stable when s+e-bs>d>s  and (a)
ifb<1, then (1-b)/2<n2<min(1,be/s-b),(b)
ifb≥1, then 0≤n2<min(1,be/s-b).Proof.
The stability of the equilibrium state is determined based on the nature of the eigenvalues in the Jacobian matrixJ(Ei). Substituting E0 into (19), we can determine the eigenvalues of the Jacobian matrix J(E0), 1 and s-d. Apparently s-d<0, so E0 is unstable. The stability of E1 can be established in a similar manner as E2, so we only discuss the stability of E2.
We can useλi to determine the eigenvalues of the operator -Δ on Ω, which satisfy the homogeneous Neumann boundary condition and ⋯>λn>⋯>λ4>λ3>λ2>λ1=0. Γi is the eigenvector relative to λi  (i=1,2,3,…), so
(22)λiΓi+ΔΓi=0inΩ,∂Γi∂n=0on∂Ω.
By substitutingE2=(n2,p2) into (18) and (19), and expanding the solution U of (18),
(23)U=∑i=1∞τi(t)Γi,τi(t)∈R2.
Substituting (23) into (18) and balancing the coefficients about each Γi, we can obtain dτi(t)/dt=Tiτi  where the matrix Ti satisfies Ti=J(E2)-λiD.
Thus,E2=(n2,p2) is stable if and only if each τi(t) decays to zero, that is, there are two eigenvalues with negative real parts for each matrix Ti. Based on simple algebraic computations, the characteristic polynomial of Ti is determined as follow,
(24)γ2-tr(Ti)γ+det(Ti)=0,tr(Ti)=J11+J22-λi(D1+D2),det(Ti)=λi2D1D2-λi(J11D2+J22D1)+J11J22-J12J21.
Iftr(Ti)<0 and det(Ti)>0 is satisfied, the above conditions are met. Therefore, with a random λi≥0, tr(Ti)<0 and det(Ti)>0 hold if J11<0 and J12J21<0.
Thus, we can obtain:(a)
ifb<1, then (1-b)/2<n2<min(1,be/s-b),(b)
ifb≥1, then 0≤n2<min(1,be/s-b). This completes the proof.Theorem 3.
If1<b<(e-2s+e2-4es)/2s, e≥4s and  (B-s)b>s hold, the steady state E2=(n2,p2) of the model (3) is globally asymptotically stable, where B=(en2)/(b+n2)+s(1-n2).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,
(25)V(N,P)=∫n2Ne+(d(1-u)-en2b+n2-d(1-n2))×××(1+bu)du+∫p2Pu-p2udu,
assuming  0<b<(e-2s+e2-4es)/2s and e≥4s are satisfied.
(26)W(t)=∬ΩV(N,P)dΩV(N,P) is positive for all (N,P) in the positive quadrant, except for E2=(n2,p2) where V(n2,p2)=0. If dW/dt<0 is proven, then W(t) is the Lyapunov function. Calculating the rate of change of W(t) for the solution of model (3), we find that,
(27)dWdt=∬Ω((e+(s(1-N)-en2b+n2-s(1-n2))(1+bN))×NN××(N(1-N)-NPb+N+D1∇2N)+P-p2P×NN××(eNPb+N+sP(1-N)-dP+D2∇2P)bN)dΩ=∬Ω((e+(s(1-N)-en2b+n2-s(1-n2))(1+bN))×NN××(N(1-N)-NPb+N)+P-p2P×NN××(eNPb+N+sP(1-N)-dP)×NN×+(e+(s(1-N)-en2b+n2-s(1-n2))(1+bN))×NN×+P-p2PD2∇2Pen2b+n2)dΩ=H1+H2.
Based on simple algebraic computations,H1 becomes
(28)H1=∬Ω(f(N)-f(n2))(g(N)-g(n2))dΩ,
where f(x)=ex/(b+x)+s(1-x) and g(x)=(1-x)(b+x).
Obviously,f(x) increases strictly monotonically, so H1 is negative for all (N,P) in the positive quadrant if and only if g(x) decreases strictly monotonically, and we find that 1<b<(e-2s+e2-4es)/2s and e≥4s holds.
In order to determineH2, we can use the Green formula. If it satisfies the homogeneous Neumann boundary conditions, we find
(29)H2=-∬Ω(D1|∇N|2ddN×(e+(s(1-N)-en2b+n2-s(1-n2))(1+bN))+D2|∇P|2ddP(P-p2P))dΩ=-∬Ω(D1|∇N|2((B-s)bN2-s)+D2|∇P|2p2P2)dΩ,
where B=en2/(b+n2)+s(1-n2). Based on Theorem 1, we know 0<N<1. If (B-s)b>s is found, H2<0. If 1<b<(e-2s+e2-4es)/2s, e≥4s and (B-s)b>s hold, so dW/dt is negative. Therefore, E2=(n2,p2) is globally asymptotically stable.


### 3.2. Traveling Wave Solutions

In the previous section, we discussed the stability of the equilibrium pointsEi  (i=0,1,2) for model (3). It can be got that E0=(0,0) is unstable, E1=(1,0) is unstable for model (3), and the steady state E2=(n2,p2) 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 E1=(1,0) and E2=(n2,p2) or E0=(0,0) and E2=(n2,p2).To find the heteroclinic orbit betweenE1=(1,0) and E2=(n2,p2), we need to show that model (3) has a solution with the special form N(x,t)=N(x+ct), P(x,t)=P(x+ct), where the wave speed parameter c is positive. After substituting N(x,t)=N(x+ct), P(x,t)=P(x+ct), and s=x+ct into (3), we get the following wave equations
(30)cN′=N(1-N)-Nb+NP+D1N′′,cP′=eNb+NP+sP(1-N)-dP+D2P′′.IfN(s) and P(s) are nonnegative and they satisfy the following boundary conditions
(31)N(-∞)=1,N(+∞)=n2,P(-∞)=0,P(+∞)=p2,
which shows there is a traveling wave solution between E1=(1,0) and E2=(n2,p2).Substitutingu=N′ and w=P′ into (30),
(32)N′=u,u′=cD1u-N(1-N)D1+NPD1(b+N),P′=w,w′=cD2w-eNPD2(b+N)-sP(1-N)D2+dPD2.After computing the Jacobian matrix for model (32) at (1,0,0,0), we get the following matrix:
(33)J(1,0,0,0)=[01001D1cD11D1(b+1)0000100dD2-1D2(b+1)cD2].Therefore, the characteristic equation of(1,0,0,0) is given by |ηE-J(1,0,0,0)|=0, and ηi are its eigenvalues, as follows
(34)η1=c+c2+4D1D1,η2=c-c2+4D1D1,η3=c+c2-4(1/(b+1)-d)D2,η4=c-c2-4(1/(b+1)-d)D2.If0<c<c*=2D2(e/(b+1)-d), then η3 and η4 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 (1,0,0,0). If s→-∞, the trajectory approaching (1,0,0,0) must have P(s)<0 for some s. However, this is contradictory if 0<P<b+1.This discussion shows that because0<c<c*=2D2(e/(b+1)-d), traveling wave solutions cannot be found for model (32), whereas if c>2D2(e/(b+1)-d, model (32) has non-negative solutions that satisfy the boundary conditions. In some suitable conditions, therefore, there is a traveling wave solution from (1,0,0,0) to (n2,0,p2,0). Next, we discuss the Jacobian matrix of model (32) at (n2,0,p2,0),
(35)J(n2,0,p2,0)=[0100a21cD1a2300001a4100cD2],
where a21=-(1-2n2)/D1+b(1-n2)/D1(b+n2), a23=n2/D1(b+n2), and a41=-eb(1-n2)/D2(b+n2)+s(b+n2)(1-n2)/D2. Thus, the characteristic equation of (n2,0,p2,0) is given by |ηE-J(n2,0,p2,0)|=0 and η is its eigenvalues, as follows: T(η)=η4-(c/D1+c/D2)η3+(c2/D1D2-a21)η2+(ca21/D2)η-a23a41.Based on the Routh-Hurwitz stability criterion, whena41<0 and c>D1D2a21, the characteristic equation of (n2,0,p2,0) 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 (n2,0,p2,0) and a trajectory from (1,0,0,0) to (n2,0,p2,0) where s→∞. Thus, if a41<0 and c>max(2D2(e/(b+1)-d),D1D2a21) hold, there is a traveling wave solution.


## 3.1. Stability Analysis

Theorem 2.
(i)E0=(0,0) is unstable;
(ii)  ifd>e/(b+1), then  E1=(1,0) for model (3) is asymptotically stable;
(iii)E2=(n2,p2)  is locally asymptotically stable when s+e-bs>d>s  and (a)
ifb<1, then (1-b)/2<n2<min(1,be/s-b),(b)
ifb≥1, then 0≤n2<min(1,be/s-b).Proof.
The stability of the equilibrium state is determined based on the nature of the eigenvalues in the Jacobian matrixJ(Ei). Substituting E0 into (19), we can determine the eigenvalues of the Jacobian matrix J(E0), 1 and s-d. Apparently s-d<0, so E0 is unstable. The stability of E1 can be established in a similar manner as E2, so we only discuss the stability of E2.
We can useλi to determine the eigenvalues of the operator -Δ on Ω, which satisfy the homogeneous Neumann boundary condition and ⋯>λn>⋯>λ4>λ3>λ2>λ1=0. Γi is the eigenvector relative to λi  (i=1,2,3,…), so
(22)λiΓi+ΔΓi=0inΩ,∂Γi∂n=0on∂Ω.
By substitutingE2=(n2,p2) into (18) and (19), and expanding the solution U of (18),
(23)U=∑i=1∞τi(t)Γi,τi(t)∈R2.
Substituting (23) into (18) and balancing the coefficients about each Γi, we can obtain dτi(t)/dt=Tiτi  where the matrix Ti satisfies Ti=J(E2)-λiD.
Thus,E2=(n2,p2) is stable if and only if each τi(t) decays to zero, that is, there are two eigenvalues with negative real parts for each matrix Ti. Based on simple algebraic computations, the characteristic polynomial of Ti is determined as follow,
(24)γ2-tr(Ti)γ+det(Ti)=0,tr(Ti)=J11+J22-λi(D1+D2),det(Ti)=λi2D1D2-λi(J11D2+J22D1)+J11J22-J12J21.
Iftr(Ti)<0 and det(Ti)>0 is satisfied, the above conditions are met. Therefore, with a random λi≥0, tr(Ti)<0 and det(Ti)>0 hold if J11<0 and J12J21<0.
Thus, we can obtain:(a)
ifb<1, then (1-b)/2<n2<min(1,be/s-b),(b)
ifb≥1, then 0≤n2<min(1,be/s-b). This completes the proof.Theorem 3.
If1<b<(e-2s+e2-4es)/2s, e≥4s and  (B-s)b>s hold, the steady state E2=(n2,p2) of the model (3) is globally asymptotically stable, where B=(en2)/(b+n2)+s(1-n2).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,
(25)V(N,P)=∫n2Ne+(d(1-u)-en2b+n2-d(1-n2))×××(1+bu)du+∫p2Pu-p2udu,
assuming  0<b<(e-2s+e2-4es)/2s and e≥4s are satisfied.
(26)W(t)=∬ΩV(N,P)dΩV(N,P) is positive for all (N,P) in the positive quadrant, except for E2=(n2,p2) where V(n2,p2)=0. If dW/dt<0 is proven, then W(t) is the Lyapunov function. Calculating the rate of change of W(t) for the solution of model (3), we find that,
(27)dWdt=∬Ω((e+(s(1-N)-en2b+n2-s(1-n2))(1+bN))×NN××(N(1-N)-NPb+N+D1∇2N)+P-p2P×NN××(eNPb+N+sP(1-N)-dP+D2∇2P)bN)dΩ=∬Ω((e+(s(1-N)-en2b+n2-s(1-n2))(1+bN))×NN××(N(1-N)-NPb+N)+P-p2P×NN××(eNPb+N+sP(1-N)-dP)×NN×+(e+(s(1-N)-en2b+n2-s(1-n2))(1+bN))×NN×+P-p2PD2∇2Pen2b+n2)dΩ=H1+H2.
Based on simple algebraic computations,H1 becomes
(28)H1=∬Ω(f(N)-f(n2))(g(N)-g(n2))dΩ,
where f(x)=ex/(b+x)+s(1-x) and g(x)=(1-x)(b+x).
Obviously,f(x) increases strictly monotonically, so H1 is negative for all (N,P) in the positive quadrant if and only if g(x) decreases strictly monotonically, and we find that 1<b<(e-2s+e2-4es)/2s and e≥4s holds.
In order to determineH2, we can use the Green formula. If it satisfies the homogeneous Neumann boundary conditions, we find
(29)H2=-∬Ω(D1|∇N|2ddN×(e+(s(1-N)-en2b+n2-s(1-n2))(1+bN))+D2|∇P|2ddP(P-p2P))dΩ=-∬Ω(D1|∇N|2((B-s)bN2-s)+D2|∇P|2p2P2)dΩ,
where B=en2/(b+n2)+s(1-n2). Based on Theorem 1, we know 0<N<1. If (B-s)b>s is found, H2<0. If 1<b<(e-2s+e2-4es)/2s, e≥4s and (B-s)b>s hold, so dW/dt is negative. Therefore, E2=(n2,p2) is globally asymptotically stable.


## 3.2. Traveling Wave Solutions

In the previous section, we discussed the stability of the equilibrium pointsEi  (i=0,1,2) for model (3). It can be got that E0=(0,0) is unstable, E1=(1,0) is unstable for model (3), and the steady state E2=(n2,p2) 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 E1=(1,0) and E2=(n2,p2) or E0=(0,0) and E2=(n2,p2).To find the heteroclinic orbit betweenE1=(1,0) and E2=(n2,p2), we need to show that model (3) has a solution with the special form N(x,t)=N(x+ct), P(x,t)=P(x+ct), where the wave speed parameter c is positive. After substituting N(x,t)=N(x+ct), P(x,t)=P(x+ct), and s=x+ct into (3), we get the following wave equations
(30)cN′=N(1-N)-Nb+NP+D1N′′,cP′=eNb+NP+sP(1-N)-dP+D2P′′.IfN(s) and P(s) are nonnegative and they satisfy the following boundary conditions
(31)N(-∞)=1,N(+∞)=n2,P(-∞)=0,P(+∞)=p2,
which shows there is a traveling wave solution between E1=(1,0) and E2=(n2,p2).Substitutingu=N′ and w=P′ into (30),
(32)N′=u,u′=cD1u-N(1-N)D1+NPD1(b+N),P′=w,w′=cD2w-eNPD2(b+N)-sP(1-N)D2+dPD2.After computing the Jacobian matrix for model (32) at (1,0,0,0), we get the following matrix:
(33)J(1,0,0,0)=[01001D1cD11D1(b+1)0000100dD2-1D2(b+1)cD2].Therefore, the characteristic equation of(1,0,0,0) is given by |ηE-J(1,0,0,0)|=0, and ηi are its eigenvalues, as follows
(34)η1=c+c2+4D1D1,η2=c-c2+4D1D1,η3=c+c2-4(1/(b+1)-d)D2,η4=c-c2-4(1/(b+1)-d)D2.If0<c<c*=2D2(e/(b+1)-d), then η3 and η4 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 (1,0,0,0). If s→-∞, the trajectory approaching (1,0,0,0) must have P(s)<0 for some s. However, this is contradictory if 0<P<b+1.This discussion shows that because0<c<c*=2D2(e/(b+1)-d), traveling wave solutions cannot be found for model (32), whereas if c>2D2(e/(b+1)-d, model (32) has non-negative solutions that satisfy the boundary conditions. In some suitable conditions, therefore, there is a traveling wave solution from (1,0,0,0) to (n2,0,p2,0). Next, we discuss the Jacobian matrix of model (32) at (n2,0,p2,0),
(35)J(n2,0,p2,0)=[0100a21cD1a2300001a4100cD2],
where a21=-(1-2n2)/D1+b(1-n2)/D1(b+n2), a23=n2/D1(b+n2), and a41=-eb(1-n2)/D2(b+n2)+s(b+n2)(1-n2)/D2. Thus, the characteristic equation of (n2,0,p2,0) is given by |ηE-J(n2,0,p2,0)|=0 and η is its eigenvalues, as follows: T(η)=η4-(c/D1+c/D2)η3+(c2/D1D2-a21)η2+(ca21/D2)η-a23a41.Based on the Routh-Hurwitz stability criterion, whena41<0 and c>D1D2a21, the characteristic equation of (n2,0,p2,0) 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 (n2,0,p2,0) and a trajectory from (1,0,0,0) to (n2,0,p2,0) where s→∞. Thus, if a41<0 and c>max(2D2(e/(b+1)-d),D1D2a21) 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(Lx,Ly)=(500,800). 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 N, whereas the snapshot with red parts represents the high density value of phytoplankton N. To investigate the interplay among several critical parameters, we determine the pattern formation using model (3) where other parameters were fixed as  d=0.1,   s=0.01,  e=0.4,  D1=1, and D2=0.01.We only changed the value ofb 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 b ranged from 0.1 to 0.5. Figures 1(a) and 1(b) show that the value of b is 0.1, 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 b. In Figures 1(c) and 1(d), the value of b is 0.2 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 b increases to 0.5 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 b. 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.Spatiotemporal dynamics of model (3). In ((a), (b)) b=0.1, ((c), (d)) b=0.2, and ((e), (f)) b=0.5. The iteration numbers are 140 in ((a), (c), (e)) and 3000 in ((b), (d), (f)).
(a)
(b)
(c)
(d)
(e)
(f)Spatiotemporal dynamics of model (3). ((a), (c), (e)) phytoplankton (2-D), ((b), (d), (f)) phase diagram. This analysis used the same set of parameter values as Figure 1.
(a)
(b)
(c)
(d)
(e)
(f)Figure3 shows the spatial patterns of the phytoplankton when the value of b 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.Spatiotemporal dynamics of model (3) where b=1.2. The number of iterations are as follows (a) 1, (b) 200, and (c) 3000. (d) The phase diagram.
(a)
(b)
(c)
(d)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 b 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 E1=(1,0) to E2=(n2,p2) with spatial variation x and time variation t, before finally reaching a steady state. We also used different time traveling waves L1, L2, L3, L4 at times t=3500, t=4000, t=4500, t=5000, respectively. Figure 4 shows that the trajectories of the phytoplankton population are different at different times, although all trajectories started at a steady state E1=(1,0), before finally moving to a steady state E2=(n2,p2). Furthermore, it is important to note that model (3) has a traveling wave solution connecting the steady state E1=(1,0) and the coexisting equilibrium points E2=(n2,p2). These results further show that the coexisting equilibrium points E2=(n2,p2) becomes more stable as the value of b increase.Traveling wave solutions for phytoplankton and the change in phytoplankton density at specific times.
(a)
(b)


## 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 parameterb has an important effect on the spatiotemporal dynamics of model (3). Earlier in the article, m1=1-b1/bk1. When the carrying capacity k1 and half-saturation constant b1 are fixed, the value of b is only determined by the refuge protecting the prey m1, which is changed readily by different biological and chemical factors. These result shows that the prey refuge m1 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.



---

*Source: 608073-2013-05-09.xml*