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

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: 𝑑𝑋𝑑𝜏=𝑟₁𝑋𝑋1−𝑘₁−𝑎₁1−𝑚₁𝑋𝑏₁+1−𝑚₁𝑋𝑌,𝑑𝑌=𝑎𝑑𝜏¹𝑒₁1−𝑚₁𝑋𝑏₁+1−𝑚₁𝑋𝑌+𝑠₁𝑌𝑋1−𝑘₁−𝑑₁𝑌, (1) 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 𝑚₁∈[0,1) is constant [29–32].

For simplicity, using the scaling 𝑋=𝑘₁𝑁, 𝑌=(𝑟₁𝑘₁/𝑎₁)𝑃, and 𝜏=𝑡/𝑟₁, we have 𝑑𝑁()−𝑁𝑑𝑡=𝑁1−𝑁𝑏+𝑁𝑃,𝑑𝑃=𝑑𝑡𝑒𝑁()𝑏+𝑁𝑃+𝑠𝑃1−𝑁−𝑑𝑃, (2) where 𝑏=𝑏₁/(1−𝑚₁)𝑘₁,  𝑒=𝑎₁𝑒₁/𝑟₁,  𝑠=𝑠₁/𝑟₁, 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: 𝜕𝑁()−𝑁𝜕𝑡=𝑁1−𝑁𝑏+𝑁𝑃+𝐷₁∇²𝑁,𝜕𝑃=𝜕𝑡𝑒𝑁()𝑏+𝑁𝑃+𝑠𝑃1−𝑁−𝑑𝑃+𝐷₂∇²𝑃, (3) 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: 𝑁₀𝑃=𝑁𝑥,𝑦⋅0≥0,₀[]×[]=𝑃𝑥,𝑦⋅0≥0,𝑥,𝑦∈Ω=0,𝐿0,𝐿 (4) and the homogeneous Neumann boundary condition: 𝜕𝑁=𝜕𝑛𝜕𝑃𝜕𝑛=0,𝑥,𝑦∈𝜕Ω, (5) 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 lim_𝑡→∞supmax_𝑥∈Ω𝑁(𝑥,𝑡)≤1,lim_𝑡→∞supmax_𝑥∈Ω𝑃(𝑥,𝑡)≤𝑏+1. (6)

