Effect of Prey Refuge on the Spatiotemporal Dynamics of a Modified Leslie-Gower Predator-Prey System with Holling Type III Schemes

In this paper, the spatiotemporal dynamics of a diffusive Leslie-Gower predator-prey model with prey refuge are investigated analytically and numerically. Mathematical theoretical works have considered the existence of global solutions, population permanence and the stability of equilibrium points, which depict the threshold expressions of some critical parameters. Numerical simulations are performed to explore the pattern formation of species. These results show that the prey refuge has a profound effect on predator-prey interactions and they have the potential to be useful for the study of the entropy theory of bioinformatics.


Introduction
The dynamic relationship between predators and their prey has fascinated mathematical biologists for a long time.A variety of mathematical models are devoted to exploring the predator-prey interaction [1][2][3][4].To understand well the population dynamics, many biological factors are included such as time delay, impulsive effect, seasonal perturbation [5][6][7][8][9].Recently, many authors [10,11] have focused on the dynamics of a class of the semi-ratio-dependent predator-prey models, in which one of the salient features is that the carrying capacity of predator is proportional to the number of prey and such models OPEN ACCESS were initially introduced by Leslie and Gower [12,13].In 2003, Aziz-Alaoui and Okiye [14] analyzed the dynamics of the following model: where u and w represent the densities of prey and predator, respectively.Furthermore, it is assumed that the prey grows logistically with the limited factor k of considering realistic surroundings and innate growth rate r .In Equation ( 1), 1 k is the average saturation rate, which indicates the quality of the food that provides prey to predator, 2 k indicates the quality of the alternative that provides the environment, s is the intrinsic growth rate of predator, e is the maximum reduction of prey due to predation and h measures the ration of prey to support one predator.Here the functional response of predator is Holling type II schemes, which usually depicts the uptake of substrate by the microorganisms in microbial kinetics [15].Oftentimes Holling type III schemes is used to describe the dynamical behavior of the invertebrate feeding on the prey and this functional response of predator has been widely included in mathematic ecological models [16][17][18][19].In fact, if the predator is the invertebrate, Holling type III functional response can fit better [20].On the other hand, the effect of a constant proportion of prey refuge on predator-prey models has become a pretty hot issue in mathematical ecology in the recent years.By investigating the theoretical models, most of theoretical conclusions show that the prey refuge has a stabilizing effect on predator-prey systems, but the dynamics of the Kolmogorov type model incorporating a constant proportion of prey refuge is qualitatively equivalent to the original system [21][22][23][24][25][26].Thus, we consider the following system: )) is a constant and u m) 1 (  reflects the prey available to the predator, a is the half-saturation constant for the predator and b indicates the quality of the alternative that provides the environment. On the other hand, all living beings live in a spatial world, which can cause that the spatial component of ecological interactions exhibits ranging from individual behavior to species abundance, diversity and population dynamics.Therefore, the spatial factor is one of the most important elements in ecosystem.Lately, Camera [27] has specified the spatiotemporal dynamics of Equation (1) with diffusion of species.Meanwhile, a large amount of literatures mainly study this theme in reaction-diffusion systems since Turing [28] pointed out that this kind of system could yield many complex patterns, which are usually consistent with a wide variety of phenomena that have been observed in chemistry, physics and biology [29][30][31].Thus, Equation (2) with the spatial factor can be described as following:  all the parameters are assumed to be positive.The rests of the paper are structured as follows: in Section 2, the existence of the global solutions and the population permanence of Equation (3) are proved.In Section 3, the local stability of the equilibrium points and the global stability of the interior equilibrium point are investigated.Furthermore, the Turing instability and the conditions of its occurrence are analyzed.In Section 4, under the condition of Turing instability, numerical simulations are illustrated to show how the prey refuge affects spatiotemporal dynamics of Equation (3).In the end, some discussions are given.

Existence of Global Solutions
, there is a unique global solution of Equation ( Proof: Equation ( 3) is mixed quasi-monotone since: is the unique solution of: (1 ) (0) , (0) is any solution of Equation (3), then there is: , then there is is any solution of: * ( 1), ( 1)( ) ( ) Similarly, there exists . This completes the proof. where then it has . By Theorems 1 and 2, there exists 0 2  T such that: then it has . This completes the proof.

Remark 1.
From Theorem 2 and Theorem 3, it is clear that Equation (3) is permanent.

Stability
It is clear that Equation (3) into Equation ( 3) and picking up all the terms which are linear in  , there is: where ) ( ( 1) ) ( 1) Proof: From above, the linearized result of Equation (3) around 0 E is: then it needs to consider the largest eigenvalue of: Assume that  is an eigenvalue of Equation ( 22) with the eigenfunction ) , ( 21 z z and 0 1  z , then  is an eigenvalue of r d   1 with homogeneous Neumann boundary condition.Furthermore, it follows that  must be real.In the same way,  is also real provided that 0 2  z .Then all eigenvalues of Equation ( 22) must be real.Let max  denote the largest eigenvalue.Consider the principal eigenvalue  ˆ of: then it shows that its principal eigenvalue  ˆ is positive and the associated eigenfunction 22), then it satisfies Equation ( 22) with    .Thus, it is clear that 0 ˆ  is an eigenvalue of Equation ( 22), and there is 0 max 

  
. This exhibits that 0 E is unstable.
Proof: From Equation ( 20), the linearized result of Equation ( 3) around 1 E is: (1 ) ( 1) then it needs to consider the largest eigenvalue of: ) As the previous case, all eigenvalues of Equation ( 25) are real.Assume that max  is the largest eigenvalue of Equation (25).Consider the principal eigenvalue  ˆ of: then it shows that its principal eigenvalue  ˆ is positive and the associated eigenfunction 0 ˆ2  z .Furthermore, assume that 1 ẑ which is positive, is the solution of: ) satisfies Equation ( 25) with 0 ˆ

  
. Thus there is 0 max 

  
. This exhibits that 1 E is unstable.
Similarly, it can be concluded that E 2 is unstable.
, then * E is locally asymptotically stable.
Proof: Assume that on  with homogeneous Neumann boundary condition and i  denote the associated eigenfunction corresponding to i  , then there is: Furthermore, the linearized result of Equation ( 3) around * E is: where ) , ( , where where: , and is the solution of the above system.Finally, , then * E is globally asymptotically stable.
Proof: Consider the Lyapunov function: . The orbital derivative of V along the solutions of Equation ( 3) is: , then: Using the result of Theorem 2, there exist Applying Lemma 1, there is: By Theorem 3, there exist 3  and T such that for T t  : From Equation (34), dt d is bounded, where: Using the Poincare inequality, there exists: In fact, there exist: From Equations ( 39) and (41), it result in: From Inequality (34), there exists a subsequence of } { m t which is also denoted } { m t , and nonnegative functions Combined with (42), we obtain: This above result and the local stability conditions can yield that 2 E is globally asymptotically stable.

Turing Instability
In order to investigate the transition of the equilibrium state, we consider small space-and time-dependent perturbations for any solution of Equation ( ) where  ,  are small enough, k is the wave number.Substituting Equation (44) into Equation (3), we linearize the system around * E and further obtain its characteristic equation: where: From Equation (45), the dispersion relation of Equation ( 3) is: Turing instability requires that the stable interior equilibrium point is driven unstable by the local dynamics and diffusion of species.The conditions for the homogeneous state of Equation ( 2) to be stable is 0 . It is clear that 0 tr tr k  .Then the stability of the homogeneous state simultaneously changes the sign of k  .From Equation (46), it easily finds that there is ,where:  have positive values, we can obtain the range of instability for a local stable equilibrium, which is called as the Turing space.In order to show the Turing space, the dispersion relation is plotted corresponding to several values of the bifurcation parameter m while in Figure 1  It should be stressed from Figure 1 that the available Turing modes are further reduced when the value of prey refuge m is increasing.Nonetheless, it is interesting to notice that Equation (3) will occur the Turing instability when the value of m less than 2512032 .0 .

Turing Pattern Formation
To better investigate how the prey refuge affects the spatiotemporal dynamics of Equation ( 3 space units per time unit.As the initial perturbation propagates, Equation (3) under the condition of Turing instability evolves a steady state, which is stationary in time and oscillatory in space.Moreover, it should be stressed that the spatial patterns of predator and prey under the condition of Turing instability are always the same type, this is because that it is assumed that the carrying capacity of predator is proportional to the number of prey, and the steady state of predator is equal to this carrying capacity.Thus, only the spatial patterns of prey are shown.
It is interesting to note from Figure 2 that some snapshots have been taken of numerical simulations when the value of m increases from 0 to 35 .0 .It should be pointed out that in these snapshots the enclosed color bars denote the range of the changing densities of prey, where higher values correspond to higher prey densities.Figure 1  stand for the stable spatiotemporal behavior.By comparing the first four diagrams, it can be observed that the spatiotemporal dynamical behaviors of Equation (3) are very rich and complex.When the value of m is 0 , the spatial distribution of prey is mainly some interconnected strips and nonuniform, which shows that the habits of prey are the main type of community survival, so it is easy to evade predator-capturing.When the values of m are 08 .0 and 15 .0 , the collective survival population expands gradually and the spatial distribution of prey tends to be uniform.When the value of m is 22 .0 , the spatial distribution of prey is almost uniform and the prey can survive in any space.On the other hand, from Figure 2 the maximum values on color bars exhibit decreasing states as the effect of prey refuge is strengthened.Inversely, the interior equilibrium density value of prey will increase as the increase of m .In order to relieve the crowed space, the competitive pressure between individuals of prey is intensified.From the biological point of view, the effect of prey refuge may be to help prey relieve the pressure of predation during diffusion.Thus, the patches of high density prey diffuse into the low.Finally, the distributions of prey tend to be uniform as the effect of prey refuge increases.However, when the value of m is more than 2512032 .0 , the prey and the predator will be involved into a stable state, so the prey can live in any space, which can be shown in behind two diagrams of Figure 2.These results show that the prey refuge not only promotes an increase in the number of prey, but also is conducive to their living space extension.
For further analysis of the effect of prey refuge on the dynamical behavior of one population, the spatiotemporal evolutions of prey have been obtained at 100  x , which correspond to Figure 2. It should be stressed from Figure 3 that these results are consistent with Figure 2, which show the accuracy and effectiveness of numerical simulations.Moreover, the comparison of the first four diagrams in Figure 3 suggests that when the value of m gradually increases and is close to 2512032 .0 , oscillations in space diminish gradually.These results show that a suitable prey refuge has a positive effect on predator-prey interactions.It is easy to see that if the effect of prey refuge is strengthened in living surroundings, predation risk is relatively reduced in the habitat and consequently the density of prey is bound to increase.And the densities of predator and prey will obtain the new balance.
Based on the above analysis, it can be seen that a suitable prey refuge can enhance the specie biomass level and promote the uniformness of the population distribution, which agree with some results of the real world.Furthermore, it is interesting to point out that the lower value of prey refuge can come into rich spatiotemporal dynamics.Moreover, the use of mathematical model with a prey refuge and diffusion is considered to explore some biological problems, and the numerical simulation can provide an approximation of the real biological behaviors.Hence, these results can promote the study of ecological patterns.  .Other parameters are fixed as Equation (49).

Conclusions
In this paper, a diffusive predator-prey system with Holling type III scheme has been studied analytically and numerically.Mathematical theoretical works have considered the existence of global solutions and the stability of equilibrium points and population permanence.On the basis of these results, we obtain the threshold expressions of some critical parameters which in turn provide a theoretical basis for the numerical simulation.Numerical simulations indicate that the prey refuge has a strong and positive effect on the spatiotemporal dynamics according to the spatial patterns and spatiotemporal evolution of prey.Furthermore, it should be stressed that the spatial pattern diagrams show that the prey refuge has a profound effect on predator-prey interactions.Using the spatiotemporal evolution of prey, the spatial distribution of prey and the accuracy effectiveness of numerical simulation can be further confirmed.All these results are expected to be of significance in the exploration of the entropy theory of bioinformatics.

Figure 1 .
Figure 1.Variation of dispersion relation of Equation (3) around the interior equilibrium point.The red line corresponds to 08 .0  m , the green is 25 .0  m and the blue is 35 .0  m .
), the spatial distribution diagrams are obtained as change of m .All numerical simulations are carried out in a discrete two-dimensional domain with 200 200  lattice sites.The step between each lattice point is defined as 25 .0   .The time evolution of Equation (3) is resorted to the forward Euler integration with a
clearly shows that Equation (3) leads to the Turing instability for