# A Stage-Structured Predator-Prey Model in a Patchy Environment

Authors: Xuejuan Lu; Yuming Chen; Shengqiang Liu
Journal: Complexity
 (2020)

---

## Abstract

In this paper, we propose a stage-structured predator-prey model with migrations among patches in ann-patch environment. The net reproduction number for each patch in isolation is obtained along with the net reproduction number of the system of patches, ℛ0. Inequalities describing the relationship among these numbers are also given. Furthermore, threshold dynamics determined by ℛ0 is established: the predator dies out if ℛ0<1 while the predator persists if ℛ0>1. Focusing on the case with two patches, we obtain that the dispersal decreases the net reproduction number ℛ0. By numerical simulations, we find that the dispersal may be a good thing or a bad thing because the dispersal could make the predator population thrive or extinct, and hence we might seek steady state in the ecological environment by controlling parameters related to the prey and the predator.


---

## Body

## 1. Introduction

The effect of dispersal of organisms on population dynamics is one of the central topics in ecology and evolutionary biology [1–5]. Each species has a characteristic pattern of variation in abundance over space. Understanding the factors responsible for manifest diversity in distributional patterns is an important prerequisite for studying survival and extinction of species. For example, many types of birds and mammals migrate from cold regions to warm regions in search of a better habitat or a breeding site [6], and for ecological communities of insects, dispersal of a predator is usually driven by its nonrandom foraging behavior which can often respond to prey-contact stimuli [7, 8].During the past couple of decades, due to the prevalence of the predator-prey relationship in the nature, predator-prey models with diffusion in patchy environment have attracted significant attention from ecologists, biologists, and biomathematicians. Levin [9] proposed two-species competition and prey-predator models with population dispersal among patches. Takeuchi [10] and Kuang and Takeuchi [11] found that stabilizing and destabilizing effects could be induced by prey dispersal. Gramlich et al. [12] considered a 2-patch 2-species system and found that dispersal between different habitats influences the dynamics and stability of populations considerably. Moreover, these effects depend on the local interactions of a population with others. Kang et al. [8] formulated a two-patch Rosenzweig–MacArthur prey-predator model with mobility only in the predator, and they found that dispersal may stabilize or destabilize the coupled system; dispersal may generate multiple interior equilibria or may destroy interior equilibria. In a few patch models, the common findings for most studies on interactions of prey and predator are that dispersal makes the dynamics very complicated. This means that dispersals of prey and predator play an important role in regulating, stabilizing, or destabilizing population dynamics of both prey and predator. In an n-patch environment, due to the intricacies that arise from dispersal, not much has been done yet. Qiu and Mitsui [13] considered a predator-prey system with diffusion and time delay in an n-patch environment, where prey disperses between n patches of a heterogeneous environment with barriers between patches where a predator cannot cross. Al-Darabsah et al. [14] proposed a prey-predator model in multiple patches through the stage-structured maturation time delay with migrations among patches. However, they only discussed the existence of equilibrium points and the uniform persistence for the special case of two patches.To our best acknowledgment, little attention has been paid to the relation between the long-term dynamics and dispersal rate of the multipatch predator-prey model. In two recent references [15, 16], net reproduction number ℛ0 was introduced into the nondiffusion predator-prey models; it was shown that the net reproduction number determines its threshold dynamics; here ℛ0 has the similar definition to the basic reproduction number of general compartmental epidemic model, as discussed in [17]. Since there are many excellent results on investigating the effects of mobility rate on basic reproduction number and the threshold dynamics of multipatch epidemic model, for reference, we refer to Hsieh et al. [18], Wang et al. [19], and Wang et al. [20]. This motivates us to study an n-patch predator-prey model with prey and predator diffusing simultaneously.In this paper, we focus on a model with prey and predator dispersal in ann-patch environment. In addition, the predator is structured as juvenile and adult. We assume that juvenile predators may diffuse as they can follow adult one. We also need to pay attention to acquisition probability of the juvenile predator, which can be zero and less than unity. Moreover, the dispersal among patches may occur only in a single species. The main focus of our study of such prey-predator interactions in an n-patch heterogeneous environment is to explore the following ecological questions:(1)
What is the net reproduction numberℛ0? And whether there is a limit on ℛ0?(2)
In a two-patch model, what is the relationship between the constant dispersal rates of prey and predator andℛ0?(3)
How does dispersal of prey and predator affect the extinction and persistence of predator in all patches?(4)
How may dispersal promote the coexistence of prey and predator when predator goes extinct in a single patch?The paper is organized as follows. In the next section, we propose the model and show the positiveness and ultimate boundedness of solution under some conditions. In Section3, we derive a general formula for the net reproduction number and give a bound for ℛ0, show the local asymptotical stability of the predator-free equilibrium when ℛ0<1, and also obtain the global attractivity of the predator-free equilibrium under some conditions. Moreover, we consider the special case with the linear predation function in a two-patch environment. In Section 4, we prove the uniform persistence when ℛ0>1, and hence the system has at least one positive equilibrium. Numerical simulations are provided in Section 5 to demonstrate our results. The paper concludes with a brief discussion.


## 2. The Model