Proof. Based on the first equation in model (3), we derive 𝜕𝑁𝜕𝑡−𝐷₁∇²𝑁≤𝑁(1−𝑁). (7)
Based on the comparison principle, let 𝑧(𝑡) be the solution of the equation 𝑧^′(𝑡)(𝑡(𝑡,𝑧(0)=𝑧)(1−𝑧))=max_𝑥∈Ω𝑁(𝑥,0). (8)
Thus, we can obtain 𝑁(𝑥,𝑡)≤𝑧(𝑡). For any 𝛿>0, there is a 𝑇_𝛿>0, such that 𝑧(𝑡)<1+𝛿 for any 𝑡≥𝑇_𝛿. Therefore, 𝑁(𝑥,𝑡) is defined for all 𝑡>0 and lim_𝑡→∞supmax_𝑥∈Ω𝑁(𝑥,𝑡)≤1.
Next, we will discuss the boundedness of 𝑃(𝑥,𝑡) based on the second equation of model (3). We let 𝛽=𝑑−𝑠>0 and we find that 𝜕𝑃𝜕𝑡−𝐷₂∇²𝑃≤𝑒𝑁𝑃𝑏+𝑁−𝛽𝑃. (9)
Assuming that there is a 𝑇=𝑇(𝑀₀,𝜑)>0 such that 𝑃(𝑥,𝑡)≥𝑀₀,∀𝑡≥𝑇,𝑥∈Ω,where𝑀₀=𝑏+1, (10)we can determine 𝑃/(𝑏+𝑁)≥1, for all 𝑡≥𝑇^∗, 𝑥∈Ω, where 𝑇^∗=max{𝑇_𝛿,𝑇}.
Based on the formulae above and the first equation in model (3), we can obtain 𝜕𝑁𝜕𝑡−𝐷₁∇²𝑁=−𝑁²𝑝+𝑁1−𝑏+𝑁≤−𝑁²,∀𝑡≥𝑇^∗ (11) so we find that 0≤𝑁(𝑥,𝑡)≤𝑧(𝑡)⟶0,𝑡⟶∞,𝑥∈Ω, (12) where 𝑧(𝑡) is the solution of the equation 𝑧^′(𝑡)=−𝑧²(𝑡) and 𝑧(𝑇^∗)=max_𝑥∈Ω𝑁(𝑥,𝑇^∗).
Then, 𝑁(𝑥,𝑡)→0 as 𝑡→∞,   𝑥∈Ω. Thus, there is a 𝑡₀≥𝑇^∗such that 𝑒𝑁≤𝛽𝑏+𝑁2,∀𝑡≥𝑡₀,𝑥∈Ω. (13)
Substituting inequality (13) into the second equation of model (3), we can get 𝜕𝑃𝜕𝑡−𝐷₂∇²𝛽𝑃≤2𝛽𝑃−𝛽𝑃=−2𝑃,∀𝑡≥𝑡₀,𝑥∈Ω. (14)
Based on the comparison principle, we can find that 𝑃(𝑥,𝑡)→0 as 𝑡→∞, 𝑥∈Ω, although this is clearly contradictory.
Let Ψ(𝑥,𝑡)=𝑃(𝑥,𝑡)−𝑀₀, 𝑡>0, 𝑥∈Ω. If there is a 𝑡₀>0, such that Ψ(𝑥,𝑡)≠0, for all 𝑡>𝑡₀, 𝑥∈Ω, then it is considered to be bounded and obviously we can find that 0<𝑃(𝑥,𝑡)<𝑀₀ for for all 𝑡>𝑡₀, 𝑥∈Ω.
If Ψ(𝑥,𝑡)=0 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 𝐿_𝑛, 𝑛>0. Given this assumption, Ψ(𝑥,𝑡)≥0, for all 𝑡∈𝐿_𝑛, 𝑛=1,3,5,…. For the intervals 𝐿_2𝑛−1, arguing the inequality (10) is not possible, so we find that 𝑃(𝑥,𝑡)≥𝑀₀∃𝐴∈𝑁⁺ if 𝑛>𝐴, for all 𝑡∈𝐿_2𝑛−1,   𝑥∈Ω,  𝜕𝑃/𝜕𝑡−𝐷₂∇²𝑃≤−(𝛽/2)𝑃.
Based on the comparison principle, 𝑃(𝑥,𝑡)≤𝑧(𝑡) for 𝑡∈𝐿_2𝑛−1, where 𝑧(𝑡) is the solution of the equation 𝑧^′(𝑡)=−(𝛽/2)𝑧(𝑡) given suitable initial conditions. Thus, 𝑃(𝑥,𝑡)≤𝑧(𝑡)≤𝑧𝑡_2𝑛−1=𝑀₀,𝑛>𝐴,∀𝑡∈𝐿_2𝑛−1,𝑥∈Ω. (15)
Therefore, we can find lim_𝑡→∞supmax_𝑥∈Ω𝑃(𝑥,𝑡)≤𝑀₀=𝑏+1.
Among the conditions that satisfy Theorem 1, we mainly address the steady states and their stabilities for (3), as follows.
(i)𝐸₀=(0,0), which corresponds to the extinction of the phytoplankton and zooplankton, which is apparently an unstable saddle point. (ii)𝐸₁=(1,0), which corresponds to the extinction of the zooplankton.(iii)Their coexistence will be observed in different states 𝐸₂=(𝑛₂,𝑝₂) and 𝐸₃=(𝑛₃,𝑝₃), where 𝑛₂=𝑠+𝑒−𝑑−𝑏𝑠−(𝑠+𝑒−𝑑−𝑏𝑠)²+4𝑏𝑠(𝑠−𝑑),𝑛2𝑠₃=𝑠+𝑒−𝑑−𝑏𝑠+(𝑠+𝑒−𝑑−𝑏𝑠)²+4𝑏𝑠(𝑠−𝑑)2𝑠 (16)and 𝑝_𝑖=(𝑛_𝑖+𝑏)(1−𝑛_𝑖)  (𝑖=2.3). Obviously, if 𝑠+𝑒−𝑑−𝑏𝑠<0,   𝑛₂ 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 𝐸_𝑖(𝑖=0,1,2). Next, we determine the stability of the equilibrium points 𝐸_𝑖(𝑖=0,1,2) for model (3). Substituting 𝑁𝑃𝑥,𝑦⋅𝑡𝑥,𝑦⋅𝑡=𝐸^𝑇_𝑖+𝑈₁𝑈𝑥,𝑦⋅𝑡₂𝑥,𝑦⋅𝑡 (17) into model (3), where |𝑈_𝑖|≪1(𝑖=1,2), and selecting all terms that are linear about 𝑈, 𝜕𝑈𝜕𝑡=𝐷∇²𝑈+𝐽𝐸_𝑖𝑈, (18) where 𝐷𝐷=¹00𝐷₂,𝐽𝐸_𝑖=𝐽₁₁𝐽₁₂𝐽₂₁𝐽_22𝐸𝑖, (19) where 𝐽₁₁=1−2𝑁−𝑏𝑃(𝑏+𝑁)²,𝐽₁₂𝑁=−,𝐽𝑏+𝑁₂₁=𝑏𝑒𝑃(𝑏+𝑁)²𝐽−𝑠𝑃,₂₂=𝑒𝑁𝑏+𝑁+𝑠−𝑑−𝑠𝑁. (20) Substituting (19) into (18), we obtain: 𝜕𝑈₁𝜕𝑡=𝐷₁∇²𝑈₁+𝐽₁₁𝑈₁+𝐽₁₂𝑈₂,𝜕𝑈₂𝜕𝑡=𝐷₂∇²𝑈₂+𝐽₂₁𝑈₁+𝐽₂₂𝑈₂. (21)

