A two patch prey-predator model with multiple foraging strategies in predators: Applications to Insects

We propose and study a two patch Rosenzweig-MacArthur prey-predator model with immobile prey and predator using two dispersal strategies. The first dispersal strategy is driven by the prey-predator interaction strength, and the second dispersal is prompted by the local population density of predators which is referred as the passive dispersal. The dispersal strategies using by predator are measured by the proportion of the predator population using the passive dispersal strategy which is a parameter ranging from 0 to 1. We focus on how the dispersal strategies and the related dispersal strengths affect population dynamics of prey and predator, hence generate different spatial dynamical patterns in heterogeneous environment. We provide local and global dynamics of the proposed model. Based on our analytical and numerical analysis, interesting findings could be summarized as follow: (1) If there is no prey in one patch, then the large value of dispersal strength and the large predator population using the passive dispersal in the other patch could drive predator extinct at least locally. However, the intermediate predator population using the passive dispersal could lead to multiple interior equilibria and potentially stabilize the dynamics; (2) For symmetric patches (i.e., all the life history parameters are the same except the dispersal strengths), the large predator population using the passive dispersal can generate multiple interior attractors; (3) The dispersal strategies can stabilize the system, or destabilize the system through generating multiple interior equilibria that lead to multiple attractors; and (4) The large predator population using the passive dispersal could lead to no interior equilibrium but both prey and predator can coexist through fluctuating dynamics for almost all initial conditions.


Introduction
The dispersal of an individual has consequences not only for individual fitness, but also for population dynamics and genetics, and species' distributions (Bowler and Benton, 2005;Clobert et al., 2001;Gilpin and Hanski, 1991;Hanski, 1999). As the impact of dispersal on population dynamics has been increasingly recognized, understanding the link between dispersal and population dynamics is vital for population management and for predicting how population responses to changes in the environment. For many animals and insects, the costs and benefits of dispersal will vary in space and time, and among individuals. Thus, the profit of the dispersal ability as a life-history strategy will vary as a result, and a plastic dispersal strategy is typically expected to respond to this variation (Bowler and Benton, 2005;Ims and Hjermann, 2001;Massot et al., 2002;Ronce et al., 2001). The varied dispersal driving forces include population density, kin selection relatedness, conspecific attraction, interspecific interactions, food availability, patch size and qualities, etc. There has been a large number of empirical studies supporting the effects of various parameters on dispersal mechanisms and strengths (Bowler and Benton, 2005). For example, the field work by Kiester and Slatkin (1974) showed evidence of Iguanid lizards that encompass two or more dispersal strategies as foraging movements. Kummel et al. (2013) showed through their field work that the foraging behavior of Coccinellids are governed not only by the conspecific attraction but also through the passive diffusion and retention on plants with high immobile aphids number. The main purpose of this article is to investigate the effects of the combinations of different strategies on population dynamics of a prey-predator interaction model when prey is immobile.
Due to the practical difficulties associated with the field study of dispersal, theoretical studies play a particularly important role in predicting the effects of varied dispersal strategies in population dynamics (Bowler and Benton, 2005). The patchy prey-predator population models with different dispersal forms have been proposed and studied in a fair amount of literature. For example, the work of (Fraser and Cerri, 1982;Hansson, 1991;Jánosi and Scheuring, 1997;Namba, 1980;Nguyen-Ngoc et al., 2012;Savino and Stein, 1989;Silva et al., 2001) explored the effects of dispersal on population dynamics of prey-predator models when local population density is a selecting factor for dispersal. The work of Huang and Diekmann (2001) and Ghosh and Bhattacharyya (2011) studied the population dynamics of a two patch model with dispersal in predator driving by local population density of prey through Holling searching-handling time budget argument. The work of Kareiva and Odell (1987) studied dynamics when the dispersal of predator is carried out due to the concentrated food resources. Cressman and Křivan (2013) investigated a two patch population-dispersal dynamics for predator-prey interactions with dispersal directed by the fitness. Recent work of (Kang et al., 2014) studied a two patch prey-predator model where predator is dispersed to the patch with the stronger strength of prey-predator interaction. These theoretical work provide useful insights on the link of dispersal strategies and prey-predator population dynamics.
Many empirical work of animal and insects show that dispersal strategies vary among species according to their life history and how they interact with the environment (Bowler and Benton, 2005). However, there is a limited theoretical work on studying how combinations of different dispersal strategies affect population dynamics of prey-predator models in the patchy environment. This paper presents an extended version of a Rosenzweig-MacArthur two patch prey-predator model studied in (Kang et al., 2014) where prey is immobile and the dispersal of predator is attracted by the strength of prey-predation interaction. Our proposed model is motivated by the field experiments of (Kiester and Slatkin, 1974;Kummel et al., 2013;Stamps, 1988). The current model integrates the two dispersal strategies of predator: (1) the passive dispersal, i.e., the classical foraging behavior where predator is driven to the patch with the lower predator population density (e.g. (Jansen, 1995)); (2) the density dependent dispersal measured through the predation attraction (Kang et al., 2014). The linear combination of these two strategies is linked through a parameter whose value is between 0 and 1, and measures the proportion of the predator population using these two dispersal strategies. We aim to use our model to explore how the combinations of these two dispersal strategies of predator affect population dynamics of prey-predator interaction.
The paper is organized as follows: Section 2 introduces the proposed model along with its biological derivation, and provide a brief summary on the dynamics of the related subsystems. Section 3 presents mathematical analysis of the local and global dynamics of the proposed model. Section 4 Investigates the effects of dispersal strategies through bifurcation diagrams. Section 5 concludes our findings along with the related potential biological interpretations.