In this section, we formulate a predator-prey model inn patches by taking into consideration diffusion among patches and a stage structure in the predator. In each patch i, there are prey individuals, juvenile predator individuals, and adult predator individuals, denoted by Ni, PiJ, and PiA, respectively.For patchi, let ri be the intrinsic growth rate in the absence of predation and Ki be the carrying capacity of the prey; diN, diJ, and diA denote the death rates of prey, juvenile predators, and adult predators, respectively; we always assume ri>diN; σi stands for the acquisition probability of juvenile predators and σi∈0,1; we always assume ei to be the probability of the captured prey being converted into juvenile predators; αi is the progression rate of juvenile predators. All parameters are assumed to be positive except that σi can be zero.We assume that the predation functionfiNi satisfies the following basic assumptions for Ni∈0,∞.(1)
fiNi>0, i∈ℕn=1,2,…,n.(2)
fiNi=NigiNi, giNi is a continuous function about Ni and giNi>0, i∈ℕn.(3)
fi′Ni≥0, i∈ℕn.The following three types of predation functionsfN in [21] satisfy the above assumptions.(1)
fN=bN, b>0.(2)
fN=bN/1+cN, b,c>0.(3)
fN=bN2/1+cN2, b,c>0.We assume that there is no capture and captured behavior in the diffusion process, andmijK for K=N,J,A are the constant dispersal rates from patch j to patch i for i≠j of prey, juvenile predators, and adult predators, respectively, with miiK=0. It is assumed that the matrices MK=mijK, i,j∈ℕn for K=N,J,A are irreducible. The flowchart of diffusion is shown in Figure 1.Figure 1
Flowchart of the dispersal process between patchi and j.Based on the above assumptions, our model is as follows:(1)N˙i=riNi1−NiKi−fiNiPiA+σiPiJ−diNNi+∑j∈ℕnmijNNj−∑j∈ℕnmjiNNi,P˙iJ=eifiNiPiA+σiPiJ−diJ+αiPiJ+∑j∈ℕnmijJPjJ−∑j∈ℕnmjiJPiJ,P˙iA=αiPiJ−diAPiA+∑j∈ℕnmijAPjA−∑j∈ℕnmjiAPiA.Similar to [19], for a vector x∈ℝn, we use diagx to denote the n×n diagonal matrix, whose diagonal elements are the components of x. We use the ordering in ℝn generated by the cone ℝ+n, that is, x≤y if y−x∈ℝ+n, x<y if x≤y and x≠y, and finally x≪y means xi<yi for any index i.Denote the vector(2)Nt,PJt,PAt=N1t,…,Nnt,P1Jt,…,PnJt,P1At,…,PnAtT∈ℝ+3n,(3)A0=N0,PJ0,PA0.Moreover, we assume that each component ofA0 is nonnegative with the following initial conditions:(4)∑i∈ℕnNi0>0,∑i∈ℕnPiA0+σiPiJ0>0.Then, it follows from the standard existence and uniqueness theorem for ordinary differential equations that there is a unique solution to system (1).Theorem 1.
Consider system (1) with nonnegative initial conditions (4). Then, for each i∈ℕn and t>0, Nt, PJt, and PAt are all positive and ultimately bounded.Proof of Theorem1.
Firstly, we prove that fori∈ℕn, Nit>0 for t>0. For convenience, we rewrite the first equation of (1) as(5)N˙it=PitNit+Qit,i∈ℕn,where(6)Pit=ri−riKiNit−giNitPiAt+σiPiJt−diN−∑j∈ℕnmjiN,Qit=∑j∈ℕnmijNNjt.
LetI1=i∈ℕnNi0>0 and I2=ℕn\I1. Since ∑i∈ℕnNi0>0, we know that I1≠∅. If i0∈I1, then from (5) with i=i0, we easily see that Ni0t>0 for t>0. If I2=∅, then we are done. Now suppose that I2≠∅. Since mijN is irreducible, there exists i1∈I2 such that mi1jN≠0 for some j∈I1. Then, Qi1t>0 for t>0. This and (5) with i=i1 imply that Ni1t>0 for t>0. Denote J1=I1∪i1 and J2=I2\i1. If J2=∅, then we are done. Otherwise, continuing this way, after a finite number of steps, we have Nit>0 for t>0 and i∈ℕ.
Next, we prove that fori∈ℕn, PiJt>0 and PiAt>0 for t>0. Since ∑i∈ℕnPiA0+σiPiJ0>0, there exists i2∈ℕn such that Pi2A0+σi2Pi2J0>0. First, we assume that Pi2A0>0. By noting(7)P˙iAt=−diA−∑j∈ℕnmjiAPiAt+αiPiJt+∑j∈ℕnmijAPjAt,which is similar to (5), we can use similar arguments as those of showing Nit>0 for t>0 and i∈ℕn to obtain PiAt>0 for t>0 and i∈ℕn. Then, from(8)P˙iJt≥−diJ+αi+∑j∈ℕnmjiJPiJt+eifiNitPiAt,we see that PiJt>0 for t>0 and i∈ℕn. Now, assume that σi2Pi2J0>0, which implies that Pi2J0>0. Then, from the second equation of (1), we get Pi2Jt>0 for t>0. This, combined with the third equation of (1), yields Pi2At>0 for t>0. Then as before, we can show that PiAt>0 and PiJt>0 for t>0 and i∈ℕn.
Finally, we show the boundedness of solutions. On the one hand, letNt=∑i∈ℕnNit. Then, by the Cauchy–Schwartz inequality, we have(9)N˙t≤∑i∈ℕnri−diNNit−∑i∈ℕnriKiNi2t≤maxi∈ℕnri−diNNt−mini∈ℕnriKi∑i∈ℕnNi2t≤maxi∈ℕnri−diNNt−mini∈ℕnriKiN2tn.
It follows thatlimsupt⟶∞Nt≤nmaxi∈ℕnri−diN/mini∈ℕnri/Ki. On the other hand, let Gt=∑i∈ℕneiNit+PiJt+PiAt. Then,(10)G˙t=∑i∈ℕneiriNit1−NitKi−eidiNNit+ei∑j∈ℕnmijNNjt−ei∑j∈ℕnmjiNNit−diJPiJt−diAPiAt≤∑i∈ℕn2eiriNit+eimaxj∈ℕnmijNNt−∑i∈ℕneiriNit+diJPiJt+diAPiAt≤maxi∈ℕn2eiri+∑i∈ℕneimaxj∈ℕnmijNNt−mini∈ℕnei,diJ,diAGt.
It is easy to see from this andlimsupt⟶∞Nt≤nmaxi∈ℕnri−diN/mini∈ℕnri/Ki that(11)limsupt⟶∞Gt≤maxi∈ℕn2eiri+∑i∈ℕneimaxj∈ℕnmijNmini∈ℕnei,diJ,diA×nmaxi∈ℕnri−diNmini∈ℕnri/Ki.
This completes the proof.
To show the existence of predator-free equilibria (PFE), we letPiJ=PiA=0,i∈ℕn in equation (1) to get(12)N˙i=riNi1−NiKi−diNNi+∑j∈ℕnmijNNj−∑j∈ℕnmjiNNi,∑i∈ℕnNi0>0.
In order to obtain a positive equilibrium for system (12), we assume that(13)sdiagri−diN+MN>0,where s⋅ represents the stability modulus of an n×n matrix and is defined by sH:=maxReλλ is an eigenvalue of H for an n×n matrix H. From [19, 20, 22], we have the following two lemmas.Lemma 1.
Under assumption (13), system (12) admits a unique equilibrium N0=N10,N20,…,Nn0T, which is positive and globally asymptotically stable in ℝ+n\0.Lemma 2.
Under assumption (13), system (1) admits a unique PFE P0=N0,0,0T, where N0=N10,N20,…,Nn0T∈ℝ+n is the unique equilibrium of (12) and 0=0,…,0︸nT∈ℝ+n.


## 3. Threshold Dynamics

### 3.1. The Net Reproduction Number