3.1. Stability Analysis

Theorem 2. (i)  𝐸₀=(0,0) is unstable;
(ii)  if 𝑑>𝑒/(𝑏+1), then 𝐸₁=(1,0) for model (3) is asymptotically stable;
(iii)  𝐸₂=(𝑛₂,𝑝₂) is locally asymptotically stable when 𝑠+𝑒−𝑏𝑠>𝑑>𝑠 and (a)if 𝑏<1, then (1−𝑏)/2<𝑛₂<min(1,𝑏𝑒/𝑠−𝑏),(b)if 𝑏≥1, then 0≤𝑛₂<min(1,𝑏𝑒/𝑠−𝑏).

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 𝐽(𝐸₀), 1 and 𝑠−𝑑. Apparently 𝑠−𝑑<0, 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 ⋯>𝜆_𝑛>⋯>𝜆₄>𝜆₃>𝜆₂>𝜆₁=0. Γ_𝑖 is the eigenvector relative to 𝜆_𝑖  (𝑖=1,2,3,…), so 𝜆_𝑖Γ_𝑖+ΔΓ_𝑖=0inΩ,𝜕Γ_𝑖𝜕𝑛=0on𝜕Ω. (22)
By substituting 𝐸₂=(𝑛₂,𝑝₂) into (18) and (19), and expanding the solution 𝑈 of (18), 𝑈=^∞_𝑖=1𝜏_𝑖(𝑡)Γ_𝑖,𝜏_𝑖(𝑡)∈𝑅². (23)
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, 𝛾²−tr𝑇_𝑖𝛾+det𝑇_𝑖𝑇=0,tr_𝑖=𝐽₁₁+𝐽₂₂−𝜆_𝑖𝐷₁+𝐷₂,det𝑇_𝑖=𝜆²_𝑖𝐷₁𝐷₂−𝜆_𝑖𝐽₁₁𝐷₂+𝐽₂₂𝐷₁+𝐽₁₁𝐽₂₂−𝐽₁₂𝐽₂₁. (24)
If tr(𝑇_𝑖)<0 and det(𝑇_𝑖)>0 is satisfied, the above conditions are met. Therefore, with a random 𝜆_𝑖≥0, tr(𝑇_𝑖)<0 and det(𝑇_𝑖)>0 hold if 𝐽₁₁<0 and 𝐽₁₂𝐽₂₁<0.
Thus, we can obtain: (a)if 𝑏<1, then (1−𝑏)/2<𝑛₂<min(1,𝑏𝑒/𝑠−𝑏),(b)if 𝑏≥1, then 0≤𝑛₂<min(1,𝑏𝑒/𝑠−𝑏). This completes the proof.

