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

**Jianglin Zhao 1 , Min Zhao 2,\* and Hengguo Yu 1**


*Received: 27 March 2013; in revised form: 5 June 2013 / Accepted: 6 June 2013 / Published: 19 June 2013* 

**Abstract:** 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.

**Keywords:** diffusive predator-prey system; Leslie-Gower; Holling type III schemes; refuge; stability; pattern formation

### **1. 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–4]. To understand well the population dynamics, many biological factors are included such as time delay, impulsive effect, seasonal perturbation [5–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 were initially introduced by Leslie and Gower [12,13]. In 2003, Aziz-Alaoui and Okiye [14] analyzed the dynamics of the following model:

$$\begin{cases} \frac{du}{dt} = ru(1 - \frac{u}{k}) - \frac{euw}{u + k\_1}, \\\\ \frac{dw}{dt} = sw(1 - \frac{hw}{u + k\_2}), \end{cases} \tag{1}$$

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), 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–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–26]. Thus, we consider the following system: 1 *k*

$$\begin{cases} \frac{du}{dt} = ru(1 - \frac{u}{k}) - \frac{e(1 - m)^2 u^2 w}{a + (1 - m)^2 u^2} \\\\ \frac{dw}{dt} = sw(1 - \frac{hw}{b + (1 - m)u}) \end{cases} \tag{2}$$