Model derivations and the related dynamics
Let x i (t), y 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 Rosenzweig-MacArthur prey-predator model. The dispersal of predator from Patch i to Patch j is driven by two mechanisms. The first mechanism relies on the strength of the prey-predation interaction in Patch j (also called "the predation strength"). Let ρ i represents the relative dispersal rate of predator at Patch i, then we obtained the following net predation attraction driven dispersal of predator at Patch i This assumption follows directly from the experimental work of Stamps (1988) in which he concluded that Anolis aeneus juveniles are attracted to conspecific territorial residents under natural conditions in the field. This assumption has also been supported by many field studies including (Alonso et al., 2004;Auger and Poggiale, 1996;Hassell and Southwood, 1978).
The second dispersal mechanism is termed as "the passive dispersal" in which the dispersal is driven by the local population density of predator. The effects of this dispersal strategy has been well studied by many researchers (Hastings, 1983;Jánosi and Scheuring, 1997;Jansen, 1995;Matthysen, 2005;Namba, 1980;Nguyen-Ngoc et al., 2012;Poggiale, 1998;Silva et al., 2001). For example overcrowding of predator in a patch may decrease the resource assessment that can constitute a cue for for the local predators to move. Following this inference, the net dispersal of predators from Patch i to Patch j is given by Motivated by the field work of Kiester and Slatkin (1974) on Iguanid lizards and Kummel et al. (2013) on Coccinellids, we incorporate these two dispersal strategies above into our model. After similar rescaling approach by Liu and Chen (2003), our proposed model is presented as follows with r 1 = 1, r 2 = r being the relative intrinsic growth rates, K i being the relative carrying capacity of prey at Patch i in the absence of predation, d i being the death rate of predator in Patch i, and the parameter s ∈ [0, 1] representing the proportion of predator population using the passive dispersal strategy: First, we have the following theorem regarding the basic dynamic properties of Model (1): Our main focus is to explore how the combination of two different dispersal strategies measured by the parameter s ∈ [0, 1] affect the two patch population dynamics. Before we continue, we first provide a summary of the dynamics of the subsystems of Model (1) including the cases of s = 0 and s = 1.
, 2, then in the absence of dispersal in predator, Model (1) is reduced to the following Rosenzweig and MacArthur (1963) prey-predator single patch models i = 1, 2 with r 1 = 1 and r 2 = r and its global dynamics which can be summarized from the work of (Hsu et al., 1977;Hsu, 1978;Liu and Chen, 2003) as follows: 1. Model (2) always has two boundary equilibria (0, 0), (K i , 0) where the extinction (0, 0) is always a saddle.
2. The boundary equilibria (K i , 0) is globally asymptotically stable if µ i > K i .
3. If K i −1 2 < µ i < K i , then (K i , 0) becomes saddle and the unique interior equilibria (µ i , ν i ) emerges which is globally asymptotically stable. 4. If 0 < µ i < K i −1 2 , the boundary equilibrium (K i , 0) is a saddle, and the unique interior equilibrium (µ i , ν i ) is a source where Hopf bifurcation occurs at µ i = K i −1 2 . The system (2) has a unique stable limit cycle.
The summary on the dynamics of Model (1) when the dispersal of predator foraging activities is driven by local population density (i.e., s = 1) and when the dispersal of predator foraging activities is driven by predation strength (i.e. s = 0) are briefly presented in Table 3 (see Kang et al. (2014) for more detailed summary on the global dynamics).