Theorem 3. If 1<𝑏<(𝑒−2𝑠+√𝑒²−4𝑒𝑠)/2𝑠, 𝑒≥4𝑠 and (𝐵−𝑠)𝑏>𝑠 hold, the steady state 𝐸₂=(𝑛₂,𝑝₂) of the model (3) is globally asymptotically stable, where 𝐵=(𝑒𝑛₂)/(𝑏+𝑛₂)+𝑠(1−𝑛₂).

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, 𝑉(𝑁,𝑃)=^𝑁_𝑛2𝑒+𝑑(1−𝑢)−𝑒𝑛₂𝑏+𝑛₂−𝑑1−𝑛₂×𝑏1+𝑢𝑑𝑢+^𝑃_𝑝2𝑢−𝑝₂𝑢𝑑𝑢, (25) assuming 0<𝑏<(𝑒−2𝑠+√𝑒²−4𝑒𝑠)/2𝑠 and 𝑒≥4𝑠 are satisfied. 𝑊(𝑡)=_Ω𝑉(𝑁,𝑃)𝑑Ω (26)𝑉(𝑁,𝑃) is positive for all (𝑁,𝑃) in the positive quadrant, except for 𝐸₂=(𝑛₂,𝑝₂) where 𝑉(𝑛₂,𝑝₂)=0. If 𝑑𝑊/𝑑𝑡<0 is proven, then 𝑊(𝑡) is the Lyapunov function. Calculating the rate of change of 𝑊(𝑡) for the solution of model (3), we find that, 𝑑𝑊=𝑑𝑡_Ω𝑒+𝑠(1−𝑁)−𝑒𝑛₂𝑏+𝑛₂−𝑠1−𝑛₂𝑏1+𝑁×𝑁(1−𝑁)−𝑁𝑃𝑏+𝑁+𝐷₁∇²𝑁+𝑃−𝑝₂𝑃×𝑒𝑁𝑃()𝑏+𝑁+𝑠𝑃1−𝑁−𝑑𝑃+𝐷₂∇²𝑃=𝑑Ω_Ω𝑒+𝑠(1−𝑁)−𝑒𝑛₂𝑏+𝑛₂−𝑠1−𝑛₂𝑏1+𝑁×𝑁(1−𝑁)−𝑁𝑃+𝑏+𝑁𝑃−𝑝₂𝑃×𝑒𝑁𝑃()+𝑠()−𝑏+𝑁+𝑠𝑃1−𝑁−𝑑𝑃𝑒+1−𝑁𝑒𝑛₂𝑏+𝑛₂−𝑠1−𝑛₂𝑏1+𝑁+𝑃−𝑝₂𝑃𝐷₂∇²𝑃𝑑Ω=𝐻₁+𝐻₂. (27)
Based on simple algebraic computations, 𝐻₁ becomes 𝐻₁=_Ω𝑓(𝑁)−𝑓𝑛₂𝑔(𝑁)−𝑔𝑛₂𝑑Ω, (28) where 𝑓(𝑥)=𝑒𝑥/(𝑏+𝑥)+𝑠(1−𝑥) and 𝑔(𝑥)=(1−𝑥)(𝑏+𝑥).
Obviously, 𝑓(𝑥) increases strictly monotonically, so 𝐻₁ is negative for all (𝑁,𝑃) in the positive quadrant if and only if 𝑔(𝑥) decreases strictly monotonically, and we find that 1<𝑏<(𝑒−2𝑠+√𝑒²−4𝑒𝑠)/2𝑠 and 𝑒≥4𝑠 holds.
In order to determine 𝐻₂, we can use the Green formula. If it satisfies the homogeneous Neumann boundary conditions, we find 𝐻₂=−_Ω𝐷₁|∇𝑁|²𝑑×𝑠()−𝑑𝑁𝑒+1−𝑁𝑒𝑛₂𝑏+𝑛₂−𝑠1−𝑛₂𝑏1+𝑁+𝐷₂|∇𝑃|²𝑑𝑑𝑃𝑃−𝑝₂𝑃𝑑Ω=−_Ω𝐷₁|∇𝑁|²(𝐵−𝑠)𝑏𝑁²−𝑠+𝐷₂|∇𝑃|²𝑝₂𝑃²𝑑Ω, (29) where 𝐵=𝑒𝑛₂/(𝑏+𝑛₂)+𝑠(1−𝑛₂). Based on Theorem 1, we know 0<𝑁<1. If (𝐵−𝑠)𝑏>𝑠 is found, 𝐻₂<0. If 1<𝑏<(𝑒−2𝑠+√𝑒²−4𝑒𝑠)/2𝑠, 𝑒≥4𝑠 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 𝐸_𝑖  (𝑖=0,1,2) for model (3). It can be got that 𝐸₀=(0,0) is unstable, 𝐸₁=(1,0) 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 𝐸₁=(1,0) and 𝐸₂=(𝑛₂,𝑝₂) or 𝐸₀=(0,0) and 𝐸₂=(𝑛₂,𝑝₂).