In this section, we first give definition of the net reproduction numberℛ0 for the general predator-prey dispersal model (1), which is similar to the basic reproduction number for epidemic system; for reference, we refer to [17, 23].LetP0 be the PFE of (1). From [17], we have(14)F=F11F1200=diageiσifiNi0diageifiNi000,(15)V=V110−V21V22,where(16)V11=diagdiJ+αi+∑j∈ℕnmjiJ−MJ,V21=diagαi,V22=diagdiA+∑j∈ℕnmjiA−MA.Note thatV11 and V22 are both irreducible nonsingular M-matrices with positive column sums and hence V11−1>0 and V22−1>0. Following the procedure introduced in [17], the net reproduction number ℛ0 of equation (1) is given by ℛ0=ρFV−1, where ρ represents the spectral radius of the matrix. Let M1=F−V. Then, by [17], there hold two equivalences:(17)ℛ0>1⟺sM1>0,ℛ0<1⟺sM1<0.As(18)FV−1=F11F1200V11−10V22−1V21V11−1V22−1,we have(19)ℛ0=ρF11V11−1+F12V22−1V21V11−1.In equation (19), F11V11−1 accounts for the contribution of juvenile predators to the adult predators, while F12V22−1V21V11−1 accounts for the contribution by adult predators themselves. Diffusion rates of preys influence the predator population capture rate, and diffusion rates of predators affect their consumption.The net reproduction number gives an important threshold of then-patch diffusion predator-prey system.Theorem 2.
Consider model (1). Suppose (13) holds. If ℛ0<1, then the PFE is locally asymptotically stable, and if ℛ0>1, then the PFE is unstable. Moreover, if ℛ0<1, then the PFE is globally attractive for solutions with initial conditions (2).Proof of Theorem2.
It follows from ([10], Theorem 2) that the PFE is locally asymptotically stable if ℛ0<1 and unstable if ℛ0>1.
Now, we prove the global attractivity of the PFE as follows.
LetNt,PJt,PAt be an arbitrary solution of (1). We first show that limt⟶∞PJt=limt⟶∞PAt=0.
SincesM1<0, there exists an ε>0 such that sM1+M2ε<0, where(20)M2ε=diageiσifiNi0+ε−fiNi0diageifiNi0+ε−fiNi000.
Note thatNt satisfies(21)N˙it≤riNit1−NitKi−diNNit+∑j∈ℕnmijNNjt−∑j∈ℕnmjiNNt,where i∈ℕn. By Lemma 1 and Comparison Theorem [24, 25], we have limsupt⟶∞Nt≤N0. It follows that there exists t0≥0 such that Nt≤N0+ε,…,ε︸n for t≥t0. Then, for t≥t0, PJt,PAt satisfies(22)P˙JtP˙At≤M1+M2εPJtPAt.
SincesM1+M2ε<0, we obtain limt⟶∞PJt=limt⟶∞PAt=0 by applying Comparison Theorem [24, 25].
Next, we prove thatlimt⟶∞Nt=N0. It suffices to show liminft⟶∞Nt≥N0 as limsupt⟶∞Nt≤N0. Since gi is continuous and Ni is bounded, i∈ℕn, there exists κ>0 such that giNit≤κ for t≥0. By equation (13), there exists η0>0 such that sdiagri−diN−1+σiκη+MN>0 for η∈0,η0. For any η∈0,η0, there exists t1≥0 such that 0≤PiJt,PiAt≤η for t≥t1. Then,(23)N˙it≥riNit1−NitKi−diN+1+σiηκNit+∑j∈ℕnmijNNjt−∑j∈ℕnmjiNNitfor t≥t1 and i∈ℕn. By Lemma 1 and Comparison Theorem [24, 25], we know that liminft⟶∞Nt≥N0η, where N0η is the positive equilibrium of (12) with di being replaced by diN+1+σiηκ, i∈ℕn. It is easy to see that limη⟶0N0η=N0. Therefore, we have liminft⟶∞Nt≥N0 as required.
To summarize, we have shown thatlimt⟶∞Nt,PJt,PAt=N0,0,0. This completes the proof.
Next, we consider the effect of dispersal onℛ0 and provide lower and upper bounds of ℛ0. In the rest of this paper, we let ai=diJ+αi for i∈ℕn. When there is no diffusion of prey, that is, MN=0n×n, we have N0=N10,N20,…,Nn0T with Ni0=Kri−diN/ri for i∈ℕn and the net reproduction number in patch i in isolation is given by(24)ℛ0i=eiσifiKiri−diN/riai+eiαifiKiri−diN/ridiAai.
We define a modified reproduction number in patchi:(25)ℛ˜0i=eiσifiKiri−diN/riai+∑j=1nmjiJ+eiαifiKiri−diN/ridiA+∑j=1nmjiAai+∑j=1nmjiJ.
Clearly,ℛ˜0i<ℛ0i. Then, we can obtain the following result on the bounds for the net reproduction number for system (1), in terms of the numbers defined in equations (24) and (25) for each patch.Theorem 3.
For system (1), ifMN=0n×n, then(26)ℛ0≥max1≤i≤nℛ˜0i.
Moreover, ifai=a, ei=e, αi=α, σi=σ, fi=f, diN=d1, diA=d2, and ri=r for i∈ℕn, then(27)maxmax1≤i≤nℛ˜0i,min1≤i≤nℛ0i≤ℛ0≤max1≤i≤nℛ0i.Proof of Theorem3.
From [18, 26], let Vjj1′ denote the matrix obtained by deleting the first row and the first column of Vjj (j=1,2); Y=yij and Z=zij denote V11−1 and V22−1, respectively. Let W=wij=G+H, where G=gij=F11Y and H=hij=F12ZV21Y. It follows that yij>0 and zij>0 for i,j∈ℕn. By [26], we have(28)ℛ0=ρW≥wii=gii+hii,i∈ℕn.
Here,(29)g11=e1σ1f1K1r1−d1Nr1y11=e1σ1f1K1r1−d1Nr1detV111′detV11.
By virtue of Fisher’s inequality [27],(30)detV11≤a1+∑j=1nmj1JdetV111′.
Then,(31)g11≥e1σ1f1K1r1−d1N/r1a1+∑j=1nmj1J.
Similarly,(32)h11=e1α1f1K1r1−d1Nr1z11y11+∑k=2ne1f1K1r1−d1Nr1αkz1kyk1.
Thus, we have(33)h11≥e1α1f1K1r1−d1Nr1z11y11=e1α1f1K1r1−d1Nr1detV221′detV22detV111′detV11.
Again, by Fisher’s inequality,(34)h11≥e1α1f1K1r1−d1N/r1d1A+∑j=1nmj1Aa1+∑j=1nmj1J.
Then, from these conclusions, we know(35)ℛ0≥ℛ˜01.
Similarly, it can be shown that(36)ℛ0≥ℛ˜0i,i=2,3,…,n,and (26) is proved.
If we letai=a, ei=e, αi=α, σi=σ, fi=f, diN=d1, diA=d2, and ri=r for i∈ℕn, then wij=eσfKir−d1/ryij+eαfKir−d1/r∑k=1nzikyik for i,j∈ℕn. Without loss of generality, assume that 0<K1≤K2≤…≤Kn. From the fact that the sum of each column of the matrix V11 is equal to a>0 and the sum of each column of the matrix V22 is d2>0, we see that ∑i=1nyij=1/a and ∑i=1nzij=1/d2 for j∈ℕn. Thus, for the matrix W, the sum of column j is(37)∑i=1nwij=∑i=1neσfKir−d1ryij+∑i=1neαfKir−d1r∑k=1nzikyik≤eσfKnr−d1r∑i=1nyij+eαfKnr−d1r∑i=1n∑k=1nzikyik=eσfKnr−d1r∑i=1nyij+eαfKnr−d1r∑k=1n∑i=1nzikyik=eσfKnr−d1/ra+eαfKnr−d1/rd2a=R0n.
Similarly, we can prove∑i=1nwij≥eσfK1r−d1/r/a+eαfK1r−d1/r/d2a=ℛ01. From [26], we obtain the following inequalities:(38)min1≤i≤nℛ0i≤ℛ0≤max1≤i≤nℛ0i.
Combining this with (26), we get (27). This completes the proof.