Mathematical analysis
From Theorem 2.1, we know that the set {(x 1 , y 1 , x 2 , y 2 ) ∈ R 4 + : x i = 0}, is invariant for both i = 1, 2. Assume that x j = 0, Model (1) is reduced to the following three species subsystem: whose basic dynamics are provided in the following theorem: 2. If µ i > K i , then predators in two patches go extinct, and the system (3) has global stability at (K i , 0, 0).

If
, then predators in the two patches are persistent.
Notes: Model (3) can apply to the case where Patch i is the source patch with prey population and Patch j is the sink patch without prey population. The predator in the sink patch is migrated from the source patch. Theorem 3.1 indicates the follows regarding the effects of the proportion of predator using the passive dispersal on Model (3): 1. Prey x i of Model (3) is always persistent for all r i > 0. This is different than the case of s = 1 since prey may go extinct when s = 1.
2. If µ i < K i and ρ i s is small enough, then the inequality holds, hence predators persist. This result suggests that, under the condition of µ i < K i , the large value of ρ i s could drive predator extinction in two patches at least locally.
The interior equilibria (x * 1 , y * 1 , y * 2 ) of Model (3) is determined by first solving for y * i and x * i in dx i dt = 0 and dy j dt = 0 as follows: An equation of y * j is obtained by solving the following equation from Model (3): A substitution of y * i from (4) into y * j gives y * j = The discussion above implies that the existence of the interior equilibrium requires a i > d i and Then we can conclude that x * i solving from Equation (4) is in term of y * i and y * j . Upon substitution of y * i and y * j into x * i we obtain the following nullclines: Based on the arguments above and additional analysis, we have the following proposition regarding the existence of the interior equilibria of Model (3): If, in addition, µ i < x * iℓ < K i for both ℓ = 1, 2, then Model (3) has two interior equilibria.
Notes: Proposition (3.1) implies that even if f i (x i ) has two positive real roots, Model (3) may have none or one interior equilibrium unless these two positive roots are in (µ i , K i ). Note that the interior equilibria of the subsystem Model (3) represent the boundary equilibria of Model (1) when x 1 = 0(i = 2) or x 2 = 0(i = 1). The existence of these boundary equilibria of Model (1) when x 1 = 0 or x 2 = 0 are hence guarantee by the conditions to obtain the interior equilibria E ℓ x i ,y i ,y j and E ℓ y j ,x i ,y i from Proposition (3.1).
These fixed values implies that at Patch 2, prey and predator coexist in the form of a unique stable limit cycle in the absence of dispersal since µ 2 = d 2 a 2 −d 2 = 35/105 < (K 2 −1)/2 = 3. We consider the following two typical cases regarding the population dynamics of prey and predator in the absence of dispersal: 1. d 1 = 0.85, a 1 = 1: Predator and prey are persistent and have global equilibrium dynamics at Patch 1 in the absence of dispersal since (K 1 − 1)/2 = 4.5 < µ 1 = d 1 a 1 −d 1 = 17/3 < 10 = K 1 . 2. d 1 = 2, a 1 = 2.1: Predator goes extinct globally at Patch 1 in the absence of dispersal since The fixed values of parameters and the two cases above provide the following four scenarios: 1. i = 1 (i.e., x 2 = 0 for Model (1)) with d 1 = 0.85, a 1 = 1. In this case, Patch 1 is the source patch and Model (3) can have up to two interior equilibria depending on the values of s (see Figure 1(a)).
3. i = 2 (i.e., x 1 = 0 for Model (1)) with d 1 = 0.85, a 1 = 1. In this case, Patch 2 is the source patch and Model (3) can have up to two interior equilibria depending on the values of s (see Figure 1(b)). The relative large value of s can stablize the dynamics (see the blue region of Figure 1(b)).
4. i = 2 (i.e., x 1 = 0 for Model (1)) with d 1 = 2, a 1 = 2.1. In this case, Patch 2 is the source patch and Model (3) can have up to two interior equilibria depending on the values of s (see Figure 1(c)). The relative large value of s can stablize the dynamics (see the blue region of Figure 1(c)).
The bifurcation diagrams ( Figure 1) suggest that the proportion of predators using the passive dispersal can have huge impacts on the number of interior equilibria of Model (3): For the small values of s, Model (3) can have one interior equilibrium (E 1 x 1 ,y 1 ,y 2 or E 1 y 1 ,x 2 ,y 2 ); For the intermediate values of s, Model (3) can have two interiors E l x 1 ,y 1 ,y 2 , l = 1, 2 (i = 1) or E l y 1 ,x 2 ,y 2 , l = 1, 2 (i = 2); For the large values of s, it has no interior equilibria. A more detail description of the effects of s on the interior equilibria of Model (3) is provided in Table (1).    (3) with y-axis representing the population size of predator at Patch 1 and x-axis represent the proportion of predator using the passive dispersal. Figure 1(a) describes the number of interior equilibria (x * 1 ,ŷ * 1 ,ŷ * 2 ) when x2 = 0 in Model (1) and their stability with respect to variation in s. Figure 1(b) and 1(c) describe the number of interior equilibria (y * 1 , x * 2 , y * 2 ) of the submodel x2 = 0 of Model (1) and their stability when s varies from 0 to 1. Blue represents the sink and green represents the saddle. (1) First, we have boundary equilibria and global dynamics of Model (1) in the following theorem. (1)] Assume that s ∈ (0, 1). Model (1) always has the following four boundary equilibrium