To find the heteroclinic orbit between 𝐸₁=(1,0) 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 𝑐𝑁^′()−𝑁=𝑁1−𝑁𝑏+𝑁𝑃+𝐷₁𝑁^′′,𝑐𝑃^′=𝑒𝑁()𝑏+𝑁𝑃+𝑠𝑃1−𝑁−𝑑𝑃+𝐷₂𝑃^′′. (30)

If 𝑁(𝑠) and 𝑃(𝑠) are nonnegative and they satisfy the following boundary conditions 𝑁(−∞)=1,𝑁(+∞)=𝑛₂,𝑃(−∞)=0,𝑃(+∞)=𝑝₂, (31) which shows there is a traveling wave solution between 𝐸₁=(1,0) and 𝐸₂=(𝑛₂,𝑝₂).

Substituting 𝑢=𝑁^′ and 𝑤=𝑃^′ into (30), 𝑁^′𝑢=𝑢,_′=𝑐𝐷₁𝑁()𝑢−1−𝑁𝐷₁+𝑁𝑃𝐷₁(),𝑃𝑏+𝑁_′𝑤=𝑤,_′=𝑐𝐷₂𝑤−𝑒𝑁𝑃𝐷₂()−()𝑏+𝑁𝑠𝑃1−𝑁𝐷₂+𝑑𝑃𝐷₂. (32)

After computing the Jacobian matrix for model (32) at (1,0,0,0), we get the following matrix: 𝐽_(1,0,0,0)=10100𝐷₁𝑐𝐷₁1𝐷₁()0𝑑𝑏+1000100𝐷₂−1𝐷₂()𝑐𝑏+1𝐷₂. (33)

Therefore, the characteristic equation of (1,0,0,0) is given by |𝜂𝐸−𝐽_(1,0,0,0)|=0, and 𝜂_𝑖 are its eigenvalues, as follows 𝜂₁=𝑐+𝑐²+4𝐷₁𝐷₁,𝜂₂=𝑐−𝑐²+4𝐷₁𝐷₁,𝜂₃=𝑐+𝑐²−4(1/(𝑏+1)−𝑑)𝐷₂,𝜂₄=𝑐−𝑐²−4(1/(𝑏+1)−𝑑)𝐷₂. (34)

If 0<𝑐<𝑐^∗=2𝐷₂(𝑒/(𝑏+1)−𝑑), 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 (1,0,0,0). If 𝑠→−∞, the trajectory approaching (1,0,0,0) must have 𝑃(𝑠)<0 for some 𝑠. However, this is contradictory if 0<𝑃<𝑏+1.