where *mm* ))1,0[( is a constant and )1( *um* 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–31]. Thus, Equation (2) with the spatial factor can be described as following:

**186** 

$$\begin{cases} \frac{\partial u}{\partial t} = d\_1 \Delta u + ru(1-\frac{u}{k}) - \frac{e(1-m)^2 u^2 w}{a + (1-m)^2 u^2} = d\_1 \Delta u + \phi(u, w), t > 0, \, x \in \Omega\\ \frac{\partial w}{\partial t} = d\_2 \Delta w + s\nu(1 - \frac{h\nu}{b + (1-m)u}) = d\_2 \Delta w + \phi(u, w), t > 0, \, x \in \Omega\\ \frac{\partial u}{\partial n} = \frac{\partial w}{\partial n} = 0, t > 0, \, x \in \partial\Omega,\\ u(0, \mathbf{x}) = u\_0(\mathbf{x}), w(0, \mathbf{x}) = w\_0(\mathbf{x}), \, x \in \overline{\Omega} \end{cases} \tag{3}$$

where *xtu* ),( and *xtw* ),( denote the densities of prey and predator at time *t* and position *x* , respectively. is the Laplacian operator, 1 *d* and 2 *d* are the diffusion coefficients of prey and predator, and *<sup>n</sup>* is differentiation in the direction of the outward unit normal to : . In Equation (3), 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.

#### **2. Existence of Global Solutions and Permanence**

#### *2.1. Existence of Global Solutions*

**Theorem 1.** For 0)(0 *xu* \$ , 0)(0 *xw* \$ , there is a unique global solution of Equation (3) such that *xtu*  0),( , *xtw*  0),( for *t*  0 and *x* : .

**Proof:** Equation (3) is mixed quasi-monotone since:

$$\frac{d\phi}{dw} = -\frac{e(1-m)^2 u^2}{a + (1-m)^2 u^2} \le 0\tag{4}$$

$$\frac{d\phi}{du} = \frac{hs(1-m)w^2}{\left(b + (1-m)u\right)^2} \ge 0\tag{5}$$

in }0,0{ <sup>2</sup> *wuR* \$\$ .

Consider that (ˆ, *wu* ˆ) is the unique solution of:

$$\begin{cases} \frac{du}{dt} = ru(1 - \frac{u}{k})\\ \frac{d\nu}{dt} = s\nu(1 - \frac{h\nu}{b + h(1 - m)\mu})\\ u(0) = u^\*, \ w(0) = w^\* \end{cases} \tag{6}$$

where 0 )(sup \* *xuu* : , 0 )(sup \* : *xww* .

Let *xtwxtu* )0,0()),(),,(( and ()),(),,(( ˆ, *wuxtwxtu* ˆ) . There exist:

$$\frac{\partial \overline{\boldsymbol{u}}}{\partial t} - d\_l \Delta \overline{\boldsymbol{u}}(t, \mathbf{x}) - \boldsymbol{\varphi}(\overline{\boldsymbol{u}}(t, \mathbf{x}), \underline{\mathbf{u}}(t, \mathbf{x})) = 0 \geq 0 \\ = \frac{\partial \underline{\boldsymbol{u}}}{\partial t} - d\_l \Delta \underline{\boldsymbol{u}}(t, \mathbf{x}) - \boldsymbol{\varphi}(\underline{\mathbf{u}}(t, \mathbf{x}), \overline{\boldsymbol{\omega}}(t, \mathbf{x})) \tag{7}$$

$$\frac{\partial \overline{\boldsymbol{w}}}{\partial t} - d\_2 \Delta \overline{\boldsymbol{\nu}}(t, \mathbf{x}) - \mathbf{g}(\overline{\boldsymbol{\mu}}(t, \mathbf{x}), \overline{\boldsymbol{\nu}}(t, \mathbf{x})) = 0 \geq 0 = \frac{\partial \underline{\boldsymbol{w}}}{\partial t} - d\_2 \Delta \underline{\boldsymbol{\nu}}(t, \mathbf{x}) - \mathbf{g}(\underline{\boldsymbol{\mu}}(t, \mathbf{x}), \underline{\boldsymbol{\nu}}(t, \mathbf{x})) \tag{8}$$

Clearly, the boundary conditions are satisfied. Then *xtwxtu* )),(),,(( and *xtwxtu* )),(),,(( are the lower-solution and upper-solution of Equation (3), respectively. Thus, Equation (3) has a unique global solution, which can satisfy ),(0 ## ˆ *tuxtu* )( , ),(0 ## ˆ *twxtw* )( for *t* \$ 0 .

#### *2.2. Permanence*

**Definition 1.** Equation (3) is said to be permanence if for any solution with nonnegative initial functions )(0 *xu* and )(0 *xw* )0)(,0)(( <sup>0</sup> <sup>0</sup> *xwxu* ;; , there exist positive constants *mi* and *Mi* ( *i* 2,1 ) such that:

$$\|m\_1 \le \liminf\_{t \to \infty} u(t, \mathbf{x}) \le \limsup\_{t \to \infty} u(t, \mathbf{x}) \le M\_1 \tag{9}$$

$$m\_2 \le \liminf\_{t \to \infty} \,\,\omega(t, \mathbf{x}) \le \limsup\_{t \to \infty} \,\omega(t, \mathbf{x}) \le M\_2 \tag{10}$$

**Theorem 2.** For any solution *xtwxtu* )),(),,(( of Equation (3):

$$\limsup\_{t \to \infty} u(t, \mathbf{x}) \le k, \quad \limsup\_{t \to \infty} w(t, \mathbf{x}) \le \frac{b + k(1 - m)}{h} \tag{11}$$

**Proof:** Suppose that *xtwxtu* )),(),,(( is any solution of Equation (3), then there is:

$$\begin{cases} \frac{\partial u}{\partial t} \le d\_1 \Delta u + ru(1 - \frac{u}{k}), t > 0, x \in \Omega\\ \frac{\partial u}{\partial n} = 0, t > 0, x \in \partial\Omega\\ u(0, x) = u\_0(\mathbf{x}) \le u^\* = \sup\_{\overline{\Omega}} u\_0(\mathbf{x}), \mathbf{x} \in \overline{\Omega} \end{cases} \tag{12}$$

Furthermore, it is assumed that *tU* )( is any solution of:

$$\begin{cases} d\hat{U} = r\hat{U}(1 - \frac{\hat{U}}{k}), \; t > 0\\ \hat{U}(0) = u \; ^\ast \le k \end{cases} \tag{13}$$

then there is *ktU <sup>t</sup>* # %" )(lim . By the comparison principle, there exists *kxtu <sup>t</sup>* # %" ),(suplim . As a consequence, for any <  0 , there exists *T*  0 such that ),( *kxtu* -# <for  *Tt* . Then there is:

$$\begin{cases} \frac{\partial \mathbf{w}}{\partial t} \le d\_2 \Delta \mathbf{w} + s \nu (1 - \frac{h \mathbf{w}}{b + (1 - m)(k + \varepsilon)}), t > T, \mathbf{x} \in \Omega\\ \frac{\partial \mathbf{w}}{\partial n} = \mathbf{0}, t > T, \mathbf{x} \in \partial \Omega\\ \nu(T, \mathbf{x}) \ge 0, \mathbf{x} \in \overline{\Omega} \end{cases} \tag{14}$$

Consider that *tW* )( is any solution of:

**188** 

$$\begin{cases} d\hat{W} = s\hat{W}(1 - \frac{h\hat{W}}{b + (1 - m)(k + \varepsilon)}), t > T\\ \hat{W}(T) = \boldsymbol{w}^\* \end{cases} \tag{15}$$

then it leads to *h kmb tW <sup>t</sup>* ))(1( )(lim -- # %" 

$$\limsup \,\mathsf{w}(t,\mathsf{x}) \le \frac{b + (1 - m)(k + \varepsilon)}{\iota}$$

<

.

Similarly, there exists *h t* %" . As a result, let < " <sup>0</sup> , it has *<sup>h</sup> mkb xtw t* )1( ),(suplim - # %" . This completes the proof.

**Theorem 3.** Assume that ))1(()1( <sup>2</sup> - *mkbkmehr* , then:

$$\liminf\_{t \to \infty} w(t, x) \ge \frac{b}{h} > 0, \liminf\_{t \to \infty} u(t, x) \ge k\eta > 0 \tag{16}$$

where *ahr mkbkmea* ))1(()1( <sup>1</sup> <sup>2</sup> - =.

**Proof:** Assume that *wu* ),( is any solution of Equation (3) with 0)(0 *xu* \$ and 0)( <sup>0</sup> *xw* \$ )0)(,0)(( <sup>0</sup> <sup>0</sup> *xwxu* ;; . Meanwhile, there is *xtu* \$ 0),( for all *<sup>t</sup>* \$ <sup>0</sup> . Then there is *b hw umb hw* \$ - <sup>1</sup> )1( 1 for all *<sup>t</sup>* \$ <sup>0</sup> , *x* : . Assume that 0),( <sup>1</sup> *xTw*  for 0 *<sup>T</sup>*<sup>1</sup>  . Consider that 

*W* is any solution of:

$$\begin{cases} d\breve{W} = s\breve{W}(1 - \frac{h\breve{W}}{b}), t > T\_1\\ \breve{W}(T\_1) = \inf\_{x \in \widetilde{\Omega}} \varkappa(T\_1, x) > 0 \end{cases} \tag{17}$$

then it has *h <sup>b</sup> tW <sup>t</sup>* %" )(lim . By the comparison principle, there is \$ 0),(inflim %" *<sup>h</sup> b xtw t* for *<sup>x</sup>* : . By Theorems 1 and 2, there exists 0 *T*<sup>2</sup>  such that:

$$r(1 - \frac{u}{k}) - \frac{ea(1 - m)^2 u \nu}{a + (1 - m)^2 u^2} \ge r(1 - \frac{u}{k}) - \frac{ea(1 - m)^2 k(b + k(1 - m))}{ah} \tag{18}$$

for 2 , *xTt* : \$ .

Let *ahr mkbkmea* ))1(()1( <sup>1</sup> <sup>2</sup> - =. Consider that *tU* )( is any solution of:

$$\begin{cases} d\breve{U} = r\breve{U}(\eta - \frac{\breve{U}}{k}), t > T\_2\\ \breve{U}(T\_2) = \inf\_{x \in \Omega} u(T\_2, x) > 0 \end{cases} \tag{19}$$

then it has *ktU* = *t* %" )(lim . Similarly, there is *kxtu* = *t* \$ %" ),(inflim . According to the assumption in the theorem, it can yield \$ 0),(inflim %" *kxtu* = *<sup>t</sup>* . This completes the proof.

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

#### **3. Stability Analysis of Equilibrium Points and Turing Instability**

#### *3.1. Stability*

It is clear that Equation (3) has the following equilibrium points: (a) )0,0( *E*<sup>0</sup> (total extinction), (b) )0,( <sup>1</sup> *kE* (extinction of the predator), (c) ),0( <sup>2</sup> *<sup>h</sup> <sup>b</sup> <sup>E</sup>* (extinction of the prey), (d) ),( *wuE* \*\*\* (coexistence of prey and predator), where ),( *wu* \*\* is the positive solution of 9 *wu* 0),( , 8 *wu* 0),( .

In order to investigate the linear stability of equilibrium solutions 2,1,0 )(*E ii* and *E*\* of Equation (3), we consider the corresponding eigenvalue problem of the linearized operator around every equilibrium point.

Substituting )),(),,((),()),(),,(( <sup>1</sup> <sup>2</sup> -> *xtzxtzExtExtwxtu* into Equation (3) and picking up all the terms which are linear in > , there is:

$$\frac{\partial \mathcal{Z}}{\partial t} = D \Delta \mathcal{Z} + L(E) \mathcal{Z} \tag{20}$$

where 2 1 0 0 *d d D* and

$$L(E) = \begin{pmatrix} r - \frac{2ru}{k} - \frac{2ea(1-m)^2 uw}{\left(a + (1-m)^2 u^2\right)^2} & -\frac{e(1-m)^2 u^2}{a + (1-m)^2 u^2} \\ & \frac{sh(1-m)w^2}{\left(b + (1-m)u\right)^2} & s(1 - \frac{2hw}{b + (1-m)u}) \end{pmatrix}\_E$$

**Proposition 1.** *E*0 is unstable.

**Proof:** From above, the linearized result of Equation (3) around *E*0 is:

$$\begin{cases} \frac{\partial \boldsymbol{z}\_1}{\partial t} = d\_1 \Delta \boldsymbol{z}\_1 + r \boldsymbol{z}\_1 \\ \frac{\partial \boldsymbol{z}\_2}{\partial t} = d\_2 \Delta \boldsymbol{z}\_2 + s \boldsymbol{z}\_2 \\ \frac{\partial \boldsymbol{z}\_1}{\partial n} \bigg|\_{\partial \Omega} = \frac{\partial \boldsymbol{z}\_2}{\partial n} \bigg|\_{\partial \Omega} = \mathbf{0} \end{cases} \tag{21}$$

then it needs to consider the largest eigenvalue of:

$$\begin{cases} d\_1 \Delta \nu\_1 + r \mathbf{z}\_1 = \lambda \mathbf{z}\_1 \\ d\_2 \Delta \nu\_2 + s \mathbf{z}\_2 = \lambda \mathbf{z}\_2 \\ \left. \frac{\partial \mathbf{z}\_1}{\partial n} \right|\_{\partial \Omega} = \left. \frac{\partial \mathbf{z}\_2}{\partial n} \right|\_{\partial \Omega} = 0 \end{cases} \tag{22}$$

Assume that ? is an eigenvalue of Equation (22) with the eigenfunction ),( <sup>21</sup> *zz* and *z*<sup>1</sup> ; 0 , then ? is an eigenvalue of 1 *rd* with homogeneous Neumann boundary condition. Furthermore, it follows that ? must be real. In the same way, ?is also real provided that *z*<sup>2</sup> ; 0 . Then all eigenvalues of Equation (22) must be real. Let ? max denote the largest eigenvalue. Consider the principal eigenvalue ? ˆ of:

$$\begin{cases} d\_1 \Delta z\_1 + r z\_1 = \lambda z\_1 \\ \left. \frac{\partial z\_1}{\partial n} \right|\_{\partial \Omega} = 0 \end{cases} \tag{23}$$

then it shows that its principal eigenvalue ? ˆ is positive and the associated eigenfunction *z*ˆ <sup>1</sup>  0 . Let us substitute (),( ˆ )0,121 *zzz* into Equation (22), then it satisfies Equation (22) with ?? ˆ . Thus, it is clear that 0 ˆ ?  is an eigenvalue of Equation (22), and there is 0 <sup>ˆ</sup> max ?? \$ . This exhibits that *E*0 is unstable.

**Proposition 2.** *E*1 is unstable.

**Proof:** From Equation (20), the linearized result of Equation (3) around *E*1 is:

$$\begin{cases} \frac{\partial \mathbf{z}\_1}{\partial t} = d\_1 \Delta \mathbf{z}\_1 - r \mathbf{z}\_1 - \frac{e(1-m)^2 k^2}{a + (1-m)^2 k^2} \mathbf{z}\_2\\ \frac{\partial \mathbf{z}\_2}{\partial t} = d\_2 \Delta \mathbf{z}\_2 + s \mathbf{z}\_2\\ \left. \frac{\partial \mathbf{z}\_1}{\partial n} \right|\_{\varepsilon \Omega} = \left. \frac{\partial \mathbf{z}\_2}{\partial n} \right|\_{\varepsilon \Omega} = 0 \end{cases} \tag{24}$$

then it needs to consider the largest eigenvalue of:

$$\begin{cases} d\_1 \Delta z\_1 - r z\_1 - \frac{e(1-m)^2 k^2}{a + (1-m)^2 k^2} z\_2 = \lambda z\_1 \\ d\_2 \Delta z\_2 + s z\_2 = \lambda z\_2 \\ \left. \frac{\partial z\_1}{\partial n} \right|\_{\partial \Omega} = \left. \frac{\partial z\_2}{\partial n} \right|\_{\partial \Omega} = 0 \end{cases} \tag{25}$$

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:

$$\begin{cases} d\_2 \Delta z\_2 + s z\_2 = \lambda z\_2\\ \left. \frac{\partial z\_2}{\partial n} \right|\_{\partial \Omega} = 0 \end{cases} \tag{26}$$

then it shows that its principal eigenvalue ? ˆ is positive and the associated eigenfunction ˆ 0 *z*<sup>2</sup>  .

Furthermore, assume that 1*z*ˆ which is positive, is the solution of:

$$\begin{cases} d\_1 \Delta z\_1 - (r + \lambda) z\_1 = \frac{e(1 - m)^2 k^2}{a + (1 - m)^2} z\_2\\ \left. \frac{\partial z\_1}{\partial n} \right|\_{\partial \Omega} = 0 \end{cases} \tag{27}$$

then (ˆ , ˆ )21 *zz* satisfies Equation (25) with 0 ˆ ??  . Thus there is 0 <sup>ˆ</sup> max ?? \$ . This exhibits that *E*1 is unstable.

Similarly, it can be concluded that *E*2 is unstable.

**Proposition 3.** Assume that 1 )1( 2 \* \* 2 \* <sup>2</sup>  - *uk u uma <sup>a</sup>* , then *E*\* is locally asymptotically stable. **Proof:** Assume that 0 ? ? ?210 ?*<sup>i</sup>* @@@@@ denote the eigenvalues of on : with homogeneous Neumann boundary condition and 9*<sup>i</sup>* denote the associated eigenfunction

corresponding to ?*<sup>i</sup>* , then there is:

$$\begin{cases} -\Delta \phi\_i = \lambda\_i \phi\_i\\ \left. \frac{\partial \phi\_i}{\partial n} \right|\_{\partial \Omega} = 0 \end{cases} \tag{28}$$

Furthermore, the linearized result of Equation (3) around *E*\* is:

$$\frac{\partial \mathcal{Z}}{\partial t} = D \Delta \mathcal{Z} + \Sigma \mathcal{Z} \tag{29}$$

where ),( *dddiagD* 21 and

$$\Sigma = \begin{pmatrix} r(1 - \frac{2u\_\*}{k}) - \frac{2ae(1 - m)^2 u\_\* w\_\*}{\left(a + (1 - m)^2 u\_\*^{\frac{2}{2}}\right)^2} & -\frac{e(1 - m)^2 u\_\*^{\frac{2}{2}}}{a + (1 - m)^2 u\_\*^{\frac{2}{2}}}\\ \frac{s(1 - m)}{h} & -s \end{pmatrix} = \begin{pmatrix} A & B\\ \mathcal{Q} & R \end{pmatrix}$$

Let )()(),( 0 *<sup>i</sup> xtxtZ i <sup>i</sup>* 9B % , <sup>2</sup> )( *Rt* B *<sup>i</sup>* substitute into Equation (29). Equaling the coefficient of C *<sup>i</sup>* , there is )( )( *<sup>t</sup> dt td ii i* B B D , where *<sup>i</sup>* AD ?*iD* . Thus, *E*3 is locally asymptotically stable for Equation (3) if and only if each *t*)( B*<sup>i</sup>* decays to zero as *t* -%" . Then, it follows that each D*i* has two eigenvalues with negative real parts, which are determined by:

$$
\mu^2 - tr(\Phi\_i)\mu + \det(\Phi\_i) = 0 \tag{30}
$$

where:

$$tr(\Phi\_i) = -(d\_1 + d\_2)\lambda\_i + A + R\tag{31}$$

$$\det(\Phi\_i) = d\_1 d\_2 \lambda\_i^2 - (R d\_1 + A d\_2) \lambda\_i + AR - BQ \tag{32}$$

Since \$ <sup>0</sup> ?*<sup>i</sup>* , *<sup>B</sup>* @ 0 and *<sup>Q</sup>*  0, it is clear that @D 0)( *<sup>i</sup> tr* and <sup>D</sup> 0)det( *<sup>i</sup>* if *<sup>A</sup>* @ <sup>0</sup> . Taking into account the assumption in Theorem, *A* @ 0 holds. This completes the proof.

For purpose of proving the global stability of *E*\* , let us introduce the following lemmas from [32].

**Lemma 1.** Let *c* and *d* be positive constants. Assume that ],[, <sup>1</sup> *ccgf* %- , *f* \$ 0 and *g* is bounded from below. If )( )( *tdf dt tdg* # and E# *dt tdf* )( in *<sup>c</sup>* %- ],[ for <sup>E</sup> <sup>0</sup> , then 0)(lim%" *tg <sup>t</sup>* .

**Lemma 2.** Consider the following equation:

**192** 

$$\begin{cases} \frac{\partial u\_i}{\partial t} = \Delta u\_i + f\_i(u\_1, \dots, u\_n), 1 \le i \le p \\ u\_i(0, \ge) = u\_0(\infty) \\ \frac{\partial u\_i}{\partial n} = 0, \quad \ge \epsilon \Omega \end{cases} \tag{33}$$

Suppose that *<sup>i</sup> xtu* ),( E@ for ##:F *niRxt* - 1,),( , and *<sup>i</sup> <sup>f</sup>* is of class <sup>1</sup> *c* on *ni <sup>n</sup> <sup>i</sup>* G ##EEA 1],,[1 , where },,,{ <sup>21</sup> *<sup>n</sup> uuuu* is the solution of the above system. Finally, suppose that *uf* 0)(*i* for *ui* <sup>0</sup> , then there exists H  0 depending only on :, , E and *<sup>i</sup> df* (the total differential of *<sup>i</sup> <sup>f</sup>* ) such that :)(,2 # H*<sup>c</sup> u* .

Assume that *wu* ),( is the unique positive solution of Equation (3). By Theorem 1, there is a constant I , which does not depend on *x* and *t*, such that *xtu* ),( # I , *xtw* ),( #I for *t* \$ 0 . By Lemma 2, there exists J 0 such that:

$$\left\|\boldsymbol{\mu}(t,\boldsymbol{x})\right\|\_{\boldsymbol{c}^{2,a}(\overline{\Omega})} \leq \mathcal{S} \,, \ \left\|\boldsymbol{\nu}(t,\boldsymbol{x})\right\|\_{\boldsymbol{c}^{2,a}(\overline{\Omega})} \leq \mathcal{S} \,\tag{34}$$

**Theorem 4.** Assume that 1 )1( 2 \* \* 2 \* <sup>2</sup> @ - *uk u uma <sup>a</sup>* , <sup>2</sup> <sup>22</sup> \* <sup>2</sup> - )1()1( *kmaum* , and 0)1()1()1()1( <sup>2</sup> \* 22 23 \* *kmaukmbkmbumab*  , then *E*\* is globally asymptotically stable. **Proof:** Consider the Lyapunov function:

$$W(\mu, \nu) = \int\_{\Omega} \mathcal{H}(\mu) d\mathbf{x} + \eta \int\_{\Omega} \mathcal{P}(\nu) d\mathbf{x} \tag{35}$$

where \* \* \* 2 ln )1( )1( ln))1(( <sup>2</sup> )1( )( *u u ub m umb umbu um <sup>u</sup>* - - -- M N , \* \* ln)( *w <sup>w</sup>* L *www* , *s ume* \* )1( = , *<sup>b</sup> uma aumbb* \* \* )1( ))1(( --- , *<sup>b</sup> uma* \* )1( N,

Then, there is N *a* @ 0 . Clearly, *V* is bounded below for all *t*  0 . The orbital derivative of *V* along the solutions of Equation (3) is:

: KK : - *dx t w dw dP dx t u du dH dt dV* = : : KKK : K: - - *dxwu dw dP dxwu du dH wdx dw dP ddxu du dH <sup>d</sup>* ),( ),( <sup>1</sup> = 2 9 8= *dxw w <sup>w</sup> ddxu umbum umb um <sup>d</sup> dt dV* <sup>2</sup> 2 \* 2 2 2 2 2 22 <sup>1</sup> ))1(()1( )1())1(( O O - - # K: K: = N K: K: P-Q - - *dxwuwwdxwu umb umauu* ),()(),( )1( ))1()(( \* \* 22 = *dxw w <sup>w</sup> ddxu umbum umb um <sup>d</sup>* <sup>2</sup> 2 \* 2 2 2 2 2 22 <sup>1</sup> ))1(()1( )1())1(( O O - - K: K: = N 

$$\begin{split} & -\int\_{\Omega} \frac{(u - u^\*)^2 (a + (1 - m)^2 u^2)}{b + (1 - m)u} \left( \frac{r}{k} - \frac{r(1 - m)^2 (u + u\_\*)(k - u\_\*)}{k \left( a + (1 - m)^2 u^2 \right)} \right) dx \\ & - \int\_{\Omega} \frac{e(1 - m)^2 v (u - u^\*)^2}{b + (1 - m)u} dx - e(1 - m)^2 u^\* \int\_{\Omega} \frac{(u - u^\*)(w - w\_\*)}{b + (1 - m)u} dx \\ & + \int\_{\Omega} \eta s (w - w\_\*) \big( -\frac{h(w - w\_\*)}{b + (1 - m)u} + \frac{(1 - m)(u - u\_\*)}{b + (1 - m)u} \big) dx \\ & \leq -d\_1 \int\_{\Omega} \frac{(\beta - \alpha)(1 - m)^2 k^2 + \beta b^2}{(1 - m)u^2 \left( b + (1 - m)u \right)^2} \left| \nabla u \right|^2 dx - d\_2 \eta \int\_{\Omega} \frac{w\_\*}{w^2} \left| \nabla v \right|^2 dx \\ & - \frac{r}{k} \int\_{\Omega} \frac{(u - u^\*)^2}{b + (1 - m)u} \big( a - (1 - m)^2 (k + u\_\*)(k - u\_\*) \big) dx - \int\_{\Omega} \frac{\eta h s (w - w\_\*)^2}{b + (1 - m)u} dx \end{split}$$

Thus, if 0)1()1()1()1( <sup>2</sup> \* 22 23 \* *kmaukmbkmbumab* and 0)1()1( <sup>2</sup> <sup>22</sup> \* <sup>2</sup> *kmaum* -, then:

$$\frac{dV}{dt} \le -\int\_{\Omega} \frac{\left(a - (1 - m)^2 (k + u\_\*) (k - u\_\*)\right)}{b + (1 - m)\mu} (\mu - u^\*)^2 d\mathbf{x} - \int\_{\Omega} \frac{\eta h \mathbf{s}}{b + (1 - m)\mu} (\mathbf{w} - \mathbf{w}\_\*)^2 d\mathbf{x} \le 0 \tag{36}$$

Using the result of Theorem 2, there exist 0 <sup>1</sup> E and 0 <sup>2</sup> E such that:

$$\frac{d}{dt}\int\_{\Omega} (u - u\_\*)^2 d\mathbf{x} = \int\_{\Omega} \mathcal{Q}(u - u\_\*) \frac{\partial u}{\partial t} d\mathbf{x} \le \mathbf{K}\_1 \tag{37}$$

$$\frac{d}{dt}\int\_{\Omega} \left(\boldsymbol{w} - \boldsymbol{w}\_{\*}\right)^{2} d\boldsymbol{x} = \int\_{\Omega} \mathcal{Z}(\boldsymbol{w} - \boldsymbol{w}\_{\*}) \frac{\partial \boldsymbol{w}}{\partial t} d\boldsymbol{x} \leq \mathbf{K}\_{2} \tag{38}$$

Applying Lemma 1, there is:

$$\lim\_{\iota\_1 \to \iota^\*} \int\_{\Omega} (\mu - \iota\_\*)^2 d\mathbf{x} = \lim\_{\iota \to \iota^\*} \int\_{\Omega} (\mu - \iota \kappa\_\*)^2 d\mathbf{x} = 0 \tag{39}$$

By Theorem 3, there exist E3 and *T* such that for  *Tt* :

$$\begin{split} \frac{dV}{dt} \leq & -d\_1 \int\_{\Omega} \frac{\left(\beta - \alpha\right) (1 - m)^2 \, k^2 + \beta b^2}{\left(1 - m\right) u^2 \left(b + (1 - m)u\right)^2} \left|\nabla u\right|^2 d\mathbf{x} - d\_2 \eta \int\_{\Omega} \frac{\mathcal{W}\_\*}{\mathcal{W}^2} \left|\nabla \boldsymbol{\nu}\right|^2 d\mathbf{x} \\ \leq & -\mathcal{K}\_3 \int\_{\Omega} \left(\left|\nabla \boldsymbol{\nu}\right|^2 + \left|\nabla \boldsymbol{\nu}\right|^2\right) d\mathbf{x} := -\sigma(t) \end{split} \tag{40}$$

From Equation (34), *dt d*Ris bounded, where:

$$\frac{d\sigma(t)}{dt} = \mathbf{K}\_3 \int\_{\Omega} 2 \sum\_{i=1}^{n} \frac{\hat{\mathcal{O}}}{\hat{\mathcal{O}}t} (\frac{\hat{\mathcal{O}}u}{\hat{\mathcal{O}}\mathbf{x}\_i}) d\mathbf{x} = \mathbf{K}\_3 \int\_{\Omega} 2 \sum\_{i=1}^{n} \frac{\hat{\mathcal{O}}}{\hat{\mathcal{O}}\mathbf{x}\_i} (\frac{\hat{\mathcal{O}}u}{\hat{\mathcal{O}}t}) d\mathbf{x}\_i$$

By Lemma 1, there is (lim 0) <sup>2</sup> <sup>2</sup> K O-O %" : *dxwu <sup>t</sup>* .

Using the Poincare inequality, there exists:

$$\lim\_{t \to \infty} \int\_{\Omega} (\boldsymbol{w} - \overline{\boldsymbol{w}})^2 d\mathbf{x} = \lim\_{t \to \infty} \int\_{\Omega} (\boldsymbol{u} - \overline{\boldsymbol{u}})^2 d\mathbf{x} = 0 \tag{41}$$

where : K: *<sup>w</sup> wdx* <sup>1</sup> , : K: *udxu* <sup>1</sup> .

In fact, there exist:

$$\left|\Omega\right|(\overline{\mu}-\mu\_\*)^2 = \int\_{\Omega} (\overline{\mu}-\mu\_\*)^2 \,d\mathbf{x} = \int\_{\Omega} (\overline{\mu}-\mu+\mu-\mu\_\*)^2 \,d\mathbf{x} \le \int\_{\Omega} (\overline{\mu}-\mu)^2 \,d\mathbf{x} + \int\_{\Omega} (\mu-\mu\_\*)^2 \,d\mathbf{x}$$

and

$$\left|\Omega\right|(\overline{w}-\boldsymbol{w}\_{\ast})^{2} = \int\_{\Omega} \left(\overline{w}-\boldsymbol{w}\_{\ast}\right)^{2}d\mathbf{x} = \int\_{\Omega} \left(\overline{w}-\boldsymbol{w}+\boldsymbol{w}-\boldsymbol{w}\_{\ast}\right)^{2}d\mathbf{x} \leq \int\_{\Omega} \left(\overline{w}-\boldsymbol{w}\right)^{2}d\mathbf{x} + \int\_{\Omega} \left(\boldsymbol{w}-\boldsymbol{w}\_{\ast}\right)^{2}d\mathbf{x}$$

From Equations (39) and (41), it result in:

$$
\overline{u} \xrightarrow{} u\_\*, \quad \overline{w} \xrightarrow{} w\_\* \text{, as } \ t \to +\infty \tag{42}
$$

From Inequality (34), there exists a subsequence of }{ *mt* which is also denoted }{ *mt* , and nonnegative functions )( <sup>~</sup> *xu* , )( <sup>~</sup> *xw* such that 0)( <sup>~</sup> ),(lim)( <sup>~</sup> ),(lim )( )( <sup>2</sup> <sup>2</sup> *<sup>m</sup>* %" *<sup>m</sup> <sup>c</sup>* : *<sup>m</sup>* %" *<sup>m</sup> <sup>c</sup>* : *xuxtu xwxtw* .

Combined with (42), we obtain:

$$\lim\_{m \to \infty} \left\| \mu(t\_m, \mathbf{x}) - \mu\_\* \right\|\_{c^2(\overline{\Omega})} = \lim\_{m \to \infty} \left\| \nu(t\_m, \mathbf{x}) - \nu\_\* \right\|\_{c^2(\overline{\Omega})} = 0 \tag{43}$$

This above result and the local stability conditions can yield that *E*2 is globally asymptotically stable.

#### *3.2. 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 (3):

$$\begin{cases} u(\mathbf{x},t) = u\_\* + \varepsilon \exp((\mathbf{k}\mathbf{x})i + \lambda t) \\ w(\mathbf{x},t) = w\_\* + \tau \exp((\mathbf{k}\mathbf{x})i + \lambda t) \end{cases} \tag{44}$$

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:

$$
\mathcal{X}^2 - \text{tr}\_k \mathcal{X} + \Delta\_k = 0 \tag{45}
$$

where:

$$\text{tr}\_k = A + R - (d\_1 + d\_2)\mathbf{k}^2, \quad \Delta\_k = AR - B\underline{Q} - \mathbf{k}^2(Ad\_2 + Rd\_1) + (\mathbf{k}^2)^2 d\_1 d\_2 \tag{46}$$

From Equation (45), the dispersion relation of Equation (3) is:

$$\mathcal{A}\_{1,2}(\mathbf{k}) = \frac{tr\_k \pm \sqrt{\left(tr\_k\right)^2 - 4\Delta\_k}}{2} \tag{47}$$

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 *RAtr* @- 0 , 0 *BQAR*  0 . It is clear that 0 *trtrk* @ . Then the stability of the homogeneous state simultaneously changes the sign of *<sup>k</sup>* . From Equation (46), it easily finds that there is @ 0 *<sup>k</sup>* for <sup>2</sup> 2 2 2 <sup>1</sup> **k** @@ TT,where:

$$\left|\kappa\_{\rm l,2}\right|^2 = \frac{Ad\_2 + Rd\_1 \mp \sqrt{\left(Ad\_2 + Rd\_1\right)^2 - 4d\_1d\_2\left(AR - BQ\right)}}{2d\_1d\_2} \tag{48}$$

If <sup>2</sup> T 2,1 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 the other parameters are fixed as:

$$\{r = 1.7, s = 0.5, a = 0.3, b = 0.2, e = 0.75, h = 0.55, d\_1 = 0.02, d\_2 = 0.25, k = 12\}\tag{49}$$

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 .

**Figure 1.** Variation of dispersion relation of Equation (3) around the interior equilibrium point. The red line corresponds to *m* 08.0 , the green is *m* 25.0 and the blue is *m* 35.0 .

#### **4. Turing Pattern Formation**

To better investigate how the prey refuge affects the spatiotemporal dynamics of Equation (3), the spatial distribution diagrams are obtained as change of *m* . All numerical simulations are carried out in a discrete two-dimensional domain with F 200200 lattice sites. The step between each lattice point is defined as U 25.0 . The time evolution of Equation (3) is resorted to the forward Euler integration with a step V 01.0 . The initial value of Equation (3) is placed in the stationary state ),( *wu* \*\* and the perturbation for this value is 0005.0 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

### **196**

snapshots the enclosed color bars denote the range of the changing densities of prey, where higher values correspond to higher prey densities. Figure 1 clearly shows that Equation (3) leads to the Turing instability for *m* @ 2512032.0 . The snapshots for *m* 0, *m* 08.0 , *m* 15.0 and *m* 22.0 are chosen to report the spatial (oscillatory) and temporal ( stationary) dynamics of Equation (3) around the interior equilibrium point, but the snapshots for *m* 27.0 and *m* 35.0 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 . 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. *m*

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 *x* 100 , 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.

**Figure 2.** Spatial distributions of prey obtained with Equation (3) for (**a**) *m* 0 , (**b**) *m* 08.0 , (**c**) *m* 15.0 , (**d**) *m* 22.0 , (**e**) *m* 27.0 , (**f**) *m* 35.0 . Other parameters are fixed as Equation (49).

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

#### **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 No. 11226256), by the Key Program of Zhejiang Provincial Natural Science Foundation of China (Grant No. LZ12C0300), and by the Zhejiang Provincial Natural Science Foundation of China (Grant No. LY13A010010).

### **Conflicts of Interest**

The authors declare no conflict of interest.

## **References**


Reprinted from *Entropy*. Cite as: Zhang, Y.; Xu, H.; Lv, X.; Wu, J. Diffusion in Relatively Homogeneous Sand Columns: A Scale-Dependent or Scale-Independent Process? *Entropy* 2013, *15*, 4376–4391.

## *Article*