Theorem 3.2. [Boundary equilibria and global dynamics of Model
with the first three always being saddle. E K 1 0K 2 0 is locally asymptotically stable if the follwoing two inequalities in (7) hold: And E K 1 0K 2 0 is saddle when one or both of equations (7) are not satisfied. In addition, (1) is persistent, and the predator population in each patch is persistent if µ i < K i for both i = 1, 2.

At least prey population in one patch of Model
Notes: Theorem 3.2 indicates that the global stability of the boundary equilibrium E K 1 0K 2 0 does not depend on the proportion of predator population using the passive dispersal since E K 1 0K 2 0 is globally asymptotically stable when µ i > K i , i = 1, 2 which is independent of s. However, the value of s > 0 and ρ i , i = 1, 2 can stabilize E K 1 0K 2 0 . For example, assume that µ i < K i and µ j > K j , then in the absence of dispersal, the boundary equilibrium E K 1 0K 2 0 is a saddle. In the presence of the dispersal, according to Theorem 3.2, if we choose ρ j large enough, then E K 1 0K 2 0 can be locally stable, thus the large dispersal at one patch may stabilize the boundary equilibrium E K 1 0K 2 0 . However, if s = 0, then dispersal has no such effects.
Under these parameter values, we have the following two cases that are shown in   Table 2).
We recapitulate the following dynamics regarding the effect of s on the equilibria E b 1ℓ and E b 2ℓ , ℓ = 1, 2: (1) Model (1) can have up to four boundary equilibria; (2) These boundary equilibria when exist are locally asymptotically stable or saddle; (3) Large s has a potential to destroy these equilibria. Also, observe the blue line for locally stable and green line for saddle in Figure 1(a) as oppose to only green line for saddle in Figure 3(a); this results suggest that the additional dimension from the three species Model (3) has a destabilization effect on the four species Model (1).   (1) and their stability with respect to variation in s when d1 = 0.85, a1 = 1. Figure 3 (1) and their change in stability when s varies from 0 to 1 with d1 = 0.85, a1 = 1 and d1 = 2, a1 = 2.1 respectively. Blue represents the sink and green represents the saddle. a 1 = 1 and d 1 = 0.85 a 1 = 2.1 and  (1) we have the following equations