### 3.2. The Model with Two Patches

In this section, we consider a special case of equation (1) with n=2 and assume that fiNi=biNi, bi>0 for i=1,2, namely,(39)N˙1=r1N11−N1K1−b1N1P1A+σ1P1J−d1NN1+m12NN2−m21NN1,P˙1J=e1b1N1P1A+e1b1σ1N1P1J−d1J+α1P1J+m12JP2J−m21JP1J,P˙1A=α1P1J−d1AP1A+m12AP2A−m21AP1A,N˙2=r2N21−N2K2−b2N2P2A+σ2P2J−d2NN2+m21NN1−m12NN2,P˙2J=e2b2N2P2A+e2b2σ2N2P2J−d2J+α2P2J+m21JP1J−m12JP2J,P˙2A=α2P2J−d2AP2A+m21AP1A−m12AP2A.The PFE of equation (39) is N10,N20,0,0,0,0T, where N10,N20T satisfies(40)r1N11−N1K1−d1NN1+m12NN2−m21NN1=0,r2N21−N2K2−d2NN2+m21NN1−m12NN2=0.Let(41)N0=N10,N20T.Note that(42)F11=diage1σ1b1N10,e2σ2b2N20,F12=diage1b1N10,e2b2N20,(43)V21=diagα1,α2,(44)V11=a1+m21J−m12J−m21Ja2+m12J,V22=d1A+m21A−m12A−m21Ad2A+m12A.Recall thatai=diJ+αi for i=1,2. Hence, we have(45)ℛ0=ρF11V11−1+F12V22−1V21V11−1=ρW,where(46)W=w11w12w21w22,with(47)w11=e1σ1b1N10a2+m12Ja1a2+a1m12J+a2m21J+e1b1α1N10a2+m12Jd2A+m12A+e1b1α2N10m21Jm12Aa1a2+a1m12J+a2m21Jd1Ad2A+d1Am12A+d2Am21A,w12=e1σ1b1N10m12Ja1a2+a1m12J+a2m21J+e1b1α1N10m12Jd2A+m12A+e1b1α2N10m12Aa1+m21Ja1a2+a1m12J+a2m21Jd1Ad2A+d1Am12A+d2Am21A,w21=e2σ2b2N20m21Ja1a2+a1m12J+a2m21J+e2b2α1N20m21Aa2+m12J+e2b2α2N20m21Jd1A+m21Aa1a2+a1m12J+a2m21Jd1Ad2A+d1Am12A+d2Am21A,w22=e2σ2b2N20a1+m21Ja1a2+a1m12J+a2m21J+e2b2α1N20m12Jm12A+e2b2α2N20a1+m21Jd1A+m21Aa1a2+a1m12J+a2m21Jd1Ad2A+d1Am12A+d2Am21A.It follows from equation (45) that(48)ℛ0=12w11+w22+w11−w222+4w12w21.Theorem 4.
Ife1=e2=e, σ1=σ2=σ, a1=a2=a, b1=b2=b, α1=α2=α, d1A=d2A=d, mijJ=m, and mijA=km, i,j=1,2, where k is a positive constant, then ℛ0 decreases as m increases.Proof of Theorem 4.
From equation (46) and under the assumptions on parameters, we get(49)W=w11w12w21w22=b1N10ub1N10vb2N20vb2N20u,where(50)u=eσa+maa+2m+eα2km2+d+kam+adada+2md+2km,v=eσmaa+2m+eαm2km+d+akada+2md+2km.
By equation (48), ℛ0 is the larger root of the following quadratic equation:(51)λ2−b1N10+b2N20uλ+b1b2N10N20u2−v2=0.
Note that(52)u+v=eσd+αad,u−v=eα+σd+2kma+2md+2km,(53)u2−v2=eσd+αad⋅eα+σd+2kma+2md+2km.
It follows fromu+v=eσd+α/ad that ∂u/∂m=−∂v/∂m. Taking partial derivatives with respect to m for equation (51) gives(54)2λ−b1N10+b2N20u∂λ∂m=λb1N10+b2N20∂u∂m−b1b2N10N20∂u2−v2∂m=λb1N10+b2N20−2b1b2N10N20eα+σdad∂u∂m.
Moreover, from the definition ofu, we have(55)∂u∂m=−eσaaa+2m2−eαad⋅6kmad+8k2m3+ad2+4km2d+ka2d+4k2m2aa+2m2d+2km2<0.
Note that2λλ=ℛ0>b1N10+b2N20u. Thus, to show ∂ℛ0/∂m<0, that is to say, ∂λ/∂mλ=ℛ0<0, it is sufficient to show for λ=ℛ0,(56)b1N10+b2N20λ−2b1b2N10N20⋅eσd+αad>0.
We claim thatλ>1/2b1N10+b2N20u+v. Since(57)λ=12b1N10+b2N20u+b1N10+b2N202u2−4b1b2N10N20u2−v2,it is equivalent to show that(58)b1N10+b2N202u2−4b1b2N10N20u2−v2>b1N10+b2N202v2,or(59)b1N10−b2N202u2−v2>0.
This is automatically true sinceu>v>0. It follows from λ>1/2b1N10+b2N20u+v that(60)b1N10+b2N20λ>12b1N10+b2N202u+v=12eσd+αadb1N10+b2N202.
Thus, as required,(61)b1N10+b2N20λ−2b1b2N10N20eσd+αad>12eσd+αadb1N10−b2N202>0.
This completes the proof.Remark 1.
Under the condition of Theorem4, the following result can be proved similarly; each entry of the matrix W decreases in terms of a or d, and it increases in terms of e, σ, and α. Hence, the net reproduction number decreases in terms of a or d, and it increases in terms of e, σ, and α.


