A two-patch prey-predator model with predator dispersal driven by the predation strength.

Foraging movements of predator play an important role in population dynamics of prey-predator systems, which have been considered as mechanisms that contribute to spatial self-organization of prey and predator. In nature, there are many examples of prey-predator interactions where prey is immobile while predator disperses between patches non-randomly through different factors such as stimuli following the encounter of a prey. In this work, we formulate a Rosenzweig-MacArthur prey-predator two patch model with mobility only in predator and the assumption that predators move towards patches with more concentrated prey-predator interactions. We provide completed local and global analysis of our model. Our analytical results combined with bifurcation diagrams suggest that: (1) dispersal may stabilize or destabilize the coupled system; (2) dispersal may generate multiple interior equilibria that lead to rich bistable dynamics or may destroy interior equilibria that lead to the extinction of predator in one patch or both patches; (3) Under certain conditions, the large dispersal can promote the permanence of the system. In addition, we compare the dynamics of our model to the classic two patch model to obtain a better understanding how different dispersal strategies may have different impacts on the dynamics and spatial patterns.


(Communicated by Haiyan Wang)
Abstract. Foraging movements of predator play an important role in population dynamics of prey-predator systems, which have been considered as mechanisms that contribute to spatial self-organization of prey and predator. In nature, there are many examples of prey-predator interactions where prey is immobile while predator disperses between patches non-randomly through different factors such as stimuli following the encounter of a prey. In this work, we formulate a Rosenzweig-MacArthur prey-predator two patch model with mobility only in predator and the assumption that predators move towards patches with more concentrated prey-predator interactions. We provide completed local and global analysis of our model. Our analytical results combined with bifurcation diagrams suggest that: (1) dispersal may stabilize or destabilize the coupled system; (2) dispersal may generate multiple interior equilibria that lead to rich bistable dynamics or may destroy interior equilibria that lead to the extinction of predator in one patch or both patches; (3) Under certain conditions, the large dispersal can promote the permanence of the system. In addition, we compare the dynamics of our model to the classic two patch model to obtain a better understanding how different dispersal strategies may have different impacts on the dynamics and spatial patterns. and its related dispersal patterns [57] is considered to be one of the key factors that promote the development of self-organized spatial patterns [1,42,71,72].
In nature, especially for ecological communities of insects, dispersal of a predator is usually driven by its non-random foraging behavior which can often response to prey-contact stimuli [23], including spatial variation in prey density [40] and different types of signals arising directly from prey [76]. For instances, bloodsucking insects respond to the carbon dioxide output and the visual signals of a moving animal, which in tsetse flies (Glossina spp.) lead to the formation of a "following swarm" associated with herds of grazing ungulates [14,34]. Most mosquitoes were attracted over a larger distance by the odor of the host [17,18,19]. The wood-wasp, Sirex noctilio, is attracted by the concentration of the scent [19,53]. Social ants excite "pheromone trails" to encourage other individuals to visit the same food source [6]. Plant-feeding insects commonly detect food items by gustatory signals [65,66,66]. These non-random foraging behaviors driven by prey-mediated patch attractants, prey attractants themselves, and arrestant stimuli following the encounter of a prey, can lead to predation rates that are greater in regions where prey are more abundant (i.e., density-dependent predation), thus regulate population dynamics of both prey and predator.
Recent experimental work on population dynamics of immobile Aphids and Coccinellids by [44] show that the foraging movements of predator Coccinellids are combinations of passive diffusion, conspecific attraction, and retention on plants with high aphid numbers which is highly dependent on the strength of prey-predator interaction. Their study also demonstrates that predation by coccinellids was responsible for self-organization of aphid colonies. Many ecological systems exhibit similar foraging movements of predator. For example, Japanese beetles are attracted to feeding induced plant volatiles and congregate where feeding is taking place [52]. Motivated by these field studies, we propose a two-patch prey-predator model incorporating foraging movements of predator driven by the strength of preypredator interaction, to explore how this non-random dispersal behavior of predator affect population dynamics of prey and predator.
Dispersal of predator plays an important role in regulating, stabilizing, or destabilizing population dynamics of both prey and predator. There are fair amount literature on mathematical models of prey-predator interactions in patchy environments. For example, see work of [47,25,3,24,26,27,64,74,20,21,62,38,39] and also see [41] for literature review. Many studies examine how the interactions between patches affect the synchronicity of the oscillations in each patch, e.g. see the work of [28,35], and how interactions may stabilize or destabilize the dynamics. For instances, [37,35] studied a model with two patches, each with the wellknown prey-predator Rosenzweig-McArthur dynamics, linked by the density independent dispersal (i.e., dispersal is driven by the difference of species' population densities in two patches). His study showed that this type of spatial predator-prey interactions might exhibit self-organization capable of producing stabilizing heterogeneities in prey distribution, and spatial populations can be regulated through the interplay of local dynamics and migration.
However, due to the intricacies that arise in density-dependent dispersal models, there are relatively limited work on models with non-random foraging behavior of predator or non-linear dispersal behavior [40] but see the two patch model with predator attraction to prey, e.g. [32], or predator attraction to conspecific, e.g. [58], or only predators migrate who are attracted to regions with concentrated food resources, see the work of [22,8]. [40] proposed a non random foraging PDE model through a mechanistic approach to demonstrate that area-restricted search does yield predator aggregation, and explore the the consequences of area-restricted search for predator-prey dynamics. In addition, they provided many supporting ecological examples (e.g. Coccinellids, blackbirds, etc.) that abide by their theory. [32] studied a two-patch predator-prey Rosenzweig-MacArthur model with nonlinear density-dependent migration in the predator. The migration term of the predator is derived by extending the Holling time budget argument to migration. Their study showed that the extension of the Holling time budget argument to movement has essential effects on the dynamics. By extending the model of [32], [16] formulated a similar two patch prey-predator model with density-independent migration in prey and density-dependent migration in the predator. Their study shows that several foraging parameters such as handling time, dispersal rate can have important consequences in stability of prey-predator system. [10] investigated the populationdispersal dynamics for predator-prey interactions in a two patch environment with assumptions that both predators and their prey are mobile and their dispersal between patches is directed to the higher fitness patch. They proved that such dispersal, irrespectively of its speed, cannot destabilize a locally stable predatorprey population equilibrium that corresponds to no movement at all.
In this paper, we formulates a new version of Rosenzweig-MacArthur two patch prey predator model with mobility only in predator. Our model is distinct from others as we assume that the non-random foraging movements of predator is driven by the strength of prey-predator interactions, i.e., predators move towards patches with more concentrated prey and predator. Our model can apply to many insects systems such as Aphids and Coccinellids, Japanese beetles and their host plants, etc. For instance, the experimental work of [5] demonstrated that attraction to uninfested potato plants by Colorado potato beetle does not occur when the plants are small. However, when small plants are infested by conspecific larvae they become highly attractive to adult beetles. Thus predators beetles are are more attracted toward patches with high prey-predator interaction strength. The prey-predation attraction can also be observed in the field work of [44]. The main focus of our study of such prey-predator interactions in heterogeneous environments is to explore the following ecological questions: 1. How does our proposed nonlinear density-dependent dispersal of predator stabilize or destabilize the system? 2. How does dispersal of predator affect the extinction and persistence of prey and predator in both patches? 3. How may dispersal promote the coexistence of prey and predator when predator goes extinct in the single patch? 4. What are potential spatial patterns of prey and predator? 5. How are the effects of our proposed nonlinear density-dependent dispersal of predator on population dynamics different from the effects of the traditional density-independent dispersal?
The rest of the paper is organized as follows: In Section (2), we provide the detailed derivation of our two patch prey-predator model. In Section (3), we perform completed local and global dynamics of our model, and derive sufficient conditions that lead to the persistence and extinction of predator as well as permanence of the model. In Section (4), we perform bifurcation simulations to explore the dynamical patterns and compare the dynamics of our model to the traditional model [26,27,37,35,36]. In Section (5), we conclude our study and discuss the potential future study. The detailed proofs of our analytical results are provided in the last section.
2. Model derivations. Let u i (t), v i (t) be the population of prey and predator in Patch i at time t, respectively. In the absence of dispersal, we assume that the population dynamics of prey and predator follow the well-known Rosenzweig-MacArthur prey-predator model as follows: where r i is the intrinsic growth of prey at Patch i; k i is the prey carrying capacity at Patch i; b i is the predator attacking rate at Patch i; h i is the predator handling time at Patch i; c i is the energy conversion rate at Patch i; and δ i is the mortality of predator at Patch i. We assume that the dispersal of predator from one patch to the other is driven by the strength of the prey-predator interaction in two patches which is termed as the attraction strength. More specifically, in the presence of dispersal, the dispersal rate of predators from Patch i to Patch j depends on the prey-predator interaction term bj uj vj 1+bj hj uj in Patch j, and the dispersal term of predators from Patch i to Patch j is modeled by bj uj vj 1+bj hj uj ρ ij v i which gives the net dispersal of predator at Patch i as follows: where ρ ij is the dispersal constant of predators at Patch j moving to Patch i. This assumption of deriving (2) is motivated by the fact that dispersal of a predator is usually driven by its non-random foraging behavior which can often response to prey-contact stimuli [23] which has been supported in many field studies including the recent work by [44]. Note that the dispersal constant ρ ij is directional and not sysmmetrical, i.e., ρ ij may not be equal to ρ ji . In addition, we assume that prey is immobile. This assumption fits in many prey-predator (or plant-insects) interactions in ecosystems such as Aphid and Ladybugs, Japanese beetles and its feeding plants, etc. Based on these assumptions, a two-patch prey-predator model extended from the single patch model (1) is described as the following set of nonlinear equations: We use the similar rescaling approach in [51] by letting Denote that ρ 1 = ρ a2 b2 and ρ 2 = ρ a1 b1 , then the two-patch model (3) could be rewritten as the following equations: where K i is the relative carrying capacity of prey in the absence of predation; a i is the relative predation rate at Patch i; d i is the relative death rate of predator at Patch i; ρ i is the relative dispersal rate of predator at Patch i; and r i is the relative maximum growth rate of prey at Patch i. All parameters are nonnegative.
In addition, we could scale r 1 away by rescaling the time, i.e., t → t r1 . Thus, for 848 YUN KANG, SOURAV KUMAR SASMAL AND KOMI MESSAN the future analysis or simulations, we could let r 1 = 1 and r 2 = r. In summary, the ecological assumptions of Model (4) can be stated as follows: 1. In the absence of dispersal, Model (4) is reduced to the following uncoupled Rosenzweig-MacArthur prey-predator single patch models where r 1 = 1 and r 2 = r and its ecological assumptions [63] can be stated as follows: (a) In the absence of predation, population of prey x i follows the logistic growth model. (b) Predator is specialist (i.e., predator y i goes extinct without prey x i ) and the functional response between prey and predator follows Hollying Type II functional response. 2. There is no dispersal in prey species. This assumption fits in many preypredator (or plant-insects) interactions in ecosystems such as Aphid and Ladybugs, Japanese beetles and its feeding plants, etc. 3. The dispersal of predator from Patch i to Patch j is driven by the preypredation interaction strength in Patch j termed as the attraction strength.
3. Mathematical analysis. The state space of Model (4) is . We have the following theorem regarding the dynamics properties of Model (4): Theorem 3.1. Assume all parameters are nonnegative and r, a i , K i , d i , i = 1, 2 are strictly positive. Model (4) is positively invariant and bounded in R 4 + with lim sup t→∞ x i (t) ≤ K i for both i = 1, 2. In addition, it has the following properties: 1. If there is no dispersal in predator, i.e., ρ 1 = ρ 2 = 0, then Model (4) is reduced to Model (5) whose dynamics can be classified in the following three cases: (a) Model (5) always has the extinction equilibrium (0, 0) which is a saddle.
2 , then the boundary equilibrium (K i , 0) is a saddle; the interior equilibrium (µ i , ν i ) is a source, and the system has a unique limit cycle which is globally asymptotically stable. In addition, the Hopf bifurcation occurs at µ i = Ki−1 2 . 2. The sets {(x 1 , y 1 , x 2 , y 2 ) ∈ R 4 + : x i = 0} and {(x 1 , y 1 , x 2 , y 2 ) ∈ R 4 + : y i = 0} are invariant for both i = 1, 2. If x j = 0, Model (4) is reduced to the single patch model Model (5). If y j = 0, Model (4) is reduced to the following two uncoupled models: where lim t→∞ x j (t) = K j and the dynamics of x i , y i is the same as Model (5).
Notes and biological implications. Theorem 3.1 provides a foundation on our further study of local stability and global dynamics of Model (4). More specifically, Item 2 of Theorem 3.1 implies that Model (4) has the same the invariant sets x i = 0 and y i = 0 for both i = 1, 2 as the single patch models (5). In addition, the results of the single patch models (5) indicate that prey is always persist while predator i is persist if 0 < µ i < K i hold. We would like to point out that the detailed proofs of dynamics of the single patch models (5) has been provided in the work of [30,31,51,77].
3.1. Boundary equilibria and the stability. Now we start with the boundary equilibria of Model (4). Recall that We define the following notations for all possible boundary equilibria of Model (4): .
The following theorem provides sufficient conditions on the existence and stability of these boundary equilibria: [Boundary equilibria of Model (4)] Model (4) always has the following four boundary equilibria where the first three ones are saddles while E K10K20 is locally asymptotically stable if µ i > K i and it is a saddle if (µ 1 − K 1 ) (µ 2 − K 2 ) < 0 or µ i < K i , i = 1, 2. Let i, j = 1, 2, i = j, and is always a saddle. The boundary equilibrium E b i2 is locally asymptotically stable if Ki−1 2 < µ i < K i and one of the following conditions holds: or one of the following conditions holds: In addition, if 0 < µ i < K i for both i = 1 and i = 2, the boundary equilibria E b 12 and E b 22 exist but they cannot be locally asymptotically stable at the same time while if r i = 1, a i = a, d i = d, K i = d for both i = 1, 2, the boundary equilibria E b 12 and E b 22 can not be locally asymptotically stable at all if they exist.
Notes and biological implications. Theorem 3.2 implies the following points regarding the effects of dispersal in predators: 1. Dispersal has no effects on the local stability of the boundary equilibrium E K10K20 . 2. Large dispersal of predator in its own patch may have stabilizing effects from the results of Item sd: In the absence of dispersal, the dynamics of Patch j is unstable at (K j , 0) since 0 < µ j < K j . However, in the presence of dispersal, large values of ρ j can lead to the local stability of the boundary equilibrium E b i2 where i, j = 1, 2 and i = j, under conditions of µ j < K j < di aj −di . 3. Large dispersal of predator in its own patch may have destabilizing effects from the results of Item ub: In the absence of dispersal, the dynamics of Patch j is local stable at (K j , 0) since K j < µ j . However, in the presence of dispersal, large values of ρ j can drive the boundary equilibrium E b j2 being unstable, under conditions of 0 < di aj −d1 < K j < µ j . 4. Under conditions of µ i < K i , the boundary equilibria E b 12 and E b 22 can not be asymptotically stable at the same time.
3.2. Global dynamics. In this subsection, we focus on the extinction and persistence dynamics of prey and predator of Model (4). First we show the following theorem regarding the boundary equilibrium E K10K20 : Notes and biological implications. Theorem 3.3 implies that the dispersal of predators does not effect the global stability of the boundary equilibrium E K10K20 .
To proceed the statement and proof of our results on persistence, we provide the definition of persistence and permanence as follows: Definition 3.4 (Persistence of single species). We say species z is persistent in R 4 + for Model (4) if there exists constants 0 < b < B, such that for any initial condition with z(0) > 0, the following inequality holds where z can be x i , y i , i = 1, 2 for Model (4). Definition 3.5 (Permanence of a system). We say Model (4) is permanent in R 4 + if there exists constants 0 < b < B, such that for any initial condition taken in R 4 + with x 1 (0)y 1 (0)x 2 (0)y 2 (0) > 0, the following inequality holds The permanence of Model (4) indicates that all species in the system are persistence.
Theorem 3.6. [Persistence of prey and predator] Prey x i , i = 1, 2 of Model (4) are always persistent for all r > 0. Predator y j is persistent if one of the following inequalities hold Notes and biological implications. Theorem 3.6 indicates that the dispersal of predators does not affect the persistence of preys, while small dispersal of predator j, under condition of Ki−1 can keep the persistence of predator j. This is consistent with the results of Item uc in Theorem 3.2.
Notes and biological implications. According to Theorem 3.6, we can conclude that Model (4) is permanent whenever both predators are persistent. Theorem 3.7 provides such sufficient conditions that can guarantee the coexistence of bother predator for the two patch model (4), thus provides sufficient conditions of its permanence. Item 1 of this theorem implies that if predator y j is persistent, and large dispersal of predator y i can promote its persistence, thus, promote the permanence since in the absence of dispersals in predator, predator y i goes extinct due to µ i > K i . This is consistent with the results of Item ub of Theorem 3.2 that large dispersal of predator y i can have destabilize effects on the boundary equilibria E b j2 .

YUN KANG, SOURAV KUMAR SASMAL AND KOMI MESSAN
which gives: Since lim sup t→∞ x i (t) ≤ K i for both i = 1, 2 and y i = q i (x i ), therefore, positive solutions of x i ∈ (0, K i ) for (9) determine interior equilibrium of Model (4). By substituting the explicit forms of p i , q i into (9), we obtain the following null clines: with r 1 = 1, r 2 = r and the following properties: Theorem 3.8. [Interior equilibrium] If µ i > K i for both i = 1, 2, then Model (4) has no interior equilibrium. Moreover, we have the following two cases: 1. Assume that a i > a j where i = 1, j = 2 or i = 2, j = 1. Define rj (Kj ai−Kj aj +ai) 2 hold. And it has at least one interior equilibrium (x * 1 , y * 1 , x * 2 , y * 2 ) if the following conditions hold for both i = 1, j = 2 and i = 2, j = 1 In addition, we have aidj aj riρj +aiaj −aidj < x * j < K j for both i = 1, j = 2 and i = 2, j = 1.

Assume that
Notes and biological implications. Theorem 3.8 provides sufficient conditions on the existence of no interior equilibrium when µ i > K i for either i = 1 or i = 2; and at least one interior equilibrium of Model (4) when µ i < K i for both i = 1, 2.
The results indicate follows: (4) has no interior equilibrium if the dispersal of its predator is too small. 2. If µ i < K i for both i = 1, 2, then large values of the predation rate a i , a j and small values of dispersal of both predators can lead to at least one interior equilibrium.
The question is how we can solve the explicit form of an interior equilibrium of Model (4). The following theorem provides us an example of such interior equilibrium of Model (4).
If 0 < µ i < K i for both i = 1 and i = 2, then E i = (µ 1 , ν 1 , µ 2 , ν 2 ) is an interior equilibrium of Model (4) and its stability can be classified in the following cases: 1. E i is locally asymptotically stable if Ki−1 2 < µ i < K i hold for both i = 1 and i = 2 while it is unstable if the following inequality holds 2 ) for either i = 1 or i = 2, then the large values of ρ i can make E i being locally asymptotically stable, i.e.,