Interior equilibria and stability of Model (1)
Consider (x * 1 , y * 1 , x * 2 , y * 2 ) as an interior equilibrium of Model (1), then the following conditions must be satisfied: which yields the following by substituting the expression of p i (x) and q i (x) into (8) The equation (9) gives the following nullclines: The complex form of (10) prevents us to obtain the explicit solutions of the interior equilibria of Model (1). We are going to explore the symmetric interior equilibrium for the symmetric Model (1) where we say that Model (1) is symmetric if a 1 = a 2 = a, d 1 = d 2 = d, K 1 = K 2 = K, r 1 = r 2 = r. Now we have the following theorem: Then E = (µ, ν, µ, ν) is an unique symmetric interior equilibrium for Model (1). Moreover, E is locally asymptotically stable if K−1 2 < µ < K while it is unstable if µ < K−1 2 for s ∈ [0, 1].
Notes: Theorem (3.3) implies the symmetric Model (1) has an unique symmetric interior equilibrium of the form E = (µ, ν, µ, ν). The related results imply that dispersal of predators and s has no effect on the local stability of this symmetric interior equilibria when it exist since K−1 2 < µ < K does not depend on ρ i , i = 1, 2 or s. We note that Model (1) can have two additional interior equilibria in the symmetric case which can be locally stable or saddle depending on the value of s (see green line for saddle and blue line for locally stable in Figures 4(a) which correspond to the additional two boundary equilibria of Model (1) in the symmetric case). We consider the following fixed symmetric parameters: r 1 = r 2 = r = 1, d 1 = d 2 = d = 5, K 1 = K 2 = K = 10, a 1 = a 2 = a = 6.
According to the bifurcation diagrams in Figures 4(a) and 4(b), Model (1) can have up to three interior equilibria in the symmetric case. It seems that the larger value of s can create two additional asymmetric interior equilibria which can be saddle or locally stable, thus generate bistability between two different interior attractors (See blue lines in Figure 4(a) when 0.78 ≤ s ≤ 0.92). The local stability of E = (µ, ν, µ, ν) does not depend on s as illustrated in Theorem 3.3.
Summary: In addition to the summary of our analysis listed in Table (3), we summarize the following dynamics of Model (1) base on mathematical analysis and bifurcation diagrams from our study: 1. The four basic boundary equilibria E 0000 , E K 1 000 , E 00K 2 0, , E K 1 0K 2 0 always exist where E 0000 , E K 1 000 , E 00K 2 0, are always saddle while E K 1 0K 2 0 is locally asymptotically stable if the two inequalities 7 are satisfied. Large dispersal of predators can stabilize the boundary equilibrium E K 1 0K 2 0 when s ∈ (0, 1]. However, the value of s has no effects on the global stability of the boundary equilibrium E K 1 0K 2 0 . (1) can have up to four other boundary equilibria

Model
The number of these boundary equilibria and the stability could be affected by the dispersal strength ρ i , i = 1, 2 and the values of s. For example, the large values of s can destroy these boundary equilibria.
3. In the symmetric case, Model (1) may potentially have three interior equilibria including E = (µ, ν, µ, ν) from Theorem 3.3 when a > d. Although the local stability of E does not depend on s, the large value of s can generate the two additional asymmetric interior equilibria, hence, create multiple interior attractors.

Always exist and always saddle
Always exist and always saddle E K 1 0K 2 0 Always exist; LAS and GAS if µ i > K i for both i = 1, 2 Always exist; GAS if µ i > K i for both i = 1, 2; while LAS if Equations 7 are satisfied Always exist; GAS if Do not exist One or two exist if 3β j µ j +K j < αj < (µj + Kj) 2 with i, j = 1, 2, i = j; Can be locally asymptotically stable or saddle as shown in Figures 3(a), 3(b), 3(c) (2) is satisfied