## 3.1. The Net Reproduction Number

In this section, we first give definition of the net reproduction numberℛ0 for the general predator-prey dispersal model (1), which is similar to the basic reproduction number for epidemic system; for reference, we refer to [17, 23].LetP0 be the PFE of (1). From [17], we have(14)F=F11F1200=diageiσifiNi0diageifiNi000,(15)V=V110−V21V22,where(16)V11=diagdiJ+αi+∑j∈ℕnmjiJ−MJ,V21=diagαi,V22=diagdiA+∑j∈ℕnmjiA−MA.Note thatV11 and V22 are both irreducible nonsingular M-matrices with positive column sums and hence V11−1>0 and V22−1>0. Following the procedure introduced in [17], the net reproduction number ℛ0 of equation (1) is given by ℛ0=ρFV−1, where ρ represents the spectral radius of the matrix. Let M1=F−V. Then, by [17], there hold two equivalences:(17)ℛ0>1⟺sM1>0,ℛ0<1⟺sM1<0.As(18)FV−1=F11F1200V11−10V22−1V21V11−1V22−1,we have(19)ℛ0=ρF11V11−1+F12V22−1V21V11−1.In equation (19), F11V11−1 accounts for the contribution of juvenile predators to the adult predators, while F12V22−1V21V11−1 accounts for the contribution by adult predators themselves. Diffusion rates of preys influence the predator population capture rate, and diffusion rates of predators affect their consumption.The net reproduction number gives an important threshold of then-patch diffusion predator-prey system.Theorem 2.
Consider model (1). Suppose (13) holds. If ℛ0<1, then the PFE is locally asymptotically stable, and if ℛ0>1, then the PFE is unstable. Moreover, if ℛ0<1, then the PFE is globally attractive for solutions with initial conditions (2).Proof of Theorem2.
It follows from ([10], Theorem 2) that the PFE is locally asymptotically stable if ℛ0<1 and unstable if ℛ0>1.
Now, we prove the global attractivity of the PFE as follows.
LetNt,PJt,PAt be an arbitrary solution of (1). We first show that limt⟶∞PJt=limt⟶∞PAt=0.
SincesM1<0, there exists an ε>0 such that sM1+M2ε<0, where(20)M2ε=diageiσifiNi0+ε−fiNi0diageifiNi0+ε−fiNi000.
Note thatNt satisfies(21)N˙it≤riNit1−NitKi−diNNit+∑j∈ℕnmijNNjt−∑j∈ℕnmjiNNt,where i∈ℕn. By Lemma 1 and Comparison Theorem [24, 25], we have limsupt⟶∞Nt≤N0. It follows that there exists t0≥0 such that Nt≤N0+ε,…,ε︸n for t≥t0. Then, for t≥t0, PJt,PAt satisfies(22)P˙JtP˙At≤M1+M2εPJtPAt.
SincesM1+M2ε<0, we obtain limt⟶∞PJt=limt⟶∞PAt=0 by applying Comparison Theorem [24, 25].
Next, we prove thatlimt⟶∞Nt=N0. It suffices to show liminft⟶∞Nt≥N0 as limsupt⟶∞Nt≤N0. Since gi is continuous and Ni is bounded, i∈ℕn, there exists κ>0 such that giNit≤κ for t≥0. By equation (13), there exists η0>0 such that sdiagri−diN−1+σiκη+MN>0 for η∈0,η0. For any η∈0,η0, there exists t1≥0 such that 0≤PiJt,PiAt≤η for t≥t1. Then,(23)N˙it≥riNit1−NitKi−diN+1+σiηκNit+∑j∈ℕnmijNNjt−∑j∈ℕnmjiNNitfor t≥t1 and i∈ℕn. By Lemma 1 and Comparison Theorem [24, 25], we know that liminft⟶∞Nt≥N0η, where N0η is the positive equilibrium of (12) with di being replaced by diN+1+σiηκ, i∈ℕn. It is easy to see that limη⟶0N0η=N0. Therefore, we have liminft⟶∞Nt≥N0 as required.
To summarize, we have shown thatlimt⟶∞Nt,PJt,PAt=N0,0,0. This completes the proof.
Next, we consider the effect of dispersal onℛ0 and provide lower and upper bounds of ℛ0. In the rest of this paper, we let ai=diJ+αi for i∈ℕn. When there is no diffusion of prey, that is, MN=0n×n, we have N0=N10,N20,…,Nn0T with Ni0=Kri−diN/ri for i∈ℕn and the net reproduction number in patch i in isolation is given by(24)ℛ0i=eiσifiKiri−diN/riai+eiαifiKiri−diN/ridiAai.
We define a modified reproduction number in patchi:(25)ℛ˜0i=eiσifiKiri−diN/riai+∑j=1nmjiJ+eiαifiKiri−diN/ridiA+∑j=1nmjiAai+∑j=1nmjiJ.
Clearly,ℛ˜0i<ℛ0i. Then, we can obtain the following result on the bounds for the net reproduction number for system (1), in terms of the numbers defined in equations (24) and (25) for each patch.Theorem 3.
For system (1), ifMN=0n×n, then(26)ℛ0≥max1≤i≤nℛ˜0i.
Moreover, ifai=a, ei=e, αi=α, σi=σ, fi=f, diN=d1, diA=d2, and ri=r for i∈ℕn, then(27)maxmax1≤i≤nℛ˜0i,min1≤i≤nℛ0i≤ℛ0≤max1≤i≤nℛ0i.Proof of Theorem3.
From [18, 26], let Vjj1′ denote the matrix obtained by deleting the first row and the first column of Vjj (j=1,2); Y=yij and Z=zij denote V11−1 and V22−1, respectively. Let W=wij=G+H, where G=gij=F11Y and H=hij=F12ZV21Y. It follows that yij>0 and zij>0 for i,j∈ℕn. By [26], we have(28)ℛ0=ρW≥wii=gii+hii,i∈ℕn.
Here,(29)g11=e1σ1f1K1r1−d1Nr1y11=e1σ1f1K1r1−d1Nr1detV111′detV11.
By virtue of Fisher’s inequality [27],(30)detV11≤a1+∑j=1nmj1JdetV111′.
Then,(31)g11≥e1σ1f1K1r1−d1N/r1a1+∑j=1nmj1J.
Similarly,(32)h11=e1α1f1K1r1−d1Nr1z11y11+∑k=2ne1f1K1r1−d1Nr1αkz1kyk1.
Thus, we have(33)h11≥e1α1f1K1r1−d1Nr1z11y11=e1α1f1K1r1−d1Nr1detV221′detV22detV111′detV11.
Again, by Fisher’s inequality,(34)h11≥e1α1f1K1r1−d1N/r1d1A+∑j=1nmj1Aa1+∑j=1nmj1J.
Then, from these conclusions, we know(35)ℛ0≥ℛ˜01.
Similarly, it can be shown that(36)ℛ0≥ℛ˜0i,i=2,3,…,n,and (26) is proved.
If we letai=a, ei=e, αi=α, σi=σ, fi=f, diN=d1, diA=d2, and ri=r for i∈ℕn, then wij=eσfKir−d1/ryij+eαfKir−d1/r∑k=1nzikyik for i,j∈ℕn. Without loss of generality, assume that 0<K1≤K2≤…≤Kn. From the fact that the sum of each column of the matrix V11 is equal to a>0 and the sum of each column of the matrix V22 is d2>0, we see that ∑i=1nyij=1/a and ∑i=1nzij=1/d2 for j∈ℕn. Thus, for the matrix W, the sum of column j is(37)∑i=1nwij=∑i=1neσfKir−d1ryij+∑i=1neαfKir−d1r∑k=1nzikyik≤eσfKnr−d1r∑i=1nyij+eαfKnr−d1r∑i=1n∑k=1nzikyik=eσfKnr−d1r∑i=1nyij+eαfKnr−d1r∑k=1n∑i=1nzikyik=eσfKnr−d1/ra+eαfKnr−d1/rd2a=R0n.
Similarly, we can prove∑i=1nwij≥eσfK1r−d1/r/a+eαfK1r−d1/r/d2a=ℛ01. From [26], we obtain the following inequalities:(38)min1≤i≤nℛ0i≤ℛ0≤max1≤i≤nℛ0i.
Combining this with (26), we get (27). This completes the proof.