YUN KANG, SOURAV KUMAR SASMAL AND KOMI MESSAN
If, in addition, a 1 = a 2 = a, , and E i = (µ, ν, µ, ν) is the only interior equilibrium for Model (4) which has the same local stability as the interior equilibrium (µ, ν) for the single patch model (5) Notes and biological implications. Theorem 3.9 implies Model (4) has an interior equilibrium In addition, Theorem 3.9 indicates that dispersal of predators has no effects on the local stability if Ki−1 2 < µ i < K i for both i = 1, 2 or one of the single patch models (5) However, large dispersal of predator at Patch i can stabilize the interior equilibrium when its single patch model model is unstable at 4. Effects of dispersal on dynamics. From mathematical analysis in the previous sections, we can have the following summary regarding the effects of dispersal of predators for Model (4): 1. Large dispersal of predator at Patch i can stabilize or destabilize the boundary equilibrium of Small dispersal of predator at Patch i may preserve its persistence under certain conditions. On the other hand, large dispersal of predator at Patch i may promote its persistence when the other predator is already persist even if µ i > K i . 3. Dispersal has no effects on the persistence of prey and the number of boundary equilibrium. It has also no effects on the local stability of the boundary equilibrium E K10K20 and the symmetric interior equilibrium (µ, ν, µ, ν) when it exists. 4. If d i > a i , then small dispersal of predator at Patch i prevents the interior equilibrium while if 0 < µ i < K i , large predations rates a i , a j and small dispersal of predators at both patches can lead to at least one interior equilibrium.
To continue our study, we will perform bifurcations diagrams and simulations to explore the effects on the dynamical patterns and compare dynamics of our model (4) to the classical two patch model (12).