Do not exist
Do not exist Table 3: Summary of the local and global dynamic of Model (1). LAS refers to the local asymptotical stability, and GAS refers to the global stability.
We implement one and two parameters bifurcation diagrams to obtain insights into the dynamical patterns of the asymmetric two patch Model (1) in the following way: 1. d 1 = 0.85 and a 1 = 1: In the absence of dispersal, the uncoupled two patch model is unstable at the interior equilibrium (5.67, 288.89, 0.33, 80). However, in the presence of the dispersal, Figure 5(a) (blue regions) suggest that the intermediate values of s can stabilize the dynamics while the large values of s with certain dispersal strengths could generate multiple interior equilibria (up to three interior equilibria), thus lead to multiple attractors potentially. Moreover, two dimensional bifurcation diagram shown in Figure 5(b) suggest that the large values of s combined with the small or large dispersal strength ρ 1 in Patch 1 can destroy the interior equilibria (see white regions in Figure 5(b)) with consequences that prey in one patch may go extinct but predator persists in each patch. Table ( a 1 = 1 and d 1 = 0.85 a 1 = 2.1 and d 1 = 2 Scenarios E 1 x 1 y 1 x 2 y 2 E 2 x 1 y 1 x 2 y 2 E 3 x 1 y 1 x 2 y 2 E 1  Figures 5(a), and 6(a). LAS refers to local asymptotical stability, ✗ implies the equilibrium does not exist, and E i x 1 y 1 x 2 y 2 , i = 1, 2, 3 are the three possible interior equilibria of Model (1).
2. d 1 = 2 and a 1 = 2.1: In the absence of dispersal, the uncoupled two patch model has extinction of predator in Patch 1 and is unstable at the boundary equilibrium (10, 0, 0.33, 80). However, in the presence of the dispersal, Figure 6(a) (blue regions) suggest that the intermediate values of s can stabilize the dynamics while the small values of s with certain dispersal strengths could generate multiple interior equilibria (up to three interior equilibria), thus lead to multiple attractors potentially. Moreover, two dimensional bifurcation diagram shown in Figure 6(b) suggest that the large values of s combined with the large dispersal strength ρ 1 in Patch 1 can destroy the interior equilibria (see white regions in Figure 5(b)) with consequences that prey in one patch may go extinct but predator persists in each patch. A more detail dynamic from Figure 6(b) is presented in Table (4). 3. Two parameter bifurcation diagrams of the relative dispersal rate ρ 2 versus the dispersal strategy s for both scenarios of d 1 = 0.85, a 1 = 1 (Figure 7(a)) and d 1 = 2, a 1 = 2.1 (Figure 7(b)). For both cases, the large s combined with the large dispersal strength in Patch 2, i.e., ρ 2 , can destroy the interior equilibrium (see white regions in Figures 7(a) and 7(b) for s > 0.6); while the small s (for d 1 = 0.85, a 1 = 1) and the large value of s (for d 1 = 2, a 1 = 2.1) could generate multiple interior equilibria (see black region for three interior equilibria and red region for two interior equilibria in Figure 7(a) and 7(b)). The proportion of the predators population engaging in the passive dispersal, i.e., s, has profound impacts on the population dynamics of prey and predator presented by Model (1) which generate complicated dynamics including different types of multiple attractors.

Conclusion
We propose and study a two patch prey predator model with the following assumptions: (1) Only predators can migrate and preys are immobile; (2) predators use two dispersal strategies: the passive dispersal and the predation attraction; (3) The model is reduced to the Rosenzweig-MacArthur model in the absence of dispersal. We provide boundedness and positivity of the proposed model in Theorem 2.1. The analytical results which is summarize in Table (3) along with the numerical results presented throughout the paper answer the questions regarding the dynamics of our proposed nonlinear model: When there is no prey in one of the patches, our model applies to the sink-source dynamics where no prey patch is the sink. Analytical results (Theorem 3.1) imply that predators could be driven to extinction locally if the product of the dispersal strength and the proportion of predator population using the passive dispersal (i.e. s) are large. In addition, the sink-source dynamics can process two interior equilibria (see Proposition (3.1)). Our simulations (Figure 1) suggest that the small values of s lead to permanence of the system which is supported by Theorem 3.1. For the intermediate values of s, the system can can have two interiors E l x 1 ,y 1 ,y 2 , l = 1, 2 (i = 1) or E l y 1 ,x 2 ,y 2 , l = 1, 2 (i = 2); For the large values of s, it has no interior equilibria with the consequences that predator goes extinct in two patches. In addition, the intermediate values of s can stabilize the dynamics with certain dispersal strengths (see blue line for locally stable in Figures  1(a), 1(b), and 1(c)).
Theorem (3.2) and Proposition (3.1) provide the existence of the boundary equilibria and the related local stability of our model (1). These results illustrate how s can potentially stabilize the basic boundary equilibria E K 1 0K 2 0 consequently driving predator extinct in both patches locally. Theorem (3.3) provide insights into the existence and stability of a symmetric interior equilibria when Model (1) is symmetric (i.e. in exception of the dispersal strength and dispersal strategy, all life history parameters are the same in both patches). The analytical results indicate that the dispersal strategies do not affect the existence and stability of this symmetric interior equilibria denoted E. However, bifurcation diagrams shown in Figures 4(a) and 4(b) suggest that the large predator population using the passive dispersal could generate two additional asymmetric interior equilibria which can be saddle or locally stable, thus generate bistability between two different attractors (see blue lines in Figure 4(a) when 0.78 ≤ s ≤ 0.92).
Our numerical simulations performed in Section 4 show that the dispersal strategies, i.e., the portion of predator population using the passive dispersal strategies, have huge impacts on the prey and predator populations in two patches. The intermediate predator population using the passive dispersal tends to stabilize the dynamics. Depending on the other life history parameters, the large or small predator population using the passive dispersal with certain dispersal strengths could generate multiple interior equilibria (up to three interior equilibria), thus lead to multiple attractors potentially. When Model (1) has two interior equilibria, it either converges to a boundary attractor or the interior attractor depending initial conditions (see Figures 9(a), and 9(b)); when Model (1) has three interior equilibria, it can have two interior equilibria (see Figures 10(a), and 10(b)). The large predator population using the passive dispersal combined with the large dispersal strength can destroy the interior equilibria with consequences that prey in one patch may go extinct but predator persists in each patch. However, there are situations when the two patch model has no interior equilibrium but all species coexist with fluctuating dynamics (see Figures 8(a) and 8(b)).
The summary of our finding illustrates how population dynamics of prey and predators are affected by changing their foraging behavior. This study give us a better understanding on how combinations of different foraging strategies used by predator favor or affect their coexistence or extinction. Many species tend to adapt to environmental conditions and change their foraging behavior accordingly (see example of foraging behavior of Ants in Markin (1970); Taylor (1977); Traniello et al. (1984)). It will be interesting to look at a two patch prey predator model with adaptive foraging behavior in which adaptation is driven by certain environmental conditions such as temperature or availability of local resources. Such work is on going by the authors.

Appendix: Proofs
Proof of Theorem 2.1 Proof. Observed that dx i dt x i =0 = 0 and dy i dt y i =0 = ρ i sy j ≥ 0 if y j ≥ 0 for i = 1, 2, j = 1, 2, and i = j. The model (1) is positively invariant in R 4 + by theorem A.4 (p. 423) in Thieme (2003). It follows that the set {(x 1 , y 1 , x 2 , y 2 ) ∈ R 4 + : x i = 0} is invariant for both i = 1, 2 under the same theorem. The proof of boundedness is as follow where d min = min{d 1 , d 2 } and T = max This shows that Model (1) is bounded in R 4 + which conclude the proof of theorem (2.1).

Proof of Theorem 3.1
Proof. Item 1: Model (1) is positively invariant and bounded in R 4 + according to Theorem 2.1. From this, it follows that Model (1) is attracted to a compact set C in R 4 + . Furthermore, if x j = 0, j = 1, or 2, then Model (1) is reduced to three species couple models (3). Consider the fact that lim t→∞ y i (t) = lim t→∞ y j (t) = 0 when x i = 0, we can conlcude that y 1 = y 2 = 0 is an omega limit set of Model (3). Additionally then by Theorem 2.5 of Hutson (1984) (Hutson, 1984), prey x i persists.
Item 2: Define V (y i , y j ) = ρ j y i + ρ i y j , then we have Therefore both predators go extinct if µ i > K i . Now Model (3) reduces to the following prey model since lim sup t→∞ y i (t) = lim sup t→∞ y j (t) = 0 Item 3: Now we focus on the persistence condition for predator y i . Since x i is persistent from Item 1 Theorem 3.1 then we can conclude that Model (3) is attracted to a compact set C s subset of C that exclude E 000 .Then according to Theorem 2.1 and 3.2, the omega limit set of Model (3) on the compact set C s is E K i 00 . Notice that the following inequalities, According to Theorem 2.5 of Hutson (1984) (Hutson, 1984), we can conclude that predator y i is persistent if the following inequalities hold holds, then we can conclude that predator y i is persistent. This implies that when time large enough, there exists some ǫ > 0 such that dy i dt y i =0 = ρ j sy j > ρ j sǫ > 0.
Thus, we could conclude that predator in Patch j also persists due to the persistence of predator in Patch i.

Proof of Proposition 3.1
Proof. The algebraic calculations imply that an interior equilibrium (x * i , y * i , y * j ) of Model (3) satisfies the following equations: This implies that Therefore, if a i < d i or µ i > K i or f i (x * i ) = 0 has no positive roots, then Model (3) has no interior equilibrium.
This indicates that there exist x 0 ∈ (−∞, 0) such that f i (x 0 ) = 0. Therefore, we can conclude that f i (x i ) has at least one negative root and at most two positive roots since f i (x i ) is a polynomial with degree 3. The derivative of f i (x i ) has the following form It follows that f i (x i ) has two positive roots if the following equation is satisfied: Since therefore we can conclude that f i (x i ) has two positive roots when 3β i µ i +K i < α i < (µ i + K i ) 2 . Thus for x * iℓ where ℓ = 1, 2 denote the two positive roots of the nullclines f i (x i ) and i = 1, 2 represent the prey population in patch one and two, we have: iℓ < K i , ℓ = 1, 2. From the arguments above we conclude that Model (3) can have up to two interior equilibria On the other hand, if ∆ i = (µ i + K i ) 2 + 3α i < 0 then f i (x i ) has no positive real roots and hence Model (3) has no interior equilibrium.
Based on the discussion above, we can conclude that the results on the local stability of four boundary equilibria of Theorem 3.2 holds.
Item 1: Let p i (x) = a i x 1+x and q i (x) = r i (K i −x)(1+x) a i K i then we have the following Now consider the following Lyapunov functions V 1 (x 1 , y 1 ) = ρ 2 x 1 K 1 p 1 (ξ) − p 1 (K 1 ) p 1 (ξ) dξ + ρ 2 y 1 and V 2 (x 2 , y 2 ) = ρ 1 x 2 K 2 p 2 (ξ) − p 2 (K 2 ) p 2 (ξ) dξ + ρ 1 y 2 Taking derivative of the functions (12) and (13) with respect to time t yield d dt V1(x1(t), y1(t)) = ρ2 p1(x1) − p1(K1) pi(x1) and d dt V2(x2(t), y2(t)) = ρ1 p2 ( Also, we denote V = V 1 + V 2 and adding (14) and (15), we obtain We observe that the function p i (x i ) increases as x i increases thus p i ( This implies that the expressions ρ 2 [p 1 (x 1 ) − p 1 (K 1 )] q 1 (x 1 ) and ρ 1 [p 2 (x 2 ) − p 2 (K 2 )] q 2 (x 2 ) are both negative for all x i ≥ 0 since all the parameters are assumed to be positive. Also, Assume µ i > K i . This implies that d i a i −d i > K i which is also equivalent to a i K i 1+K i = p i (K i ) < d i . Since p i (K i ) < d i then p i (K i ) − d i < 0. The derivative dV dt is therefore negative which implies that both V 1 and V 2 are Lyapunov functions, and the boundary equilibrium E K 1 0K 2 0 = (K 1 , 0, K 2 , 0) is globally stable when µ i > K i by Theorem 3.2 in Hsu (1978).
Therefore, according to Theorem 2.5 of Hutson (1984), we can conclude that prey population in two patches, i.e., x 1 + x 2 , is persistent. Moreover, if x j = 0, Model (1) is reduced to the subsystem (3) where prey x i is persistent according to Theorem 3.1. Thus, we can conclude prey population in at least one patch is persistent.
Therefore, according to Theorem 2.5 of Hutson (1984) and the proof of Proposition 3.1, we can conclude that predator population in each patch is persistent.