## 3.2. The Model with Two Patches

In this section, we consider a special case of equation (1) with n=2 and assume that fiNi=biNi, bi>0 for i=1,2, namely,(39)N˙1=r1N11−N1K1−b1N1P1A+σ1P1J−d1NN1+m12NN2−m21NN1,P˙1J=e1b1N1P1A+e1b1σ1N1P1J−d1J+α1P1J+m12JP2J−m21JP1J,P˙1A=α1P1J−d1AP1A+m12AP2A−m21AP1A,N˙2=r2N21−N2K2−b2N2P2A+σ2P2J−d2NN2+m21NN1−m12NN2,P˙2J=e2b2N2P2A+e2b2σ2N2P2J−d2J+α2P2J+m21JP1J−m12JP2J,P˙2A=α2P2J−d2AP2A+m21AP1A−m12AP2A.The PFE of equation (39) is N10,N20,0,0,0,0T, where N10,N20T satisfies(40)r1N11−N1K1−d1NN1+m12NN2−m21NN1=0,r2N21−N2K2−d2NN2+m21NN1−m12NN2=0.Let(41)N0=N10,N20T.Note that(42)F11=diage1σ1b1N10,e2σ2b2N20,F12=diage1b1N10,e2b2N20,(43)V21=diagα1,α2,(44)V11=a1+m21J−m12J−m21Ja2+m12J,V22=d1A+m21A−m12A−m21Ad2A+m12A.Recall thatai=diJ+αi for i=1,2. Hence, we have(45)ℛ0=ρF11V11−1+F12V22−1V21V11−1=ρW,where(46)W=w11w12w21w22,with(47)w11=e1σ1b1N10a2+m12Ja1a2+a1m12J+a2m21J+e1b1α1N10a2+m12Jd2A+m12A+e1b1α2N10m21Jm12Aa1a2+a1m12J+a2m21Jd1Ad2A+d1Am12A+d2Am21A,w12=e1σ1b1N10m12Ja1a2+a1m12J+a2m21J+e1b1α1N10m12Jd2A+m12A+e1b1α2N10m12Aa1+m21Ja1a2+a1m12J+a2m21Jd1Ad2A+d1Am12A+d2Am21A,w21=e2σ2b2N20m21Ja1a2+a1m12J+a2m21J+e2b2α1N20m21Aa2+m12J+e2b2α2N20m21Jd1A+m21Aa1a2+a1m12J+a2m21Jd1Ad2A+d1Am12A+d2Am21A,w22=e2σ2b2N20a1+m21Ja1a2+a1m12J+a2m21J+e2b2α1N20m12Jm12A+e2b2α2N20a1+m21Jd1A+m21Aa1a2+a1m12J+a2m21Jd1Ad2A+d1Am12A+d2Am21A.It follows from equation (45) that(48)ℛ0=12w11+w22+w11−w222+4w12w21.Theorem 4.
Ife1=e2=e, σ1=σ2=σ, a1=a2=a, b1=b2=b, α1=α2=α, d1A=d2A=d, mijJ=m, and mijA=km, i,j=1,2, where k is a positive constant, then ℛ0 decreases as m increases.Proof of Theorem 4.
From equation (46) and under the assumptions on parameters, we get(49)W=w11w12w21w22=b1N10ub1N10vb2N20vb2N20u,where(50)u=eσa+maa+2m+eα2km2+d+kam+adada+2md+2km,v=eσmaa+2m+eαm2km+d+akada+2md+2km.
By equation (48), ℛ0 is the larger root of the following quadratic equation:(51)λ2−b1N10+b2N20uλ+b1b2N10N20u2−v2=0.
Note that(52)u+v=eσd+αad,u−v=eα+σd+2kma+2md+2km,(53)u2−v2=eσd+αad⋅eα+σd+2kma+2md+2km.
It follows fromu+v=eσd+α/ad that ∂u/∂m=−∂v/∂m. Taking partial derivatives with respect to m for equation (51) gives(54)2λ−b1N10+b2N20u∂λ∂m=λb1N10+b2N20∂u∂m−b1b2N10N20∂u2−v2∂m=λb1N10+b2N20−2b1b2N10N20eα+σdad∂u∂m.
Moreover, from the definition ofu, we have(55)∂u∂m=−eσaaa+2m2−eαad⋅6kmad+8k2m3+ad2+4km2d+ka2d+4k2m2aa+2m2d+2km2<0.
Note that2λλ=ℛ0>b1N10+b2N20u. Thus, to show ∂ℛ0/∂m<0, that is to say, ∂λ/∂mλ=ℛ0<0, it is sufficient to show for λ=ℛ0,(56)b1N10+b2N20λ−2b1b2N10N20⋅eσd+αad>0.
We claim thatλ>1/2b1N10+b2N20u+v. Since(57)λ=12b1N10+b2N20u+b1N10+b2N202u2−4b1b2N10N20u2−v2,it is equivalent to show that(58)b1N10+b2N202u2−4b1b2N10N20u2−v2>b1N10+b2N202v2,or(59)b1N10−b2N202u2−v2>0.
This is automatically true sinceu>v>0. It follows from λ>1/2b1N10+b2N20u+v that(60)b1N10+b2N20λ>12b1N10+b2N202u+v=12eσd+αadb1N10+b2N202.
Thus, as required,(61)b1N10+b2N20λ−2b1b2N10N20eσd+αad>12eσd+αadb1N10−b2N202>0.
This completes the proof.Remark 1.
Under the condition of Theorem4, the following result can be proved similarly; each entry of the matrix W decreases in terms of a or d, and it increases in terms of e, σ, and α. Hence, the net reproduction number decreases in terms of a or d, and it increases in terms of e, σ, and α.