Bifurcation diagrams and simulations.
In this subsection, we perform bifurcation diagrams and simulations to obtain additional insights on the effects of dispersal on the dynamics of our proposed two patch model (4). We fix r 1 = 1, r 2 = 1.5, K 1 = 5, K 2 = 3, d 1 = 0.2, d 2 = 0.1. Then according to Theorem 3.1, we know that in the absence of dispersal, the dynamics of Patch 1 has global stability at (5, 0) if 0 < a 1 < 0.24; it has global stability at its unique interior equilibrium 24 < a 1 < 0.3; and it has a unique limit cycle if a 1 > 0.3; while the dynamics of Patch 2 has global stability at (3, 0) if 0 < a 2 < 0.133; it has global stability its unique interior equilibrium 133 < a 2 < 0.2 while it has a unique limit cycle if a 2 > 0.2. Now we consider the following cases: 1. Choose a 1 = 0.25 and a 2 = 0.15. In the absence of dispersal, the dynamics at both Patch 1 and 2 have global stability at its unique interior equilibrium (4, 1), (2, 1.5), respectively. After turning on the dispersal, the coupled two patch model can have one interior equilibrium (see the blue regions in Figure  1(a)) which can be locally stable (see the blue dots in Figure 1(b)), or be a saddle (see the green dots in Figure 1(b)) where the coupled system has fluctuated dynamics; or it can have two interior equilibria (see the red regions in Figure 1 Figure 1(a)). Bifurcation diagrams Figure 1(a),1(b), and 4(b) suggest that dispersal may destabilize system and generate fluctuated dynamics; may generate multiple interior attractors (the case of three interior equilibria), thus generate multiple attractors; or even may drive extinction of predator in one or both patches (he case of two interior equilibria, no interior equilibrium, respectively). 2. Choose a 1 = 0.25 and a 2 = 0.25. In the absence of dispersal, the dynamics of Patch 1 has global stability at its unique interior equilibrium (4, 1) while the dynamics of Patch 2 has a unique stable limit cycle around (2, 1.5). After turning on the dispersal, the coupled two patch model can have one interior equilibrium (see the blue regions in Figure 2(a)) which can be locally stable (see the blue dots in Figure 2(b)), or be a saddle (see the green dots in Figure  2(b)), or be a source (see the red dots in Figure 2(b)) where the coupled system has fluctuated dynamics for the later two cases; or it can have two interior equilibria (see the red regions in Figure 1(a)) which could be two saddles or one sink, one saddle and generate bistability between the interior attractor and the boundary attractor (see Figure 2(b)); or it can have three interior equilibria (see the black regions in Figure 1(a)) which generate multiple interior attractors (see the examples of two sinks and one saddle of Figure 4(b)) or it could have no interior equilibrium (see white regions of Figure 2(a)). Bifurcation diagrams of Figure 2(a), 2(b), and 4(b) suggest that dispersal may stabilize system and generate equilibrium dynamics; may generate multiple interior equilibria (the case of three interior equilibria), thus generate multiple attractors; or even may drive extinction of predator in one or two both patches (the case of two interior equilibria, no interior equilibrium, respectively). 3. Choose a 1 = 0.35 and a 2 = 0.25. In the absence of dispersal, the dynamics of both Patch 1 and 2 have a unique stable limit cycle. After turning on the dispersal, the coupled two patch model can have one interior equilibrium (see the blue regions in Figure 3(a)) which can be a sink (see the red dots in Figure 3(b)) where the coupled system has fluctuated dynamics; or it can have two interior equilibria (see the red regions in Figure 2 suggest that dispersal may generate bistability between the interior attractor and the boundary attractor; or even may drive the extinction of predator in one or both patches (the case of two interior equilibria, no interior equilibrium, respectively).
In summary, Figure 1, 2, 3, and 4 suggest that dispersal of predator may stabilize or destabilize interior dynamics; it may drive the extinction of predator in one or both patches; and it may generate the following patterns of multiple attractors via two or three interior equilibria: 1. Multiple interior attractors through three interior equilibria: In the presence of dispersal, Model (4) can have the following types of interior equilibria and the corresponding dynamics: • Two interior sinks and one interior saddle: Depending on the initial conditions with x 1 (0)y 1 (0)x 2 (0)y 2 (0) > 0, Model (4) converges to one of two sinks for almost all initial conditions (see examples in Figure 1(b)-4(b)). • One interior sink and two interior saddles: Depending on the initial conditions with x 1 (0)y 1 (0)x 2 (0)y 2 (0) > 0, Model (4) either converges to the sink or has fluctuated dynamics for almost all initial conditions (see examples in Figure 1(b), 4(a)). We should also expect the case of one sink v.s. one saddle v.s. one source and the case of two source v.s. one saddle when the interior sink(s) become unstable and go through Hopf-bifurcation. In addition, Model (4) seems to be permanent whenever it processes three interior equilibria. 2. Boundary attractors and interior attractors through three interior equilibria: • one interior sink and one interior saddle: Depending on the initial conditions with x 1 (0)y 1 (0)x 2 (0)y 2 (0) > 0, Model (4)  on the initial conditions, predator at one patch can go extinct when the system has two interior equilibria. In general, simulations suggest that Model (4) is permanent when it processes one or three interior equilibria while it has bistability between interior attractors and the boundary attractors whenever it processes two interior equilibria.

4.2.
Comparisons to the classic model. The dispersal of predator in our model is driven by the strength of prey-predator interactions. This is different from the classical dispersal model such as Model (12) which has been introduced in [36]: where i = 1, j = 2 or i = 2, j = 1 with r 1 = 1, r 2 = r. The symmetric case of Model (12) (i.e., r i = r j , a i = a j , K i = K j , d i = d j , and ρ i = ρ j ) has been  discussed and studied by [36] through simulations of different scenarios of local bifurcation analysis. Jansen's study shows that the classical two-patch model (12) has a rich dynamical behavior where spatial predator-prey populations can be regulated through the interplay of local dynamics and migration: (i) for very small migration rates the oscillations always synchronize; (ii) For intermediate migration rates the synchronous oscillations are unstable and there are periodic, quasi-periodic, and intermittently chaotic attractors with asynchronous dynamics; and (ii) For large predator migration rates, attractors in the form of equilibria or limit cycles exist in which one of the patches contains no prey.
Recently, [51] studied Model (12) with both dispersal in prey and predator. Liu provide global stability of the interior equilibrium for the symmetric case and performed simulations for the asymmetric cases. Here we provide rigorous results on the persistence and permanence conditions that can be used for the comparisons to our Model (4) in the following theorem: Then we have the following summary on the dynamics of Model (12) 1. Model (12) is positively invariant and bounded in its state space R 4 + with lim sup t→∞ x i (t) ≤ K i for both i = 1 and i = 2. 2. Boundary equilibria: Model (12) always has the following four boundary equilibria E 0000 , E K1000 , E 00K20 , E K10K20 where the first three ones are saddles while E K10K20 is locally asymptotically stable if and it is a saddle if one of the above inequalities does not hold. If 0 <μ i < K i , then the boundary equilibrium E b i exists which is locally asymptotically stable if Ki−1 2 <μ i < K i , r j < a jν i j . 3. Subsystem i: If x j = 0, then Model (12) reduces to the following subsystem (13) with three species x i , y i , y j : whose global dynamics can be described as follows: 3a Prey x i is persistent for Model (13) with lim sup t→∞ x i (t) ≤ K i . 3b Model (13) has global stability at (K i , 0, 0) ifμ i > K i . 3c Model (13) has global stability at Persistence of prey: Prey x i persists ifμ j < 0, orμ j > K j , or <ν i j . In addition, under these conditions, Model (12) has global stability at (0,ν j i ,μ j ,ν j ). 6. Persistence and extinction of predators: Predator y i and y j have the same persistence and extinction conditions. Predators persist if 0 < µ i < K i for both i = 1 and i = 2 while both predators go extinct if µ i > K i for both i = 1 and i = 2. In addition, Model (12) has global stability at (K 1 , 0, K 2 , 0) for the later case. 7. Permanence of Model (12): Model (12) is permanent if 0 < µ i < K i for both i = 1 and i = 2 and one of 4(a), 4(b), 4(c) hold. 8. The symmetric case: Let r 1 = r 2 = 1, a 1 = a 2 = a, b 1 = b 2 = b, d 1 = d 2 = d and K 1 = K 2 = K, then µ 1 = µ 2 = µ and ν 1 = ν 2 = ν. Therefore, we can conclude that Model (13) has global stability at (µ, ν, µ, ν) if K−1 2 < µ < K. In addition, the local stability of (µ, ν, µ, ν) for Model (12) is the same as the Scenarios Model (4) whose dispersal is driven by the strength of prey-predator interactions Classical Model (12) whose dispersal is driven by the density of predators E K 1 0K 2 0 LAS and GAS if µi > Ki for both i = 1, 2. Dispersal has no effects on its stability.
GAS if µi > Ki for both i = 1, 2; While > 0. Large dispersal may be able to stabilize the equilibrium.
< µi < Ki and one of the conditions sa, sb, sc, sd in Theorem (3.2) holds. Large dispersal has potential to either stabilize or stabilize the equilibrium.

Does not exists
< µi < Ki and r j (K j +1) 2 4a j K j < ν j i . Large dispersal of predator in Patch i will either destroy or destabilize the equilibrium while large dispersal of predator in Patch j may stabilize the equilibrium. Table 1. The comparison of boundary equilibria between Model (4) and Model (12). LAS refers to the local asymptotical stability, and GAS refers to the global stability.
Notes. Theorem 4.1 indicates follows: 1. If µ i > K i and 0 < µ j < K j , then the large dispersal of predator at Patch i stabilizes E K10K20 . 2. Proper dispersal of predators can drive the extinction of prey in one patch. 3. Dispersal has no effects on the persistence of predator. This is different from our proposed model (12). To see how different types of strategies in dispersal of predators affect population dynamics of prey and predator, we start with the comparison of the boundary equilibria of our model (4) and the classic model (12). Both Model (4) and (12) always have four boundary equilibria E 0000 = (0, 0, 0, 0), E K1000 = (K 1 , 0, 0, 0), E 00K20 = (0, 0, K 2 , 0) and E K10K20 = (K 1 , 0, K 2 , 0) among which first three are saddles in both the cases. If 0 < µ i < K i for i = 1, 2, then Model (4) has another two boundary equilibria E b i1 and E b i2 where E b i1 is a saddle. If 0 <μ i < K i , then Model (4) has another boundary equilibrium E b i . We summarize and compare the dynamics of our model (4) with dispersal in predator driven by the strength of predation and the classical model (12) with dispersal in predator driven by the difference of predator densities in Table 1-Table 3. We highlight effects of dynamical outcomes due to different dispersal strategies in predators between Model (4) and (12) as follows: 1. The boundary equilibria: E K10K20 , E b i2 and E b i . The comparisons listed in Table 1 suggest that dispersal of predator has larger effects on the boundary equilibrium of the classic model than ours. 2. Persistence and extinction of prey. According to the comparison of sufficient conditions leading either persistence or extinction of prey in a patch listed in  Extinction of predator Simulations suggestions (see the yellow regions of Figure 1(a) and Figure 3(a)) that the large dispersal of predator in Patch i may lead to the its own extinction.
Predators in both patches have the same extinction conditions. They go extinct if µi > Ki or µi < 0 for i = 1, 2. Table 3. The comparison of predator persistence and extinction between Model (4) and Model (12).
huge impact on the prey for the classical model (12) but not for our model (4). 3. Persistence and extinction of predator. Simulations and the comparison of sufficient conditions leading either persistence or extinction of predators in a patch listed in Table 3, suggest that the strength of dispersal ability of predator has profound impacts on the persistence of predator for our model (4) while it has no effects on the persistence of predator for the classical model (12). 4. Permanence of a system depends on the persistence of each species involved in the system. Our comparisons of sufficient conditions leading to the persistence of prey and predator listed in Table 2-3, indicate that dispersal of predator has important impacts in the persistence of predator in our model (4) while it has significant effects on the persistence of prey of the classical model (12). We can include that (i) the large dispersal of predator in a patch has potential lead to the extinction of prey (the classical model (12)) or predator (our model (12)) in that patch, thus destroy the permanence of the system; (ii) the small dispersal of predator in Patch i with the large dispersal in Patch j can promote the persistence of prey (the classical model (12)) or predator (our model (12)) in Patch i, thus promote the permanence of the system. 5. Interior equilibria: Both our model (4) and the classical model (12) have the maximum number of three interior equilibria. However, for the symmetric case, our model (4) can have the unique interior equilibrium (see Theorem 3.9) while the classical model can potentially process three interior equilibria [36].
6. Multiple attractors: Both our model (4) and the classical model (12) have two types bi-stability: (a) The boundary attractors where one of prey or predator can not sustain and the interior attractors where all four species can co-exist; and (b) Two distinct interior attractors. One big difference we observed is that for the symmetric case when each single patch model has global stability at its unique interior equilibrium, our model (4) can have only one interior attractor while the classical model can potentially have two distinct interior attractors. This is due to the fact that Model (4) has unique interior equilibrium while Model (12) can potentially process three interior equilibria as we mentioned earlier.

5.
Conclusion. The idea of "metapopulation" originated from [48] where R. Levins used the concept to study the dynamics of pests in agricultural field in which insect pests move from site to site through migrations. Since Levin's work, many mathematical models have been applied to study prey-predator interactions between two or multiples patches that are connected through random dispersion, see examples in [37,36,7,43,35,46,47,59,9,26,2,55,29]. The study of these metapopulation models help us get a better understanding of the dynamics of species interacting in a heterogeneous environment, and allow us to obtain a useful insight of random dispersal effects on the persistence and permanence of these species in the ecosystem. Recently, there has been increasing empirical and theoretical work on the non-random foraging movements of predators which often responses to preycontact stimuli such as spatial variation in prey density [11,40], or different type of signals arising directly from prey [76]. See more related examples of mathematical models in [49,45,12,8,13,22,10,44,16,32,28,56]. Kareiva [41] provided a good review on varied mathematical models that deal with dispersal and spatially distributed populations and pointed out the needs of including non-random foraging movements in meta-population models. Motivated by this and the recent experimental work of immobile Aphids and Coccinellids by [44], we formulate a two patch prey-predator model (4) with the following assumptions: (a) In the absence of dispersal the model reduced to the two uncoupled Rosenzweig-MacArthur preypredator single patch models (5); (b) Prey is immobile; and (c) Predator foraging movements are driven by the strength of prey-predator interaction. We provide basic dynamical properties such positivity and boundedness of our model in Theorem 3.1. Based on our analytic results and bifurcation diagrams, we list our main findings regarding the following questions stated in the introduction how our proposed nonlinear density-dependent dispersal of predator stabilizes or destabilizes the system; how it affects the extinction and persistence of prey and predator in both patches; how it may promote the coexistence ; and how it can generate spatial population patterns of prey and predator: 1. Theorem (3.2) provides us the existence and local stability features of the eight boundary equilibria of our model (4). This result indicates that large dispersal of predator in its own patch may have both stabilizing and destabilizing effects on the boundary equilibrium depending on certain conditions. Theorem (3.3) gives sufficient conditions on the extinction of predator in both patches, which suggest that predator can not survive in the coupled system if predator is not able to survive at its single patch. In this case, dispersal of predator has no effect on promoting the persistence of predator but dispersal may drive predator extinct even if predator is able to persist at the single patch state (see white regions of Figure 1(a), 2(a), and 3(a)). 2. Theorem (3.6) provides sufficient conditions of the persistence of prey and predator while Theorem (3.7) provides sufficient conditions of the permanence of our two patch model. These results imply that under certain conditions, large dispersal of predator can promote its persistence, thus, promote the permanence of the coupled system while predator in that patch goes extinct in the absence of dispersal. Our numerical studies also suggests that large dispersal can also drive the extinction of predators in both patches (see white regions of Figure 1(a), 2(a), and 3(a)). 3. Theorem (3.8) and Theorem (3.9) provide sufficient conditions on the existence and the local stability of the interior equilibria under certain conditions. Our analytic study shows that large dispersal of predator may be able to stabilize the interior equilibrium when one of the single patch has global stable interior equilibrium while the other one has limit cycle dynamics. At the mean time, our bifurcation diagrams (see Figure 1(b), 2(b), and 2(b)) suggest that the stabilizing or destabilizing effects of predator's dispersal are not definite, i.e., dispersal can either stabilize or destabilize the system depending on other life history parameters. Moreover, our simulations also suggest that the dispersal of predator can either generate multiple interior equilibria or destroy the interior equilibrium which leads to the extinction of predator in one patch or predators in both patches.
Comparisons to the classic model (12). We provide detailed comparison between the dynamics of our model (4) to the classic model (12). These comparisons suggest that the mode of forging movement of predator has profound impacts on the dynamics of the coupled two patch model. Here we highlight two significant differences: (1) the strength of dispersal ability of predator has profound impacts on the persistence of predator for our model (4) while it has no effects on the persistence of predator for the classical model (12). However, the dispersal of predator has huge impacts on the persistence of prey for the classical model (12) while it has little or no effects on the persistence of prey for our model (4). And (2) for the symmetric case, our model (4) has a unique interior equilibrium while the classical model (12) can have up to three interior equilibria thus it is able to generate different spatial patterns.
Proof. Notice that both dxi dt xi=0 = 0 and dyi dt yi=0 = 0 for i = 1, 2, thus according to Theorem A.4 (p.423) in [73], we can conclude that the model (4) is positive invariant in R 4 + . Now we can go ahead to show the boundedness of the system. First, we have the following inequalities due to the property of positive invariance: Therefore, we have (5) is bounded in R 4 + . If there is no dispersal in predator, i.e., ρ i = 0, i = 1, 2, we can easily check that Model (4) is reduced to the two uncoupled Rosenzweig-MacArthur prey-predator single patch models (5) with r 1 = 1 and r 2 = r. The global dynamics of the single patch model (5) can be summarized from the work of [50,51,30,31]. Thus, we omit the detailed proof here.
Summarizing the discussions above, we can conclude that the statement of Theorem 3.1 holds.
YUN KANG, SOURAV KUMAR SASMAL AND KOMI MESSAN (14) After substituting the boundary equilibria E 0000 , E K1000 , E 00K20 , E µ1ν100 and E 00µ2ν2 into the Jacobian Matrix (14), we can conclude that these equilibria are saddles since they have both positive and negative eigenvalues.
The eigenvalues of (14) evaluated at E K10K20 are as follows: Therefore, E K10K20 is locally asymptotically stable if µ i > K i , i = 1, 2 while it is a saddle if either (µ 1 − K 1 ) (µ 2 − K 2 ) < 0 or µ i < K i , i = 1, 2 holds. Now we focus on the local stability of E µ1ν1K20 and E K10µ2ν2 when they exist. After substituting the boundary equilibrium E µ1ν1K20 to (14), we can obtain the eigenvalues of the Jacobian matrix evaluated at this boundary equilibrium as follows: and Notice that the eigenvalues of λ 3 and λ 4 being negative is equivalent to the case that the unique interior equilibrium (µ 1 , ν 1 ) being locally asymptotically stable for the single patch model (5) when i = 1. Thus, we can conclude that K1−1 2 < µ 1 < K 1 are sufficient conditions for λ 3 and λ 4 being negative. Now we explore sufficient conditions for λ 2 being negative. First, we have µ 1 < K 1 due to the existence of E µ1ν1K20 . We have the following three cases: 1. If µ 2 > K 2 ⇔ K 2 (a 2 − d 2 ) − d 2 < 0, then the first term of λ 2 is negative. This also implies that Model (4) has no boundary equilibria of E 00µ2ν2 and E K10µ2ν2 . Since µ 1 < K 1 , therefore, we have λ 2 < 0 for all Therefore, we can conclude that λ 2 is negative if either a 2 ≤ d 1 , a2−d1 > 0), then λ 2 < 0 if ρ 2 is small enough, i.e., satisfies the condition of In this case, we can conclude that λ 2 is negative if 2. If µ 2 < K 2 ⇔ K 2 (a 2 − d 2 ) − d 2 > 0, then the first term of λ 2 is positive. This also implies that Model (4) has two boundary equilibria of E 00µ2ν2 and E K10µ2ν2 . In this case, sufficient conditions for λ 2 being negative are K 2 < d1 a2−d1 and ρ 2 large enough. More specifically, ρ 2 has to satisfy the following inequality: . Summarizing the discussions above, we can conclude that the boundary equilibrium E µ1ν1K20 is locally asymptotically stable if K1−1 2 < µ 1 < K 1 and one of the following conditions holds: 1.
or one of the following conditions holds: . Similarly, we can obtain sufficient conditions for the local stability of the boundary equilibrium E K10µ2ν2 as the statement.
If µ i < K i , then Model (4) has the boundary equilibria E µ1ν1K20 and E K10µ2ν2 according to Theorem 3.2 and the discussions above. If both E µ1ν1K20 and E K10µ2ν2 are locally stable, then the following inequalities are satisfied: which are contradiction. Therefore, E µ1ν1K20 and E K10µ2ν2 can not be local stable at the same time.
This implies that one of the eigenvalues of the Jacobian matrix of Model (4) evaluated at E b i2 is positive, i.e., can not be stable for both i = 1 and i = 2. Proof of Theorem 3.3.
aiKi , then we have We construct the following Lyapunov functions and Now taking derivatives of the functions (15) and (16) with respect to time t, we get and . (18) Let V = V 1 + V 2 . Now adding (17) and (18), we get positive for x i > K i and it is negative for x i < K i . At the mean time, we have q i (x i ) is positive for x i < K i and it is negative for This implies that both V 1 and V 2 are Lyapunov functions, and the boundary equilibrium E K10K20 = (K 1 , 0, K 2 , 0) is globally stable when µ i > K i according to Theorem 3.2 in [31].
Proof of Theorem 3.6.
Proof. According to Theorem 3.1, we know that Model (4) is attracted to a compact set C in R 4 + . Moreover, if y j = 0, Model (4) is reduced to the two uncoupled models (7) while if x j = 0 it is reduced to a single patch model (5).
First we focus on the persistence conditions for prey x 1 . Model (4) is reduced to a single patch model (5) when x 1 (0) = 0, i.e., we have x 1 = y 1 = 0. Notice that According to Theorem 2.5 of [33], we can conclude that prey x 1 is persistent. Similarly, we can show that prey x 1 is persistent for all r > 0.
Since both x 1 and x 2 are persistent, then we can conclude that Model (4) is attracted to a subcompact set C s of C that excludes E 0000 , E K1000 and E 00K20 . Therefore, we can restrict the dynamics of Model (4) on the compact set C s . Now we focus on the persistence conditions for predator y 1 . According to Theorem 3.1, if y 1 = 0, Model (4) is reduced to the two uncoupled models (7). In this case, according to both Theorem 3.1 and 3.2, the omega limit sets of (4) on the compact set C s are E K1000 , , E K10K20 , E K10µ2ν2 if K2−1 2 < µ 2 < K 2 while they are E K1000 , E K10K20 if µ 2 > K 2 . Now we consider the following two cases: 1. If µ 2 > K 2 , according to Theorem 2.5 of [33], we can conclude that predator y 1 is persistent if all of the following equations are strictly positive: Since a1K1 1+K1 − d 1 > 0 ⇔ µ 1 < K 1 , therefore, we can conclude that predator y 1 is persistent if µ 1 < K 1 and µ 2 > K 2 . 2. If K2−1 2 < µ 2 < K 2 , according to Theorem 2.5 in [33] and discussions above, we can conclude that predator y 1 is persistent if µ 1 < K 1 and the following equation is strictly positive: According to the proof of Theorem 3.2, we can see that sufficient condition that dy1 y1dt E K 1 0µ 2 ν 2 > 0 holds is the same as sufficient condition for the boundary equilibrium E K10µ2ν2 being unstable when µ 1 < K 1 . Therefore, we can conclude that predator y 1 is persistent if one of the following inequalities hold (a) µ j < K j , µ i > K i . Or . where i = 1, j = 2 or i = 2, j = 1.one of the following inequalities hold.
Based on the discussion above, we can conclude that the statement of Theorem 3.6 holds.
Proof of Theorem 3.7.

