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

The spatiotemporal dynamics of a phytoplankton-zooplanktonmodelwith an alternative prey and refuge effect is investigatedmathematically and numerically.The stability of the equilibriumpoint and the travelingwave 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.


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][5][6][7][8][9][10][11][12][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 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

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:  1 represents the inherent growth rate of the phytoplankton without any environment limitations;  1 is the carrying capacity of the phytoplankton in the presence of predators and harvesting;  1 represents the consumption rate;  1 is the half-saturation constant that determines how quickly the maximum is attained when the phytoplankton density increases;  1 is the conversion factor denoting the number of new herbivores produced for each capture;  1 is the growth rate of micrograzers due to the alternative prey;  1 denotes the foodindependent zooplankton mortality rate where  1 >  1 ;  1 is the refuge protect the prey, where  1 ∈ [0, 1) is constant [29][30][31][32].For simplicity, using the scaling  =  1 ,  = ( 1  1 /  1 ), and  = / 1 , we have where 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 ∇ 2 =  2 / 2 + 2 / 2 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  1 and  2 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  2 + .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].

Mathematical Analysis
First, we will prove that all solutions of model ( 3) are bounded.
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   () = − 2 () and ( * ) = max ∈Ω (,  * ).
Among the conditions that satisfy Theorem 1, we mainly address the steady states and their stabilities for (3), as follows.
(i)  0 = (0, 0), which corresponds to the extinction of the phytoplankton and zooplankton, which is apparently an unstable saddle point.
(iii) Their coexistence will be observed in different states  2 = ( 2 ,  2 ) and  3 = ( 3 ,  3 ), where and 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 into model (3), where , and selecting all terms that are linear about , where where Substituting ( 19) into (18), we obtain: (21)

Stability Analysis
Proof.The stability of the equilibrium state is determined based on the nature of the eigenvalues in the Jacobian matrix (  ).Substituting  0 into (19), we can determine the eigenvalues of the Jacobian matrix ( 0 ), 1 and  − .Apparently  −  < 0, so  0 is unstable.The stability of  1 can be established in a similar manner as  2 , so we only discuss the stability of  2 .
We can use   to determine the eigenvalues of the operator −Δ on Ω, which satisfy the homogeneous Neumann boundary condition and By substituting  2 = ( 2 ,  2 ) 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   = ( 2 ) −   .
Thus,  2 = ( 2 ,  2 ) 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 tr (  ) < 0 and det (  ) > 0 is satisfied, the above conditions are met.Therefore, with a random   ≥ 0, tr (  ) < 0 and det (  ) > 0 hold if  11 < 0 and  12  21 < 0.

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 = 1, and  2 = 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.
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  1 = (1, 0) to  2 = ( 2 ,  2 ) with spatial variation  and time variation , before finally reaching a steady state.We also used different time traveling waves  1 ,  2 ,  3 ,  4 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 = (1, 0), before finally moving to a steady state  2 = ( 2 ,  2 ).Furthermore, it is important to note that model (3) has a traveling wave solution connecting the steady state  1 = (1, 0) and the coexisting equilibrium points  2 = ( 2 ,  2 ).These results further show that the coexisting equilibrium points  2 = ( 2 ,  2 ) becomes more stable as the value of  increase.

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 = 1 −  1 / 1 .When the carrying capacity  1 and half-saturation constant  1 are fixed, the value of  is only determined by the refuge protecting the prey  1 , which is changed readily by different biological and chemical factors.These result shows that the prey refuge  1 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.

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