## 4. Persistence and Existence of Positive Equilibria

Under certain conditions, we use the techniques of persistence theory [28, 29] to show the uniform persistence of the predator population and the existence of at least one positive equilibrium (or coexistence equilibrium) when ℛ0>1.Theorem 5.
If equation (13) holds and ℛ0>1, then system (1) is uniformly persistent, that is, there exists a positive constant ϵ such that every solution ΦtA0≡Nt,PJt,PAt with A0∈ℝ+n×Intℝ+2n satisfies(62)liminft⟶∞PiJt,PiAt≥ε,ε,i∈ℕn.
Hence, it has at least one positive equilibrium.Proof of Theorem5.
Let(63)X˜=A0∈ℝ+3nNi0>0,i∈ℕn,X˜0=A0∈X˜PiJ0>0,PiA0>0,i∈ℕn,(64)∂X˜0=X˜∖X˜0.
It then suffices to prove that (1) is uniformly persistent with respect to X˜0,∂X˜0.
It is easy to see from Theorem1 that X˜0 is positively invariant and system (1) is dissipative. Define(65)M∂:=A0∈∂X˜0ΦtA0∈∂X˜0,t≥0.
We claim that(66)M∂=A0∈X˜PiJ0=PiA0=0,i∈ℕn.
Suppose on the contrary that there existsi0∈ℕn such that either Pi0J0>0 or Pi0A0>0, then similar arguments as those in the proof of Theorem 1 gives PiJt>0 and PiAt>0 for t>0 and i∈ℕn. This contradicts with the definition of M∂, and hence the claim has been proved.
Next, chooseε>0 such that sM1+M2−ε>0. Then, choose η>0 such that N0η≥N0−ε/2,…,ε/2︸n and sdiagri−diN−1+σiκη+MN<0. Here, M2ε, N0η, and κ are those defined in the proof of Theorem 2. We claim that(67)limsupt⟶∞maxi∈ℕnPiJt,PiAt>η2.
It is similar to the proof of Theorem2. By way of contradiction, suppose there exists A0 such that limsupt⟶∞maxi∈ℕnPiJt, PiAt≤η/2. Then, there exists t0≥0 such that(68)maxi∈ℕnPiJt,PiAt≤η,fort≥t0.
Then,Nt satisfies (23) for t≥t0. As sdiagri−diN−1+σiκη+MN<0, we have liminft⟶∞Nt≥N0η by Lemma 1 and Comparison Theorem. Thus, there exists t1≥t0 such that Nt≥N0−ε,…,ε︸n for t≥t1. It follows that PJt,PAt satisfies(69)P˙JtP˙At≥M1+M2−εPJtPAt,for t≥t1. Consider the following auxiliary system:(70)P¯˙JtP¯˙At=M1+M2−εP¯JtP¯At.
SincesM1+M2−ε>0, we know λ:=sM1+M2−ε is an eigenvalue of M1+M2−ε with a positive eigenvector, say v¯. As PJt1,PAt1≫0,0, there exists ξ>0 such that PJt1,PAt1≥ξv¯. Again, by Comparison Theorem, we have PJt,PAt≥ξeλt−t1v¯ for t≥t1. It follows that none of PiJt and PiAt, i∈ℕn, is bounded, a contradiction to the boundedness of solutions (see Theorem 1).
Note thatN0 is globally asymptotically stable in ℝ+n∖0 for equation (17). By the above claim, we know that both 0,0,0 and P0 are isolated invariant subsets in X, Ws0,0,0∩A0=∅, where Ws0,0,0 is stable manifold of 0,0,0. Clearly, every orbit in M∂ converges to either 0,0,0 or P0, and 0,0,0 and P0 are acyclic in M∂. By ([30], Theorem 5 and Remark 1), we conclude that system (1) is uniformly persistent with respect to X˜0,∂X˜0. By [31], system (1) has at least one equilibrium N∗,PJ∗,PA∗∈A0 with PJ∗≫0 and PA∗≫0. We further claim that N∗≠0. Otherwise, if N∗=0, by summing up the equilibrium equations for PJ, we get ∑i∈ℕndiJ+αiPiJ∗=0. Thus, PJ∗=0, a contradiction. This completes the proof.


## 5. Simulations