This discussion shows that because 0<𝑐<𝑐^∗=2𝐷₂(𝑒/(𝑏+1)−𝑑), traveling wave solutions cannot be found for model (32), whereas if 𝑐>2𝐷₂(𝑒/(𝑏+1)−𝑑, 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 (𝑛₂,0,𝑝₂,0). Next, we discuss the Jacobian matrix of model (32) at (𝑛₂,0,𝑝₂,0), 𝐽_{(𝑛2,0,𝑝2,0)}=𝑎0100₂₁𝑐𝐷₁𝑎₂₃0𝑎0001₄₁𝑐00𝐷₂, (35) where 𝑎₂₁=−(1−2𝑛₂)/𝐷₁+𝑏(1−𝑛₂)/𝐷₁(𝑏+𝑛₂), 𝑎₂₃=𝑛₂/𝐷₁(𝑏+𝑛₂), and 𝑎₄₁=−𝑒𝑏(1−𝑛₂)/𝐷₂(𝑏+𝑛₂)+𝑠(𝑏+𝑛₂)(1−𝑛₂)/𝐷₂. Thus, the characteristic equation of (𝑛₂,0,𝑝₂,0) is given by |𝜂𝐸−𝐽_{(𝑛2,0,𝑝2,0)}|=0 and 𝜂 is its eigenvalues, as follows: 𝑇(𝜂)=𝜂⁴−(𝑐/𝐷₁+𝑐/𝐷₂)𝜂³+(𝑐²/𝐷₁𝐷₂−𝑎₂₁)𝜂²+(𝑐𝑎₂₁/𝐷₂)𝜂−𝑎₂₃𝑎₄₁.

Based on the Routh-Hurwitz stability criterion, when 𝑎₄₁<0 and 𝑐>𝐷₁𝐷₂𝑎₂₁, the characteristic equation of (𝑛₂,0,𝑝₂,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 (𝑛₂,0,𝑝₂,0) and a trajectory from (1,0,0,0) to (𝑛₂,0,𝑝₂,0) where 𝑠→∞. Thus, if 𝑎₄₁<0 and 𝑐>max(2𝐷₂(𝑒/(𝑏+1)−𝑑),𝐷₁𝐷₂𝑎₂₁) 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 (𝐿𝑥,𝐿𝑦)=(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 𝑁, 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 𝑑=0.1,   𝑠=0.01,  𝑒=0.4,  𝐷₁=1, and 𝐷₂=0.01.

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 0.1 to 0.5. Figures 1(a) and 1(b) show that the value of 𝑏 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 𝑏. In Figures 1(c) and 1(d), the value of 𝑏 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 𝑏 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 𝑏. 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.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1 

Spatiotemporal dynamics of model (3). In ((a), (b)) 𝑏=0.1, ((c), (d)) 𝑏=0.2, and ((e), (f)) 𝑏=0.5. The iteration numbers are 140 in ((a), (c), (e)) and 3000 in ((b), (d), (f)).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2 

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.

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.

(a)

(b)

(c)

(d)

Figure 3 

Spatiotemporal dynamics of model (3) where 𝑏=1.2. The number of iterations are as follows (a) 1, (b) 200, and (c) 3000. (d) The phase diagram.

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 𝐸₁=(1,0) to 𝐸₂=(𝑛₂,𝑝₂) with spatial variation 𝑥 and time variation 𝑡, before finally reaching a steady state. We also used different time traveling waves 𝐿₁, 𝐿₂, 𝐿₃, 𝐿₄ at times 𝑡=3500, 𝑡=4000, 𝑡=4500, 𝑡=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 𝐸₁=(1,0), 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 𝐸₁=(1,0) and the coexisting equilibrium points 𝐸₂=(𝑛₂,𝑝₂). These results further show that the coexisting equilibrium points 𝐸₂=(𝑛₂,𝑝₂) becomes more stable as the value of 𝑏 increase.

(a)

(b)

Figure 4 

Traveling wave solutions for phytoplankton and the change in phytoplankton density at specific times.

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, 𝑚₁=1−𝑏₁/𝑏𝑘₁. 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).