Proof. If
Kj −1 2 < µ j < K j , µ i > K i , then according to Theorem 3.6, we can conclude that prey x i for both i = 1, 2 and predator y j is persistent. This implies that Model (4) is permanent if predator y i is persistent. Since Kj −1 2 < µ j < K j , then Theorem 3.2 indicates that the omega limit set of Model (4) when y i = 0 is E µ1ν1K20 when i = 2, j = 1 while its omega limit set is E K10µ2ν2 when i = 2, j = 1. Now let i = 1, j = 2, then according to Theorem 2.5 of [33], we can conclude that predator y 1 is persistent if the following equation is strictly positive:
Similarly, we can show that predator y 2 is persistent when i = 2, j = 1. Therefore, Model (4) is permanent if the following inequalities hold for either i = 2, j = 1 or i = 1, j = 2, According to Theorem 3.6, we can conclude that prey x i for both i = 1, 2 and predator y i is persistent if the following inequalities hold Therefore, Model (4) is permanent if the above inequalities hold for both i = 1, j = 2 and i = 2, j = 1. On the other hand, predator y j is persistent if the following inequalities hold .
Therefore, both predator y i and y j are persistent if the following inequalities hold for either i = 1, j = 2 or i = 2, j = 1, Based on the discussion above, we can conclude that the statement of Theorem 3.7 holds.
Proof of Theorem 3.8.
Proof. If µ i > K i for both i = 1, 2, then Model (4) has global stability at (K 1 , 0, K 2 , 0) according to Theorem 3.3. This implies that Model (4) has no interior equilibrium when µ i > K i for both i = 1, 2.
Notice that the nullclines x 1 = F (x 2 ) = ft(x2) f b (x2) and x 2 = G(x 1 ) = gt(x1) g b (x1) has the following properties: The discussion so far also indicates that we have both f b (x 2 ) > 0 for x 2 ∈ [0, K 2 ] and g b (x 1 ) > 0 for x 1 ∈ [0, K 1 ] if the following inequalities hold for i = 1, j = 2 or i = 2, j = 1 Now assume that these conditions hold, then we have F (x 2 ) and G(x 1 ) are positive on their restricted domain. By algebraic calculations, if a i > max{d 1 , d 2 } for both i = 1, 2, then both F (x 2 ) and G(x 1 ) have its unique critical points x c i , i = 1, 2 in their restricted domain where If F (x c 2 ) < K 1 and G(x c 1 ) < K 2 , then we can conclude that both maps x 1 = F (x 2 ) and x 2 = G(x 1 ) are unimode and the skew product of F × G maps [0, K 2 ] × [0, K 1 ] to its compact subset. Since both F and G are continuous and differentiable, therefore, x 1 = F (x 2 ) and x 2 = G(x 1 ) has at least one positive intersection for x 2 ∈ [0, K 2 ], x 1 ∈ [0, K 1 ]. Now we focus on sufficient condition that leads to F (x c 2 ) < K 1 and G(x c 1 ) < K 2 when a 1 > a 2 . Since , therefore, we have F (x c 2 ) < K 1 and G(x c 1 ) < K 2 when a 1 > a 2 if the following inequalities hold and Therefore, we can conclude that Model (4) has at least one interior equilibrium (x * 1 , y * 1 , x * 2 , y * 2 ) if the following inequalities hold In addition, since both F (x 2 ) and G(x 1 ) are unimode maps in their domain with unique local maximum, thus, we have Therefore, we have aidj aj riρj +aiaj −aidj < x * j < K j for both i = 1, j = 2 and i = 2, j = 1.
Applying the similar arguments for the case a i > a j , we have Therefore, we can conclude that Model (4) has at least one interior equilibrium (x * 1 , y * 1 , x * 2 , y * 2 ) if the following inequalities hold a 1 = a 2 = a > max{d 1 , d 2 }, ρ i < 4(K i a − K i d i − d i ) K j r j for both i = 1, j = 2 and i = 2, j = 1. In addition, dj riρj +a−dj < x * j < K j hold for both i = 1, j = 2 and i = 2, j = 1.
Therefore, both predators go extinct if µ i > K i for both i = 1 and i = 2. Since both lim sup t→∞ y i (t) = 0 for both i = 1, 2. Then we have Model (12) reduced to the following uncoupled prey model which converges to x i = K i . Thus, Model (13) has global stability at (K 1 , 0, K 2 , 0) when µ i > K i for both i = 1, 2.