In this section, numerical simulations are used to analyse the effects of the diffusion ratesmijN, mijJ, and mijA on the net reproduction number ℛ0 in a two-patch environment. Here, for having concrete observation, we let mijJ=mijA=m and mijN=mN, i,j=1,2.First, we fix the parameters as follows:(71)r1=0.45,b1=1,σ1=0.01,e1=0.3,α1=0.4,d1J=d1A=0.3,r2=0.5,b2=1,σ2=0.01,e2=0.3,α2=0.4,d2J=d2A=0.3.We change the reproduction number of each isolated patch by changing the carrying capacity and the death rate of prey. Table1 summarizes all the possibilities of the reproduction numbers in both patches.Table 1
The relationship ofℛ01 with ℛ02.
ℛ01<1ℛ01>1ℛ02<1Figure2, ℛ01<1, ℛ02<1, K1<K2Figure4, ℛ01>1, ℛ02<1, K1<K2Figure3, ℛ01<1, ℛ02<1, K1>K2Figure5, ℛ01>1, ℛ02<1, K1>K2ℛ02>1Figure6, ℛ01<1, ℛ02>1, K1<K2Figure8, ℛ01>1, ℛ02>1, K1<K2Figure7, ℛ01<1, ℛ02>1, K1>K2Figure9, ℛ01>1, ℛ02>1, K1>K2WhenmN⟶0, from Figure 2, we find that the net reproduction number ℛ0 increases and then decreases as m increases, and the final value of the threshold is greater than unity, but it is not difficult to observe that ℛ0 is almost unchanged if the predator population hardly diffuses (m⟶0). When the prey population and the predator population simultaneously diffuse, ℛ0 increases and then decreases, and the final value of the threshold is less than unity. In Figure 3, if we reduce the carrying capacity and the death rate of the second isolated patch in Figure 2, that is, K2 decreases from 3 to 2.45 and d2N decreases from 0.23 to 0.15, we draw the conclusion that ℛ0 increases and the final value is greater than unity as m increases no matter how mN changes, and ℛ0 does not change if m⟶0.Figure 2
The effects of the diffusion rates on the net reproduction numberℛ0 where ℛ01=0.96, ℛ02=0.9257, K1=2.8, K2=3, d1N=0.18, d2N=0.23, and mN, mJ=mA=m vary from 0.01 to 1; all other parameters are given in equation (71).Figure 3
The effects of the diffusion rates on the net reproduction numberℛ0 where ℛ01=0.96, ℛ02=0.98, K1=2.8, K2=2.45, d1N=0.18, d2N=0.15, and mN, mJ=mA=m vary from 0.01 to 1; all other parameters are given in equation (71).When we increase the production number of the first isolated patch and haveℛ01>1, ℛ02<1 in Figure 4, it is easy to find that the threshold ℛ0 always decreases and is close to unity as m increases and mN⟶0. On the contrary, ℛ0 first decreases and then increases and is greater than unity as mN increases when m⟶0. However, ℛ0 always decreases and is less than unity as the prey population and the predator population simultaneously diffuse. In Figure 5, we increase the parameters including the carrying capacity and the death rate in the first isolated patch in Figure 4, that is, K1 increases from 2.8 to 3.2 and d1N increases from 0.15 to 0.18. Under these circumstances, we find that the curve of ℛ0 changes a lot. The curve of ℛ0 is waved as m increases when mN⟶0 and is greater than unity when m⟶1. ℛ0 always decreases and tends to stabilize as mN increases when m⟶0. When the dispersal of the prey and the predator happens simultaneously, the curve of ℛ0 is waved and is finally less than unity.Figure 4
The effects of the diffusion rates on the net reproduction numberℛ0 where ℛ01=1.0667, ℛ02=0.9257, K1=2.8, K2=3, d1N=0.15, d2N=0.23, and mN, mJ=mA=m vary from 0.01 to 1; all other parameters are given in equation (71).Figure 5
The effects of the diffusion rates on the net reproduction numberℛ0 where ℛ01=1.0971, ℛ02=0.9257, K1=3.2, K2=3, d1N=0.18, d2N=0.23, and mN, mJ=mA=m vary from 0.01 to 1; all other parameters are given in equation (71).The situations of Figures6 and 7 are similar to Figures 4 and 5 owing to the differences of the carrying capacities and the death rates of the preys. Contrary to circumstances of Figures 4 and 5, we know that the final values of ℛ0 in Figures 6 and 7 are greater than unity as mN,mJ,mA⟶1.Figure 6
The effects of the diffusion rates on the net reproduction numberℛ0 where ℛ01=0.96, ℛ02=1.1657, K1=2.8, K2=3, d1N=0.18, d2N=0.16, and mN, mJ=mA=m vary from 0.01 to 1; all other parameters are given in equation (71).Figure 7
The effects of the diffusion rates on the net reproduction numberℛ0 where ℛ01=0.96, ℛ02=1.008, K1=2.8, K2=2.45, d1N=0.18, d2N=0.14, and mN, mJ=mA=m vary from 0.01 to 1; all other parameters are given in equation (71).The situations of Figures8 and 9 are very interesting; ℛ0 increases and is greater than unity if mN⟶0 or mJ,mA⟶0, but ℛ0 gradually decrease and is even less than unity when mN,m⟶1.Figure 8
The effects of the diffusion rates on the net reproduction numberℛ0 where ℛ01=1.0286, ℛ02=1.12, K1=3, K2=3.5, d1N=0.18, d2N=0.22, and mN, mJ=mA=m vary from 0.01 to 1; all other parameters are given in equation (71).Figure 9
The effects of the diffusion rates on the net reproduction numberℛ0 where ℛ01=1.0286, ℛ02=1.12, K1=3, K2=2.8, d1N=0.18, d2N=0.15, and mN, mJ=mA=m vary from 0.01 to 1; all other parameters are given in equation (71).


## 6. Discussion

In this paper, we studied a predator-prey model to describe the dynamics of predation with stage structure in the predator and dispersal in ann-patch environment. The population dispersal among patches can be interpreted as migration or predation in which the predator has to seek prey by virtue of food reduction. In addition, predation function fN includes three typical types having biological significance.In Section2, we first proposed the predator-prey model with diffusion in an n-patch environment. And we proved that the solutions of system (1) are all positive and ultimately bounded with nonnegative initial conditions (Theorem 1). In Section 3, we derived the net reproduction number and obtained the stability of the PFE through the relationship between ρFV−1 and 1 (Theorem 2). Moreover, we provided the upper bound and the lower bound of ℛ0 (Theorem 3). For a two-patch model (39), the constant dispersal rate of prey and predator decreases ℛ0 (Theorem 4). If ρFV−1>1 and the dispersal rate matrices of prey and predator are assumed to be irreducible, system (1) persists globally (i.e., in all patches) under some conditions in Section 4 (Theorem 5).By numerical simulations, we found some interesting phenomena. From Figures2 and 3, we first know that the net reproduction number ℛ0 may be greater than unity even if ℛ01<1, ℛ02<1, and it may also be less than unity. This phenomenon is related with the carrying capacities K1,K2 and the death rates of prey d1N,d2N. In Figures 4–7, if one of the isolated patch thresholds is greater than unity and the other is less than unity, when the diffusion happens, the curve of ℛ0 changes in a complex way. Here, we still obtain all kinds of figures through changing parameters K1,K2,d1N,d2N. It is not hard to think that the curve of ℛ0 becomes complicated when we change the other parameters. In addition, because dispersal leads to extinction of predators, we believe that dispersal may be a bad thing, but sometimes dispersal may be a good thing if we want to eliminate some pests.There are some interesting questions worthy of further study; for instance, the corresponding conclusions for the above two functional response types still remain open due to the analytical challenge, which remain as future questions.



---

*Source: 6028019-2020-08-08.xml*