1. Introduction

Spatial heterogeneity, dispersal patterns, and biotic interactions influence the distribution of species within a landscape [4, 47, 60, 68, 70]. Spatial self-organization results from local interactions between organisms and the environment, and emerges at patch-scales [61, 75]. For example, limited dispersal ability 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 prey-predator 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 population-dispersal 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 predator-prey 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:

$ $$\begin{aligned} \frac{du_i}{dt}&=r_i u_i\left(1 - \frac{u_i}{k_i}\right) - \frac{b_i u_i v_i}{1 + b_i h_iu_i} \\ \frac{dv_i}{dt}&= \frac{c_ib_i u_i v_i}{1 + b_ih_iu_i} -\delta_i v_i \label{Onepatch-o} \end{aligned}$$ $

(1)

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 $\delta_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 $\frac{b_j u_j v_j}{1 + b_jh_ju_j}$ in Patch $j$, and the dispersal term of predators from Patch $i$ to Patch $j$ is modeled by $\frac{b_j u_j v_j}{1 + b_jh_ju_j} \rho_{ij} v_i$ which gives the net dispersal of predator at Patch $i$ as follows:

$ $$\begin{aligned}\frac{b_j u_j v_j}{1 + b_jh_ju_j} \rho_{ij} v_i-\frac{b_iu_iv_i}{1 + b_ih_iu_i} \rho_{ji} v_j\label{dispersal} \end{aligned}$$ $

(2)

where $\rho_{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 $\rho_{ij}$ is directional and not sysmmetrical, i.e., $\rho_{ij}$ may not be equal to $\rho_{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:

$ $$\begin{aligned} \frac{du_1}{dt}&=r_1 u_1\left(1 - \frac{u_1}{k_1}\right) - \frac{b_1 u_1 v_1}{1 + b_1 h_1u_1} \\ \frac{dv_1}{dt}&= \frac{c_1b_1 u_1v_1}{1 + b_1h_1u_1} -\delta_1 v_1 + \left(\frac{b_1 u_1v_1}{1 + b_1h_1u_1} \rho_{21} v_2-\frac{b_2 u_2v_2}{1 + b_2h_2u_2} \rho_{12} v_1 \right) \\ \frac{du_2}{dt}&=r_2 u_2\left(1 - \frac{u_2}{k_2}\right) - \frac{b_2 u_2 v_2}{1 + b_2 h_2u_2} \\ \frac{dv_2}{dt}&= \frac{c_2b_2 u_2 v_2}{1 + b_2h_2u_2} -\delta_2 v_2+ \left(\frac{b_2 u_2v_2}{1 + b_2h_2u_2} \rho_{12} v_1-\frac{b_1 u_1v_1}{1 + b_1h_1u_1} \rho_{21} v_2 \right). \label{Twopatch-o} \end{aligned}$$ $

(3)

We use the similar rescaling approach in [51] by letting

$x_i=b_ih_i u_i,\,\, y_i=\frac{b_ih_i}{c_i}v_i, \,\, K_i=b_ih_ik_i,\,\, a_i=\frac{c_i}{h_i},\,\, d_i=\frac{c_i\delta_i}{b_ih_i},$

then we have

$ $$\begin{aligned}&v_1\left(\frac{b_1 u_1}{1 + b_1h_1u_1} \rho_{21} v_2-\frac{b_2 u_2v_2}{1 + b_2h_2u_2} \rho_{12} \right)\\=&v_1\left(\frac{b_1 \frac{x_1}{b_1h_1}}{1 + x_1} \rho_{21} \frac{c_2}{b_2h_2}y_2-\frac{b_2 \frac{x_2}{b_2h_2}\frac{c_2}{b_2h_2}y_2}{1 + x_2} \rho_{12} \right)\\ &=v_1\left(\frac{\frac{\rho_{21}a_1a_2 }{c_1b_2}x_1y_2}{1 + x_1}-\frac{\frac{\rho_{12}a_2^2}{c_2b_2}x_2y_2}{1 + x_2} \right) \end{aligned}$$ $

and

$ $$\begin{aligned}&v_2\left(\frac{b_2 u_2}{1 + b_2h_2u_2} \rho_{12} v_1-\frac{b_1 u_1v_1}{1 + b_1h_1u_1} \rho_{21} \right)\\=&v_2\left(\frac{b_2 \frac{x_2}{b_2h_2}}{1 + x_2} \rho_{12} \frac{c_1}{b_1h_1}y_1-\frac{b_1 \frac{x_1}{b_1h_1}\frac{c_1}{b_1h_1}y_1}{1 + x_1} \rho_{21} \right)\\ &=v_2\left(\frac{\frac{\rho_{12}a_1a_2 }{c_2b_1}x_2y_1}{1 + x_2}-\frac{\frac{\rho_{21}a_1^2}{c_1b_1}x_1y_1}{1 + x_1} \right) \end{aligned}$$ $

For mathematical convenience, we assume that $\frac{\rho_{12}}{c_2}=\frac{\rho_{21}}{c_1}=\rho$, then we have

$v_1\left(\frac{\frac{\rho_{21}a_1a_2 }{c_1b_2}x_1y_2}{1 + x_1}-\frac{\frac{\rho_{12}a_2^2}{c_2b_2}x_2y_2}{1 + x_2} \right)=v_1 y_2\rho\frac{a_2 }{b_2}\left( \frac{a_1x_1}{1 + x_1}-\frac{a_2x_2}{1 + x_2}\right)$

and

$v_2\left(\frac{\frac{\rho_{12}a_1a_2 }{c_2b_1}x_2y_1}{1 + x_2}-\frac{\frac{\rho_{21}a_1^2}{c_1b_1}x_1y_1}{1 + x_1} \right)=y_1v_2\rho\frac{a_1 }{b_1}\left(\frac{a_2x_2}{1 + x_2}-\frac{a_1x_1y_1}{1 + x_1} \right).$

Denote that $\rho_1=\rho\frac{a_2 }{b_2}$ and $ \rho_2=\rho\frac{a_1 }{b_1}$, then the two-patch model (3) could be rewritten as the following equations:

$ $$\begin{aligned} \frac{dx_1}{dt}&= r_1x_1\left(1 - \frac{x_1}{K_1}\right) - \frac{a_1 x_1 y_1}{1 + x_1} \\ \frac{dy_1}{dt}&= \frac{a_1 x_1 y_1}{1 + x_1} -d_1 y_1\\ &+ \rho_1\left( \underbrace{ \frac{a_1 x_1 y_1}{1 + x_1}}_\text{attraction strength to Patch 1} y_2- \underbrace{\frac{a_2 x_2 y_2}{1 + x_2}}_\text{attraction strength to Patch 2} \rho_{12}y_1\right)\\ \frac{dx_2}{dt}&= r_2 x_2\left(1 - \frac{x_2}{K_2}\right) - \frac{a_2 x_2 y_2}{1 + x_2} \\ \frac{dy_2}{dt}&= \frac{a_2 x_2 y_2}{1 + x_2} -d_2 y_2+\rho_2\left( \underbrace{\frac{a_2 x_2 y_2}{1 + x_2}}_\text{attraction strength to Patch 2}y_1- \underbrace{\frac{a_1 x_1 y_1}{1 + x_1}}_\text{attraction strength to Patch 1} y_2\right) \label{2patch} \end{aligned}$$ $

(4)

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$; $\rho_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\rightarrow \frac{t}{r_1}$. Thus, for 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

$ $$\begin{aligned} \frac{dx_i}{dt}&=r_i x_i\left(1 - \frac{x_i}{K_i}\right) - \frac{a_i x_i y_i}{1 + x_i} \\ \frac{dy_i}{dt}&= \frac{a_i x_i y_i}{1 + x_i} -d_i y_i \label{Onepatch} \end{aligned}$$ $

(5)

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 Ⅱ functional response.

2.There is no dispersal in prey species. 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.

3.The dispersal of predator from Patch $i$ to Patch $j$ is driven by the prey-predation interaction strength in Patch $j$ termed as the attraction strength.

---

3. Mathematical analysis

The state space of Model (4) is $\mathbb R^4_+$. Let $\mu_i=\frac{d_i}{a_i-d_i},\nu_1=\frac{(K_1-\mu_1)(1+\mu_1)}{a_1K_1}, \nu_2=\frac{r(K_2-\mu_2)(1+\mu_2)}{a_2K_2}$. 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 $\mathbb R^4_+$ with $\limsup_{t\rightarrow\infty} x_i(t)\leq K_i$ for both $i=1,2.$ In addition, it has the following properties:

1.If there is no dispersal in predator, i.e., $\rho_1=\rho_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.

(b) If $\mu_i>K_i$ or $\mu_i<0$, then the boundary equilibrium $(K_i,0)$ is globally asymptotically stable.

(c) If $\frac{K_i-1}{2}<\mu_i<K_i$, then the boundary equilibrium $(K_i,0)$ is a saddle while the interior equilibrium $(\mu_i, \nu_i)$ is globally asymptotically stable.

(d) If $0<\mu_i<\frac{K_i-1}{2}$, then the boundary equilibrium $(K_i,0)$ is a saddle; the interior equilibrium $(\mu_i, \nu_i)$ is a source, and the system has a unique limit cycle which is globally asymptotically stable. In addition, the Hopf bifurcation occurs at $\mu_i=\frac{K_i-1}{2}$.

2.The sets $\{(x_1,y_1,x_2,y_2)\in\mathbb R^4_+:x_i=0 \}$ and $\{(x_1,y_1,x_2,y_2)\in\mathbb 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:

$ $$\begin{aligned} \frac{dx_i}{dt}&=r_i x_i\left(1 - \frac{x_i}{K_i}\right) - \frac{a_i x_i y_i}{1 + x_i} \\ \frac{dy_i}{dt}&= \frac{a_i x_i y_i}{1 + x_i} -d_i y_i \end{aligned}$$ $

(6)

$ $$\begin{aligned}\frac{dx_j}{dt}&=r_j x_j\left(1 - \frac{x_j}{K_j}\right) \label{yj0} \end{aligned}$$ $

(7)

where $\lim_{t\rightarrow\infty}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<\mu_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

$\mu_i=\frac{d_i}{a_i-d_i},\nu_1=\frac{(K_1-\mu_1)(1+\mu_1)}{a_1K_1}, \nu_2=\frac{r(K_2-\mu_2)(1+\mu_2)}{a_2K_2}.$

We define the following notations for all possible boundary equilibria of Model (4):

$ $$\begin{array}{cccc} E_{0000}=(0,0,0,0), &E_{K_1000}=(K_1,0,0,0),&E_{\mu_1\nu_100}=(\mu_1,\nu_1,0,0),\\ E_{K_10\mu_2\nu_2}=(K_1,0,\mu_2,\nu_2),&E_{K_10K_20}=(K_1,0,K_2,0), &E_{00K_20}=(0,0,K_2,0),\\ E_{00\mu_2\nu_2}=(0,0,\mu_2,\nu_2),&E_{\mu_1\nu_1K_20}=(\mu_1,\nu_1,K_2,0) \end{array}$$ .$

The following theorem provides sufficient conditions on the existence and stability of these boundary equilibria:

Theorem 3.2.  [Boundary equilibria of Model (4)] Model (4) always has the following four boundary equilibria

$E_{0000},\,E_{K_1000},\,E_{00K_20},\,E_{K_10K_20}$

where the first three ones are saddles while $E_{K_10K_20}$ is locally asymptotically stable if $\mu_i>K_i$ and it is a saddle if $\left(\mu_1-K_1\right)\left(\mu_2-K_2\right)<0$ or $\mu_i<K_i, i=1,2$. Let $i,j=1,2,\,i\neq j$, and

$E_{11}^b=E_{\mu_1\nu_100},\,E_{12}^b=E_{\mu_1\nu_1K_20},\,E_{21}^b=E_{00\mu_2\nu_2} \;and\;E_{22}^b=E_{K_10\mu_2\nu_2}.$

Then if $0<\mu_i<K_i$, then Model (4) has additional two boundary equilibria $E^b_{i1}$ and $E_{i2}^b$ where $E_{i1}^b$ is always a saddle. The boundary equilibrium $E_{i2}^b$ is locally asymptotically stable if $\frac{K_i-1}{2}<\mu_i<K_i$ and one of the following conditions holds:

sa: $a_j\leq d_i, K_j<\mu_j$.

sb: $0<K_j<\min\Big\{\mu_j,\frac{d_i}{a_j-d_i}\Big\}.$

sc: $0<\frac{d_i}{a_j-d_i}<K_j<\mu_j \; and\;\rho_j<\frac{d_j-K_j(a_j-d_j)}{\nu_i\left[K_j(a_j-d_i)-d_i\right]}.$

sd: $0<\mu_j<K_j<\frac{d_i}{a_j-d_i} \; and \;\rho_j>\frac{K_j(a_j-d_j)-d_j}{\nu_i\left[d_i-K_j(a_j-d_i)\right]}.$

And $E_{i2}^b$ is a saddle if $0<\mu_i<\frac{K_i-1}{2}$ or one of the following conditions holds:

ua: $K_j>\max\Big\{\mu_j,\frac{d_i}{a_j-d_i}\Big\}.$

ub: $0<\frac{d_i}{a_j-d_i}<K_j<\mu_j\; and \;\rho_j>\frac{d_j-K_j(a_j-d_j)}{\nu_i\left[K_j(a_j-d_i)-d_i\right]}.$

uc: $\mu_j<K_j<\frac{d_i}{a_j-d_i} \; and\; \rho_j<\frac{K_j(a_j-d_j)-d_j}{\nu_i\left[d_i-K_j(a_j-d_i)\right]}.$

In addition, if $0<\mu_i<K_i$ for both $i=1$ and $i=2$, the boundary equilibria $E_{12}^b$ and $E_{22}^b$ 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_{12}^b$ and $E_{22}^b$ 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_{K_10K_20}$.

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<\mu_j<K_j$. However, in the presence of dispersal, large values of $\rho_j$ can lead to the local stability of the boundary equilibrium $E^b_{i2}$ where $i,j=1,2$ and $i\neq j$, under conditions of $\mu_j<K_j<\frac{d_i}{a_j-d_i}$.

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<\mu_j$. However, in the presence of dispersal, large values of $\rho_j$ can drive the boundary equilibrium $E^b_{j2}$ being unstable, under conditions of $0<\frac{d_i}{a_j-d_1}<K_j<\mu_j$.

4.Under conditions of $\mu_i<K_i$, the boundary equilibria $E_{12}^b$ and $E_{22}^b$ 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_{K_10K_20}$:

Theorem 3.3.  Model (4) has global stability at $E_{K_10K_20}$ if $\mu_i>K_i, i=1,2$.

Notes and biological implications.  Theorem 3.3 implies that the dispersal of predators does not effect the global stability of the boundary equilibrium $E_{K_10K_20}$.

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 $\mathbb 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

$b\leq\liminf\limits_{\tau\rightarrow\infty} z(\tau)\leq \limsup\limits_{\tau\rightarrow\infty} z(\tau)\leq B.$

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 $\mathbb R^4_+$ if there exists constants $0<b<B$, such that for any initial condition taken in $\mathbb R^4_+$ with $x_1(0)y_1(0)x_2(0)y_2(0)>0$, the following inequality holds

$b\leq\liminf\limits_{\tau\rightarrow\infty}\min\{ x_1(\tau), y_1(\tau), x_2(\tau),y_2(\tau)\}$

$ \leq \limsup\limits_{\tau\rightarrow\infty} \max\{x_1(\tau), y_1(\tau), x_2(\tau),y_2(\tau)\}\leq B. $

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

1.$\mu_j<K_j, \mu_i>K_i$. Or

2.$\frac{K_i-1}{2}<\mu_i<K_i \;and\; K_j>\max\Big\{\mu_j,\frac{d_i}{a_j-d_i}\Big\}$. Or

3.$\frac{K_i-1}{2}<\mu_i<K_i,\,\,\mu_j<K_j<\frac{d_i}{a_j-d_i} \;and\;\rho_j<\frac{K_j(a_j-d_j)-d_j}{\nu_i\left[d_i-K_j(a_j-d_i)\right]}$.

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 $\frac{K_i-1}{2}<\mu_i<K_i,\,\,\mu_j<K_j<\frac{d_i}{a_j-d_i}$, can keep the persistence of predator $j$. This is consistent with the results of Item uc in Theorem 3.2.

Theorem 3.7  [Permanence of the two patch dispersal model] Model (4) is permanent if one of the following inequalities hold

1.$\frac{K_j-1}{2}<\mu_j<K_j, \,0<\frac{d_j}{a_i-d_j}<K_i<\mu_i \;and\; \rho_i>\frac{d_i-K_i(a_i-d_i)}{\nu_j\left[K_i(a_i-d_j)-d_j\right]}$ where $i=1,j=2$ or $i=2,j=1$. Or

2.$\frac{K_j-1}{2}<\mu_j<K_j,\, \mu_i>\frac{K_i-1}{2} \;and\; K_i>\max\Big\{\mu_i,\frac{d_j}{a_i-d_j}\Big\}$ for both $i=1,j=2$ and $i=2,j=1$. Or

3.$\frac{K_i-1}{2}<\mu_i<K_i,\,K_i>\max\Big\{\mu_i,\frac{d_j}{a_i-d_j}\Big\},\,\frac{K_j-1}{2}<\mu_j<K_j<\frac{d_i}{a_j-d_i} \;and\;\rho_j<\frac{K_j(a_j-d_j)-d_j}{\nu_i\left[d_i-K_j(a_j-d_i)\right]}$ where $i=1,j=2$ or $i=2,j=1$.

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 $\mu_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}$.

---

3.3. Interior equilibrium and stability

Let $p_i(x)=\frac{a_i x}{1+x}$ and $q_i(x)=\frac{r_i(K_i-x)(1+x)}{a_iK_i}$, then we have

$ $$\begin{array}{lcl} \frac{dx_i}{dt}&=&r_ix_i\left(1-\frac{x_i}{K_i}\right)-\frac{a_i x_iy_i}{(1+x_i)}\\&=&\frac{a_i x_i}{1+x_i}\left[\frac{r_i(K_i-x_i)(1+x_i)}{a_iK_i}-y_i\right]=p_i(x_i)\left[q_i(x_i)-y_i\right]\\\\ \frac{dy_i}{dt}&=&y_i\left[\frac{a_i x_i}{1+x_i}-d_i+\rho_i y_j\left(\frac{a_i x_i}{1+x_i}-\frac{a_j x_j}{1+x_j}\right)\right]\\&=&y_i\left[p_i(x_i)-d_i+\rho_i y_j\left(p_i(x_i)-p_j(x_j)\right)\right]\end{array}$$ $

If $(x_1^*,y_1^*,x_2^*, y_2^*)$ is an interior equilibrium of Model (4), then it satisfies the following equations:

$ \label{interior-eq1} $$\begin{array}{lcl} q_i(x_i)-y_i=0\Leftrightarrow q_i(x_i)=y_i,\,\\ p_i(x_i)-d_i+\rho_i y_j\left(p_i(x_i)-p_j(x_j)\right)=0 \Leftrightarrow p_i(x_i)=\frac{a_i x_i}{1+x_i}\\=\frac{\rho_i y_jp_j(x_j)+d_i}{1+\rho_i y_j}=\frac{\rho_i q_j(x_j)p_j(x_j)+d_i}{1+\rho_i q_j(x_j)} \end{array}$$ $

(8)

which gives:

$ x_i=\frac{\rho_i q_j(x_j)p_j(x_j)+d_i}{a_i(1+\rho_i q_j(x_j))-(\rho_i q_j(x_j)p_j(x_j)+d_i)} $

(9)

$ =\frac{\rho_i q_j(x_j)p_j(x_j)+d_i}{\rho_i q_j(x_j)\left[a_i-p_j(x_j)\right]+a_i-d_i} $

(10)

Since $\limsup_{t\rightarrow\infty} x_i(t)\leq K_i$ for both $i=1,2$ and $y_i=q_i(x_i)$, therefore, positive solutions of $x_i\in (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:

$ \label{interior-eq3} $$\begin{array}{lcl} x_1&=&\frac{a_2\left[r_2\rho_1x_2\left(K_2-x_2\right)+K_2d_1\right]}{r_2\rho_1x_2\left(K_2a_1-K_2a_2-a_1\right)-r_2\rho_1x_2^2(a_1-a_2)+K_2(a_1r_2\rho_1+a_1a_2-a_2d_1)}\\&=&\frac{f_t(x_2)}{f_b(x_2)}=F(x_2)\\\\ x_2&=&\frac{a_1\left[r_1\rho_2x_1\left(K_1-x_1\right)+K_1d_2\right]}{r_1\rho_2x_1\left(K_1a_2-K_1a_1-a_2\right)-r_1\rho_2x_1^2(a_2-a_1)+K_1(a_2r_1\rho_2+a_1a_2-a_1d_2)}\\&=&\frac{g_t(x_1)}{g_b(x_1)}=G(x_1) \end{array}$$ $

(11)

with $r_1=1, r_2=r$ and the following properties:

1.$F(0)=\frac{f_t(0)}{f_b(0)}=\frac{a_2K_2d_1}{K_2(a_1r_2\rho_1+a_1a_2-a_2d_1)}=\frac{a_2d_1}{a_1r_2\rho_1+a_1a_2-a_2d_1}$ and $F(K_2)=\frac{f_t(K_2)}{f_b(K_2)}=\frac{d_1}{a_1-d_1}=\mu_1$.

2.$f_t(x_2)=a_2\left[r_2\rho_1x_2\left(K_2-x_2\right)+K_2d_1\right]\geq a_2K_2 d_1>0\;\mbox{ for }\;x_2\in [0, K_2]$ and $f_b(x_2)\big\vert_{a_1=a_2=a}=a\left[r_2\rho_1(K_2-x_2)+K_2(a-d_1)\right].$

3.$G(0)=\frac{g_t(0)}{g_b(0)}=\frac{a_1K_1d_2}{K_1(a_2r_1\rho_2+a_1a_2-a_1d_2)}=\frac{a_1d_2}{a_2r_1\rho_2+a_1a_2-a_1d_2}$ and $G(K_1)=\frac{g_t(K_1)}{g_b(K_1)}=\frac{d_2}{a_2-d_2}=\mu_2$.

4.$g_t(x_1)=a_1\left[r_1\rho_2x_1\left(K_1-x_1\right)+K_1d_2\right]\geq a_1 K_1d_1>0\mbox{ for } x_1\in [0, K_1]$ and $g_b(x_1)\big\vert_{a_1=a_2=a}=a\left[r_1\rho_2(K_1-x_1)+K_1(a-d_2)\right].$

Theorem 3.8.  [Interior equilibrium] If $\mu_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

$x^c_i=\frac{K_i\left(r_i\rho_j+a_i-d_j-\sqrt{(a_i-d_j)(r_i\rho_j+a_i-d_j)}\right)}{r_i\rho_j}.$

Model (4) has no interior equilibrium if $ a_i>a_j, \,\rho_i<\frac{4K_ja_j(a_i-a_j)(d_i-a_i)}{r_j\left(K_ja_i-K_ja_j+a_i\right)^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$

$a_i>\max\{a_j,d_1,d_2\}, a_j>\max\{d_1,d_2\}, \,\rho_j<\frac{4K_ia_i(a_j-a_i)(d_j-a_j)}{r_i\left(K_ia_j-K_ia_i+a_j\right)^2},$

$ \,F(x^c_2)<K_1,\,\;and\; G(x^c_1)<K_2$

where sufficient conditions for the inequalities $F(x^c_2)<K_1,\,\;and\; G(x^c_1)<K_2$ hold are

$\rho_i\leq \frac{4(K_ia_i-K_id_i-d_i)}{K_jr_j} \;and\;\rho_j<\frac{4K_ja_j(K_ia_i-K_id_i-d_i)}{a_jr_jK_j^2+r_jK_i(K_ja_j-K_ja_i-a_i)^2}.$

In addition, we have $\frac{a_id_j}{a_jr_i\rho_j+a_ia_j-a_id_j}<x_j^*<K_j$ for both $i=1,j=2$ and $i=2,j=1$.

2.Assume that $a_1=a_2=a$. Model (4) has no interior equilibrium if $d_1>a+r_2\rho_1$ or $d_2>a+r_1\rho_2$ while it 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\},\,F(x^c_2)<K_1,\,\;and\; G(x^c_1)<K_2$

where sufficient conditions for the inequalities $F(x^c_2)<K_1,\,\;and\; G(x^c_1)<K_2$ hold are

$\rho_i<\frac{4(K_ia-K_id_i-d_i)}{K_jr_j}$

for both $i=1,j=2$ and $i=2,j=1$. In addition, we have $\frac{d_j}{r_i\rho_j+a-d_j}<x_j^*<K_j$ for both $i=1,j=2$ and $i=2,j=1$.

Notes and biological implications. Theorem 3.8 provides sufficient conditions on the existence of no interior equilibrium when $\mu_i>K_i$ for either $i=1$ or $i=2$; and at least one interior equilibrium of Model (4) when $\mu_i<K_i$ for both $i=1,2$. The results indicate follows:

1.If $\mu_i>K_i$, then Model (4) has no interior equilibrium if the dispersal of its predator is too small.

2.If $\mu_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).

Theorem 3.9.   [Interior equilibrium and the stability] Suppose that $d_1=d_2=d$. Let

$\mu_1=\frac{d}{a_1-d}, \nu_1=\frac{(K_1-\mu_1)(1+\mu_1)}{a_1K_1}, \mu_2=\frac{d}{a_2-d}, \nu_2=\frac{r(K_2-\mu_2)(1+\mu_2)}{a_2K_2}.$

If $0<\mu_i<K_i$ for both $i=1$ and $i=2$, then $E^i=(\mu_1,\nu_1,\mu_2,\nu_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 $\frac{K_i-1}{2}<\mu_i<K_i$ hold for both $i=1$ and $i=2$ while it is unstable if the following inequality holds

$\frac{\mu_1(K_1a_1-a_1-K_1d-d)}{a_1K_1}+\frac{r\mu_2(K_2a_2-a_2-K_2d-d)}{a_2K_2}>0.$

2.Assume that $(K_1a_1-a_1-K_1d-d)(K_2a_2-a_2-K_2d-d)<0$ and

$\frac{\mu_1(K_1a_1-a_1-K_1d-d)}{a_1K_1}+\frac{r\mu_2(K_2a_2-a_2-K_2d-d)}{a_2K_2}<0.$

If $K_ia_i-a_i-K_id-d>0$ (i.e., $\mu_i<\frac{K_i-1}{2}$) for either $i=1$ or $i=2$, then the large values of $\rho_i$ can make $E^i$ being locally asymptotically stable, i.e.,

$ \rho_i>\max\Big\{\frac{-\nu_j-r_j\mu_i\mu_j(K_ia_i-K_id-a_i-d)(K_ja_j-K_jd-a_j-d)}{(K_iK_ja_j\nu_jd\nu_i(a_i-d)^2)},$

$\frac{-\frac{\mu_i\nu_jK_j(\nu_i\rho_j+1)(a_j-d)^2(K_ia_i-K_id-a_i-d)}{r_j\mu_j\nu_iK_i(a_i-d)^2(K_ja_j-K_jd-a_j-d)}-1}{\nu_j}\Big\} .$

If, in addition, $a_1=a_2=a, K_1=K_2=K, r=1$, then $\mu_1=\mu_2=\mu=\frac{d}{a-d},\,\nu_1=\nu_2=\nu=\frac{(K-\mu)(1+\mu)}{aK}$, and $E^i=(\mu,\nu,\mu,\nu)$ is the only interior equilibrium for Model (4) which has the same local stability as the interior equilibrium $(\mu,\nu)$ for the single patch model (5), i.e., $E^i$ is locally asymptotically stable if $\frac{K-1}{2}<\mu<K$ while it is unstable if $\mu<\frac{K-1}{2}$.

Notes and biological implications.  Theorem 3.9 implies Model (4) has an interior equilibrium $E^i=(\mu_1,\nu_1,\mu_2,\nu_2)$ if $d_1=d_2=d$ and $0<\mu_i<K_i$ for both $i=1,2$. In addition, Theorem 3.9 indicates that dispersal of predators has no effects on the local stability if $\frac{K_i-1}{2}<\mu_i<K_i$ for both $i=1,2$ or one of the single patch models (5) is unstable and $\frac{\mu_1(K_1a_1-a_1-K_1d-d)}{a_1K_1}+\frac{r\mu_2(K_2a_2-a_2-K_2d-d)}{a_2K_2}>0.$ However, large dispersal of predator at Patch $i$ can stabilize the interior equilibrium when its single patch model model is unstable at $(\mu_i,\nu_i)$ with $\frac{\mu_1(K_1a_1-a_1-K_1d-d)}{a_1K_1}+\frac{r\mu_2(K_2a_2-a_2-K_2d-d)}{a_2K_2}<0.$

---

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 $x_i^*=K_i, y_i^*=0, x_j^*=\mu_j, y_j^*=\nu_j$ depending on additional conditions.

2.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 $\mu_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_{K_10K_20}$ and the symmetric interior equilibrium $(\mu,\nu,\mu,\nu)$ when it exists.

4.If $d_i>a_i$, then small dispersal of predator at Patch $i$ prevents the interior equilibrium while if $0<\mu_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.

5.If $d_i=d_j$ and $0<\mu_i<\frac{K_i-1}{2}, \frac{K_j-1}{2}<\mu_j<K_j$, then large dispersal of predator at Patch $i$ can stabilize the interior equilibrium $(\mu_i,\nu_i,\mu_j,\nu_j)$.

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).

---

4.1. 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 $\left(\frac{0.2}{a_1-0.2},\frac{\left(5-\frac{0.2}{a_1-0.2}\right)\left(1+\frac{0.2}{a_1-0.2}\right)}{5}\right)$ if $2<\frac{0.2}{a_1-0.2}<5\Leftrightarrow 0.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 $\left(\frac{0.1}{a_2-0.1},\frac{1.5\left(3-\frac{0.1}{a_2-0.1}\right)\left(1+\frac{0.1}{a_2-0.1}\right)}{3}\right)$ if $1<\frac{0.1}{a_2-0.1}<3\Leftrightarrow 0.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)$, $\left(2,1.5\right)$, respectively. After turning on the dispersal, the coupled two patch model can have one interior equilibrium (see the blue regions in Figure (1a)) which can be locally stable (see the blue dots in Figure (1b)), or be a saddle (see the green dots in Figure (1b)) where the coupled system has fluctuated dynamics; or it can have two interior equilibria (see the red regions in Figure (1a)) which could be saddles and generate bistability between fluctuated interior dynamics and the boundary attractor $(4,1,3,0)$ (see the examples of two interior saddles of Figure (4a)); or it can have three interior equilibria (see the black regions in Figure (1a)) which generate multiple interior attractors (see the examples of two saddles, one sink and two sinks, one saddle of Figure (1b)) or it could have no interior equilibrium (see white and yellow regions of Figure (1a)). Bifurcation diagrams Figure (1a), (1b), and (4b) 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 $\left(2,1.5\right)$. After turning on the dispersal, the coupled two patch model can have one interior equilibrium (see the blue regions in Figure (2a)) which can be locally stable (see the blue dots in Figure (2b)), or be a saddle (see the green dots in Figure (2b)), or be a source (see the red dots in Figure (2b)) 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 (1a)) which could be two saddles or one sink, one saddle and generate bistability between the interior attractor and the boundary attractor (see Figure (2b)); or it can have three interior equilibria (see the black regions in Figure (1a)) which generate multiple interior attractors (see the examples of two sinks and one saddle of Figure (4b)) or it could have no interior equilibrium (see white regions of Figure (2a)). Bifurcation diagrams of Figure (2a), (2b), and (4b) 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 (3a)) which can be a sink (see the red dots in Figure (3b)) where the coupled system has fluctuated dynamics; or it can have two interior equilibria (see the red regions in Figure (2a)) which could be two saddles, one sink v.s. one saddle, one source v.s. one saddle and generate bistability between the interior attractors and the boundary attractor (see Figure (3b)); it could have no interior equilibrium (see white and yellow regions of Figure (3a)). Bifurcation diagrams Figure (2a)-(3b) 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 (1b)-(4b)).

• 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 (1b), (4a)).

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) converges either to the interior sinks or to the boundary attractors with one predator going extinct for almost all initial conditions (see examples in Figure (2b), (3b), (4b)).

• two interior saddles: Model (4) converges either to the fluctuated interior attractors or to the boundary attractors with one predator going extinct for almost all initial conditions depending on the initial conditions with $x_1(0)y_1(0)x_2(0)y_2(0)>0$ (see examples in Figure (2b), (3b), 4).

• one interior source 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 either to the fluctuated interior attractors or to the boundary attractors with one predator going extinct for almost all initial conditions (see examples in Figure (3b)).

Model (4) has bistability between interior attractors and the boundary attractors whenever it processes two interior equilibria. This implies that depending 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]:

$ $$\begin{aligned} \frac{dx_i}{dt}&=r_i x_i\left(1 - \frac{x_i}{K_i}\right) - \frac{a_i x_i y_i}{1 + x_i} \\ \frac{dy_i}{dt}&= \frac{a_i x_i y_i}{1 + x_i} -d_i y_i+ \rho_i( y_j-y_i)\\ \frac{dx_j}{dt}&=r_j x_j\left(1 - \frac{x_j}{K_j}\right) - \frac{a_jx_jy_j}{1 + x_j} \\ \frac{dy_j}{dt}&= \frac{a_j x_j y_j}{1 + x_j} -d_j y_j- \rho_j( y_j-y_i) \label{2patch2} \end{aligned}$$ $

(12)

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 $\rho_i=\rho_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: (ⅰ) for very small migration rates the oscillations always synchronize; (ⅱ) For intermediate migration rates the synchronous oscillations are unstable and there are periodic, quasi-periodic, and intermittently chaotic attractors with asynchronous dynamics; and (ⅱ) 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:

Theorem 4.1.  [Summary of the dynamics of Model (12)] Define

$\hat{\mu_i}=\frac{\hat{d_i}}{a_i-\hat{d_i}},\, \hat{\nu_i}=q_i(\hat{\mu_i})=\frac{r_i(K_i-\hat{\mu_i})(1+\hat{\mu_i})}{a_iK_i} \;and\; \hat{\nu}_j^i=\frac{\rho_j\hat{\nu_i}}{d_j+\rho_j}$

where $\hat{d_i}=d_i +\frac{\rho_id_j}{d_j+\rho_j}$. Let $E^b_1=E_{x_1y_10y_2}=\left(\hat{\mu_1},\hat{\nu_1},0,\hat{\nu}_2^1\right)$ and $E^b_2=E_{0y_1x_2y_2}=\left(0,\hat{\nu}^2_1,\hat{\mu_2},\hat{\nu}_2\right)$. Then we have the following summary on the dynamics of Model (12)

1.Model (12) is positively invariant and bounded in its state space $\mathbb R^4_+$ with $\limsup_{t\rightarrow\infty} x_i(t)\leq K_i$ for both $i=1$ and $i=2$.

2.Boundary equilibria: Model (12) always has the following four boundary equilibria

$E_{0000},\,E_{K_1000},\,E_{00K_20},\,E_{K_10K_20}$

where the first three ones are saddles while $E_{K_10K_20}$ is locally asymptotically stable if

$d_1+d_2+\rho_1+\rho_2>\frac{a_1K_1}{1+K_1}+\frac{a_2K_2}{1+K_2}$

and

$ $$\begin{array}{l} &\hat{d}_1-\frac{a_1K_1}{1+K_1} +\frac{K_2a_2\left(\frac{a_1K_1}{1+K_1}-d_1-\rho_1\right)}{(d_2+\rho_2)(1+K_2)}>0\\ &\Leftrightarrow&\\ & \left[d_1-\frac{a_1K_1}{1+K_1}\right]\left[1-\frac{K_2a_2}{(d_2+\rho_2)(1+K_2)}\right]+\frac{\rho_1}{d_2+\rho_2}\left[d_2-\frac{K_2a_2}{(1+K_2)}\right]>0.\end{array}$$ $

and it is a saddle if one of the above inequalities does not hold. If $0<\hat{\mu}_i<K_i$, then the boundary equilibrium $E^b_i$ exists which is locally asymptotically stable if $\frac{K_i-1}{2}<\hat{\mu_i}<K_i, r_j<a_j \hat{\nu}_j^i$.

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$:

$ $$\begin{aligned} \frac{dx_i}{dt}&=r_i x_i\left(1 - \frac{x_i}{K_i}\right) - \frac{a_i x_i y_i}{1 + x_i} \\ \frac{dy_i}{dt}&= \frac{a_i x_i y_i}{1 + x_i} -d_i y_i+ \rho_i( y_j-y_i)\\ \frac{dy_j}{dt}&= -d_j y_j- \rho_j( y_j-y_i) \label{2patch2s} \end{aligned}$$ $

(13)

whose global dynamics can be described as follows:

3a Prey $x_i$ is persistent for Model (13) with $\limsup_{t\rightarrow\infty} x_i(t)\leq K_i$.

3b Model (13) has global stability at $(K_i,0,0)$ if $\hat{\mu_i}>K_i$.

3c Model (13) has global stability at $(\hat{\mu_i},\hat{\nu_i},\hat{\nu}_j^i)$ if $\frac{K_i-1}{2}<\hat{\mu_i}<K_i.$

4.Persistence of prey: Prey $x_i$ persists if $\hat{\mu_j}<0$, or $\hat{\mu_j}>K_j$, or $\frac{K_j-1}{2}<\hat{\mu_j}<K_j, r_i>a_i \hat{\nu}_i^j$ hold.Both prey $x_i$ and $x_j$ persist if one of the following three conditions hold

4(a) $\hat{\mu_i}>K_i$ for both $i=1$ and $i=2$. Or $\hat{\mu_i}<0$ for both $i=1$ and $i=2$.

4(b)$\frac{K_i-1}{2}<\hat{\mu_i}<K_i, r_j>a_j \hat{\nu}_j^i$ for both $i=1,j=2$ and $i=2,j=1$.

4(c) $\hat{\mu_i}>K_i$, or $\hat{\mu_i}<0$, and $\frac{K_j-1}{2}<\hat{\mu_j}<K_j, r_i>a_i\hat{\nu}_i^j$ for either $i=1, j=2$ or $i=2,j=1$.

5.Extinction of prey $x_i$: Prey $x_i$ goes to extinction if $\frac{K_j-1}{2}<\hat{\mu_j}<K_j$ and $\frac{r_i (K_i+1)^2}{4a_iK_i}<\hat{\nu}_j^i$. In addition, under these conditions, Model (12) has global stability at $(0,\hat{\nu}_i^j,\hat{\mu_j},\hat{\nu_j})$.

6.Persistence and extinction of predators: Predator $y_i$ and $y_j$ have the same persistence and extinction conditions. Predators persist if $0<\mu_i<K_i$ for both $i=1$ and $i=2$ while both predators go extinct if $\mu_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<\mu_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 $\mu_1=\mu_2=\mu$ and $\nu_1=\nu_2=\nu$. Therefore, we can conclude that Model (13) has global stability at $(\mu,\nu,\mu,\nu)$ if $\frac{K-1}{2}<\mu<K$. In addition, the local stability of $(\mu,\nu,\mu,\nu)$ for Model (12) is the same as the local stability of $(\mu,\nu)$ for single patch models when $\rho_1=\rho_2=0$ of Model (12).

Notes. Theorem 4.1 indicates follows:

1.If $\mu_i>K_i$ and $0<\mu_j<K_j$, then the large dispersal of predator at Patch $i$ stabilizes $E_{K_10K_20}$.

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} = \left(0,0,0,0\right)$, $E_{K_1000} = \left( K_1,0,0,0\right)$, $E_{00K_20} = \left( 0,0,K_2,0\right)$ and $E_{K_10K_20} = \left( K_1,0,K_2,0\right)$ among which first three are saddles in both the cases. If $0<\mu_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<\hat{\mu}_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_{K_10K_20}$, $E_{i2}^b$ and $E_i^b$. 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 Table 2, we can conclude that the strength of dispersal ability of predator has 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 (ⅰ) 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; (ⅱ) 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 prey-contact 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 prey-predator 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 (1a), (2a), and (3a)).

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 (1a), (2a), and (3a)).

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 (1b), (2b), and (2b)) 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.

---

6. Proofs

Proof of Theorem 3.1

Proof.Notice that both $\frac{dx_i}{dt} \big\vert_{x_i=0}=0$ and $\frac{dy_i}{dt} \big\vert_{y_i=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 $\mathbb 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:

$\frac{dx_i}{dt}= r_i x_i\left(1 - \frac{x_i}{K_i}\right) - \frac{a_i x_i y_i}{1 + x_i}\leq r_i x_i\left(1 - \frac{x_i}{K_i}\right)$

which implies that

$\limsup\limits_{t\rightarrow\infty} x_i(t)\leq K_i.$

Define $V=\rho_2 (x_1+y_1)+\rho_1 (x_2+y_2)$, then we have

$ $$\begin{array}{l} \;\;\;\frac{dV}{dt}\\=rho_2\frac{d(x_1+y_1)}{dt}+\rho_1\frac{d(x_2+y_2)}{dt}\\ =\rho_2x_1\left(1 - \frac{x_1}{K_1}\right) +\rho_1 rx_2\left(1 - \frac{x_2}{K_2}\right) -\rho_2 d_1 y_1 -\rho_1 d_2 y_2\\ =\rho_2x_1\left(1 - \frac{x_1}{K_1}\right) +\rho_1 rx_2\left(1 - \frac{x_2}{K_2}\right)+\rho_2 d_1 x_1 +\rho_1 d_2 x_2 \\\;\;\;-\rho_2 d_1 (x_1+y_1) -\rho_1 d_2 (x_2+y_2)\\ \leq M-d \left[\rho_2 (x_1+y_1)+\rho_1 (x_2+y_2)\right]=M-d V \end{array}$$ $

where $d=\min\{d_1,d_2\}$ and

$M=\max\limits_{0\leq x_1\leq K_1}\Big\{\rho_2x_1\left(1 - \frac{x_1}{K_1}+ d_1 \right)\Big\} +\max\limits_{0\leq x_2\leq K_2}\Big\{\rho_1x_2\left(1 - \frac{x_2}{K_2}+d_2 \right)\Big\}.$

Therefore, we have

$\limsup\limits_{t\rightarrow\infty} V(t)=\limsup\limits_{t\rightarrow\infty} \rho_2 (x_1(t)+y_1(t))+\rho_1 (x_2(t)+y_2(t))\leq \frac{M}{d}$

which implies that Model (5) is bounded in $\mathbb R^4_+$.

If there is no dispersal in predator, i.e., $\rho_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.

Recall that both $\frac{dx_i}{dt} \big\vert_{x_i=0}=0$ and $\frac{dy_i}{dt} \big\vert_{y_i=0}=0$ for $i=1,2$, therefore, the sets $\{(x_1,y_1,x_2,y_2)\in\mathbb R^4_+:x_i=0 \}$ and $\{(x_1,y_1,x_2,y_2)\in\mathbb R^4_+: y_i=0 \}$ are invariant. This indicates that if $x_j(0)=0$, then $x_j(t)=0$ for all $t>0$. Therefore, the population of $y_j$ converges to 0 since

$\frac{dy_j}{dt} = -d_j y_j-\rho_j\frac{a_i x_i y_i}{1 + x_i} y_j \leq 0\Rightarrow \limsup\limits_{t\rightarrow\infty} y_j(t)=0.$

Applying the results in [69], we can conclude that Model (4) is reduced to the single patch model (5) when $x_j=0$.

In the case that $y_j=0$, Model (4) is reduced to Model (7) by replacing $y_j=0$ in Model (4).

Summarizing the discussions above, we can conclude that the statement of Theorem 3.1 holds.

Proof of Theorem 3.2.

Proof.According Theorem 3.1, sufficient condition for the single patch model (5) having the unique interior equilibrium $(\mu_i,\nu_i),i=1,2$ is $\mu_i<K_i$. Therefore, sufficient condition for Model (4) having boundary equilibria $E_{\mu_1\nu_100}$ and$E_{\mu_1\nu_1K_20}$ is $\mu_1<K_1$. Similarly, sufficient condition for Model (4) having boundary equilibria $E_{00\mu_2\nu_2}$ and $E_{K_10\mu_2\nu_2}$ is $\mu_2<K_2$.

The local stability of an equilibrium $(x_1^*,y_1^*,x_2^*,y_2^*)$ can be determined by the eigenvalues $\lambda_i$, $i = 1,~2,~3,~4$ of the Jacobian matrix $J_{(x_1^*,y_1^*,x_2^*,y_2^*)}$ (14) of Model (4) evaluated at the equilibrium.

$ $$\begin{array}{l} {J_{(x_1^*,y_1^*,x_2^*,y_2^*)}} = \\ \left| \begin{array}{l} \left( {1 - \frac{{2x_1^*}}{{{K_1}}}} \right) - \frac{{{a_1}y_1^*}}{{{{\left( {1 + x_1^*} \right)}^2}}}\;\;\;\;\;\;\;\;\;\;\; - \frac{{{a_1}x_1^*}}{{1 + x_1^*}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;0\\ \;\;\;\;\;\;\frac{{{a_1}y_1^*\left( {1 + {\rho _1}y_2^*} \right)}}{{{{\left( {1 + x_1^*} \right)}^2}}}\;\;\;\;\;\;\;\;\;\;\frac{{{a_1}x_1^*\left( {1 + {\rho _1}y_2^*} \right)}}{{1 + x_1^*}} - \frac{{{\rho _1}{a_2}x_2^*y_2^*}}{{1 + x_2^*}} - {d_1}\;\;\;\;\;\; - \frac{{{\rho _1}{a_2}y_1^*y_2^*}}{{{{\left( {1 + x_2^*} \right)}^2}}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rho _1}y_1^*\left( {\frac{{{a_1}x_1^*}}{{1 + x_1^*}} - \frac{{{a_2}x_2^*}}{{1 + x_2^*}}} \right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;r\left( {1 - \frac{{2x_2^*}}{{{K_2}}}} \right) - \frac{{{a_2}y_2^*}}{{{{\left( {1 + x_2^*} \right)}^2}}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \frac{{{a_2}x_2^*}}{{1 + x_2^*}}\\ \;\;\;\;\;\;\; - \frac{{{\rho _2}{a_1}y_1^*y_2^*}}{{{{\left( {1 + x_1^*} \right)}^2}}}\;\;\;\;\;\;\;\;\;\;{\rho _2}y_2^*\left( {\frac{{{a_2}x_2^*}}{{1 + x_2^*}} - \frac{{{a_1}x_1^*}}{{1 + x_1^*}}} \right)\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{{{a_2}y_2^*\left( {1 + {\rho _2}y_1^*} \right)}}{{{{\left( {1 + x_2^*} \right)}^2}}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{{{a_2}x_2^*\left( {1 + {\rho _2}y_1^*} \right)}}{{1 + x_2^*}} - \frac{{{\rho _2}{a_1}x_1^*y_1^*}}{{1 + x_1^*}} - {d_2} \end{array}$$ \right| \end{array} $

(14)

After substituting the boundary equilibria $E_{0000}$, $E_{K_1000}$, $E_{00K_20}$, $E_{\mu_1\nu_100}$ and $E_{00\mu_2\nu_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_{K_10K_20}$ are as follows:

$\lambda_1 = -1~(<0),~~\lambda_2 = \frac{a_1K_1}{1+K_1}-d_1<0\Leftrightarrow\mu_1>K_1,$

$~~\lambda_3 = -r~(<0),~~\lambda_4 = \frac{a_2K_2}{1+K_2}-d_2<0\Leftrightarrow\mu_2>K_2.$

Therefore, $E_{K_10K_20}$ is locally asymptotically stable if $\mu_i>K_i,i=1,2$ while it is a saddle if either $\left(\mu_1-K_1\right)\left(\mu_2-K_2\right)<0$ or $\mu_i<K_i, i=1,2$ holds.

Now we focus on the local stability of $E_{\mu_1\nu_1K_20}$ and $E_{K_10\mu_2\nu_2}$ when they exist. After substituting the boundary equilibrium $E_{\mu_1\nu_1K_20}$ to (14), we can obtain the eigenvalues of the Jacobian matrix evaluated at this boundary equilibrium as follows:

$\lambda_1 = -r,\,\lambda_2=\frac{K_2(a_2-d_2)-d_2}{1+K_2}+\rho_2\frac{\nu_1\left[K_2(a_2-d_1)-d_1\right]}{(1+K_2)}$

and

$\lambda_3+\lambda_4=\frac{K_1(a_1-d_1)-(a_1+d_1)}{a_1K_1(a_1-d_1)},\,\lambda_3\lambda_4=\frac{d_1(K_1(a_1-d_1)-d_1)}{a_1K_1}.$

Notice that the eigenvalues of $\lambda_3$ and $\lambda_4$ being negative is equivalent to the case that the unique interior equilibrium $(\mu_1,\nu_1)$ being locally asymptotically stable for the single patch model (5) when $i=1$. Thus, we can conclude that $\frac{K_1-1}{2}<\mu_1<K_1$ are sufficient conditions for $\lambda_3$ and $\lambda_4$ being negative. Now we explore sufficient conditions for $\lambda_2$ being negative. First, we have $\mu_1<K_1$ due to the existence of $E_{\mu_1\nu_1K_20}$. We have the following three cases:

1.If $\mu_2>K_2 \Leftrightarrow K_2(a_2-d_2)-d_2<0$, then the first term of $\lambda_2$ is negative. This also implies that Model (4) has no boundary equilibria of $E_{00\mu_2\nu_2}$ and $E_{K_10\mu_2\nu_2}$. Since $\mu_1<K_1$, therefore, we have $\lambda_2<0$ for all $\rho_2>0$ if $a_2\leq d_1$ or $K_2<\frac{d_1}{a_2-d_1}$ since

$K_2(a_2-d_1)-d_1<0 \Leftrightarrow either a_2\leq d_1 \mbox{ or }K_2<\frac{d_1}{a_2-d_1}.$

Therefore, we can conclude that $\lambda_2$ is negative if either $a_2\leq d_1, K_2<\mu_2$ or $K_2<\min\Big\{\mu_2,\frac{d_1}{a_2-d_1}\Big\}.$ Assume that $K_2(a_2-d_1)-d_1>0$ (i.e., $K_2>\frac{d_1}{a_2-d_1}>0$), then $\lambda_2<0$ if $\rho_2$ is small enough, i.e., satisfies the condition of $\rho_2<\frac{d_2-K_2(a_2-d_2)}{\nu_1\left[K_2(a_2-d_1)-d_1\right]}.$ In this case, we can conclude that $\lambda_2$ is negative if

$0<\frac{d_1}{a_2-d_1}<K_2<\mu_2 \;and\;\rho_2<\frac{d_2-K_2(a_2-d_2)}{\nu_1\left[K_2(a_2-d_1)-d_1\right]}.$

2.If $\mu_2<K_2 \Leftrightarrow K_2(a_2-d_2)-d_2>0$, then the first term of $\lambda_2$ is positive. This also implies that Model (4) has two boundary equilibria of $E_{00\mu_2\nu_2}$ and $E_{K_10\mu_2\nu_2}$. In this case, sufficient conditions for $\lambda_2$ being negative are $K_2<\frac{d_1}{a_2-d_1}$ and $\rho_2$ large enough. More specifically, $\rho_2$ has to satisfy the following inequality:

$\rho_2>\frac{K_2(a_2-d_2)-d_2}{\nu_1\left[d_1-K_2(a_2-d_1)\right]}.$

Therefore, $\lambda_2$ is negative if $\mu_2<K_2<\frac{d_1}{a_2-d_1} \;{\rm{and}}\;\rho_2>\frac{K_2(a_2-d_2)-d_2}{\nu_1\left[d_1-K_2(a_2-d_1)\right]}.$

Summarizing the discussions above, we can conclude that the boundary equilibrium $E_{\mu_1\nu_1K_20}$ is locally asymptotically stable if $\frac{K_1-1}{2}<\mu_1<K_1$ and one of the following conditions holds:

1.$a_2\leq d_1, K_2<\mu_2$.

2.$K_2<\min\Big\{\mu_2,\frac{d_1}{a_2-d_1}\Big\}.$

3.$0<\frac{d_1}{a_2-d_1}<K_2<\mu_2 \;and\;\rho_2<\frac{d_2-K_2(a_2-d_2)}{\nu_1\left[K_2(a_2-d_1)-d_1\right]}.$

4.$\mu_2<K_2<\frac{d_1}{a_2-d_1} \;and\;\rho_2>\frac{K_2(a_2-d_2)-d_2}{\nu_1\left[d_1-K_2(a_2-d_1)\right]}.$

And $E_{\mu_1\nu_1K_20}$ is a saddle if $\mu_1<\frac{K_1-1}{2}$ or one of the following conditions holds:

1.$K_2>\max\Big\{\mu_2,\frac{d_1}{a_2-d_1}\Big\}.$

2.$0<\frac{d_1}{a_2-d_1}<K_2<\mu_2 \;and\;\rho_2>\frac{d_2-K_2(a_2-d_2)}{\nu_1\left[K_2(a_2-d_1)-d_1\right]}.$

3.$\mu_2<K_2<\frac{d_1}{a_2-d_1} \;and\;\rho_2<\frac{K_2(a_2-d_2)-d_2}{\nu_1\left[d_1-K_2(a_2-d_1)\right]}.$

Similarly, we can obtain sufficient conditions for the local stability of the boundary equilibrium $E_{K_10\mu_2\nu_2}$ as the statement.

If $\mu_i<K_i$, then Model (4) has the boundary equilibria $E_{\mu_1\nu_1K_20}$ and $E_{K_10\mu_2\nu_2}$ according to Theorem 3.2 and the discussions above. If both $E_{\mu_1\nu_1K_20}$ and $E_{K_10\mu_2\nu_2}$ are locally stable, then the following inequalities are satisfied:

$\mu_1<K_1<\frac{d_2}{a_1-d_2} \Rightarrow \frac{d_1}{a_1-d_1}<\frac{d_2}{a_1-d_2}\Rightarrow d_1<d_2$

and

$\mu_2<K_2<\frac{d_1}{a_2-d_1} \Rightarrow \frac{d_2}{a_2-d_2}<\frac{d_1}{a_2-d_1}\Rightarrow d_2<d_1$

which are contradiction. Therefore, $E_{\mu_1\nu_1K_20}$ and $E_{K_10\mu_2\nu_2}$ can not be local stable at the same time.

Now if $r_i=1, a_i=a,d_i=d, K_i=d$ for both $i=1,2$, then $E_{i2}^b, i=1,2$ exist if $Ka-Kd-d>0$. This implies that one of the eigenvalues of the Jacobian matrix of Model (4) evaluated at $E_{i2}^b$ is positive, i.e.,

$\frac{K(a-d)^(Ka-Kd-d)+\rho_j (Ka-Kd-d)^2}{K(K+1)(a-d)^2}>0$

which indicates that $E_{i2}^b$ can not be stable for both $i=1$ and $i=2$.

Proof of Theorem 3.3.

Proof.Let $p_i(x)=\frac{a_i x}{1+x}$ and $q_i(x)=\frac{r_i(K_i-x)(1+x)}{a_iK_i}$, then we have

$r_ix_i\left(1-\frac{x_i}{K_i}\right)-\frac{a_i x_iy_i}{(1+x_i)}=\frac{a_i x_i}{1+x_i}\left[\frac{r_i(K_i-x_i)(1+x_i)}{a_iK_i}-y_i\right]$

$=p_i(x_i)\left[q_i(x_i)-y_i\right].$

We construct the following Lyapunov functions

$ $$\begin{aligned} V_1(x_1,y_1) =\rho_2 \int^{x_1}_{K_1} \frac{p_1(\xi)-p_1(K_1)}{p_1(\xi)}d\xi + \rho_2y_1 \end{aligned}$$ $

(15)

and

$ $$\begin{aligned} V_2(x_2,y_2) =\rho_1 \int^{x_2}_{K_2} \frac{p_2(\xi)-p_2(K_2)}{p_2(\xi)}d\xi + \rho_1y_2 \end{aligned}$$ $

(16)

Now taking derivatives of the functions (15) and (16) with respect to time $t$, we get

$ $$\begin{array}{l} &\frac{d}{dt}V_1(x_1(t),y_1(t))\\ =& \rho_2\frac{p_1(x_1)-p_1(K_1)}{p_i(x_1)}\frac{dx_1}{dt} + \rho_2\frac{dy_1}{dt} \\ =& \rho_2\left[ p_1(x_1)-p_1(K_1)\right]\left[ q_1(x_1)-y_1\right] +\rho_2 y_1\left[ p_1(x_1)-d_1\right] + \rho_1\rho_2 y_1y_2\left[ p_1(x_1)-p_2(x_2)\right]\\ =&\rho_2\left[ p_1(x_1)-p_1(K_1)\right] q_1(x_1)+ \rho_2 y_1\left[ p_1(K_1)-d_1\right]+\rho_1\rho_2 y_1y_2\left[ p_1(x_1)-p_2(x_2)\right] \end{array}$$ $

(17)

and

$ $$\begin{array}{l} \;\;\;\;\;\;\;\;\frac{d}{{dt}}{V_2}({x_2}(t),{y_2}(t))\\ = {\rho _1}\frac{{{p_2}({x_2}) - {p_2}({K_2})}}{{{p_2}({x_2})}}\frac{{d{x_2}}}{{dt}} + {\rho _1}\frac{{d{y_2}}}{{dt}}\\ = {\rho _1}\left[ {{p_2}({x_2}) - {p_2}({K_2})} \right]\left[ {{q_2}({x_2}) - {y_2}} \right] + {\rho _1}{y_2}\left[ {{p_2}({x_2}) - {d_2}} \right] + {\rho _1}{\rho _2}{y_1}{y_2}\left[ {{p_2}({x_2}) - {p_1}({x_1})} \right].\\ = - {\rho _1}\left[ {{p_2}({x_2}) - {p_2}({K_2})} \right]{q_2}({x_2}) + {\rho _1}{y_2}\left[ {{p_2}({K_2}) - {d_2}} \right] + {\rho _1}{\rho _2}{y_1}{y_2}\left[ {{p_2}({x_2}) - {p_1}({x_1})} \right] \end{array}$$ $

(18)

Let $V=V_1+V_2$. Now adding (17) and (18), we get

$\nonumber $$\begin{aligned} \frac{d}{dt}V=&\frac{d}{dt}V_1(x_1(t),y_1(t)) + \frac{d}{dt}V_2(x_2(t),y_2(t))\\ =& \rho_2\left[ p_1(x_1)-p_1(K_1)\right] q_1(x_1)+ \rho_2 y_1\left[ p_1(K_1)-d_1\right] \\ &+ \rho_1\left[ p_2(x_2)-p_2(K_2)\right] q_2(x_2)+ \rho_1 y_2\left[ p_2(K_2)-d_2\right]. \end{aligned}$$ $

Since $\mu_i> K_i\Leftrightarrow \frac{d_i}{a_i-d_i}>K_i\Leftrightarrow \frac{a_iK_i}{1+K_i}=p_i(K_i)<d_i$, thus, we have $p_i(K_i)-d_i<0$. Notice that $p_i(x_i)$ is an increasing function in $x_i$, therefore, $\left[ p_i(x_i)-p_i(K_i)\right]$ is 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 $x_i>K_i$. This implies that $\left[ p_i(x_i)-p_i(K_i)\right] q_i(x_i) \leq 0$ for all $x_i\geq0$. Therefore, we have $\frac{d}{dt}V<0$ in $\mathbb R^4_+$. This implies that both $V_1$ and $V_2$ are Lyapunov functions, and the boundary equilibrium $E_{K_10K_20} = \left( K_1,0,K_2,0\right)$ is globally stable when $\mu_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 $\mathbb 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

$\frac{dx_1}{x_1dt}\Big\vert_{x1=0,y_1=0} = \left(1 - \frac{x_1}{K_1}\right) - \frac{a_1 y_1}{1 + x_1} \Big\vert_{x1=0,y_1=0} =1>0.$

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_{K_1000}\,\;and\;\,E_{00K_20}$. 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_{K_1000},\,,\,E_{K_10K_20},\,E_{K_10\mu_2\nu_2}$ if $\frac{K_2-1}{2}<\mu_2<K_2$ while they are $E_{K_1000},\,E_{K_10K_20}$ if $\mu_2>K_2$. Now we consider the following two cases:

1.If $\mu_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:

$ $$\begin{array}{lcl} \frac{dy_1}{y_1dt} \Big\vert_{E_{K_1000}}& =& \left[ \frac{a_1 x_1 }{1 + x_1} -d_1+ \rho_1\left( \frac{a_1 x_1 y_2}{1 + x_1}- \frac{a_2 x_2 y_2}{1 + x_2}\right)\right] \Big\vert_{E_{K_1000}}=\frac{a_1 K_1 }{1 + K_1} -d_1\\ \frac{dy_1}{y_1dt} \Big\vert_{E_{K_10K_20}}& =& \left[ \frac{a_1 x_1 }{1 + x_1} -d_1+ \rho_1\left( \frac{a_1 x_1 y_2}{1 + x_1}- \frac{a_2 x_2 y_2}{1 + x_2}\right)\right] \Big\vert_{E_{K_10K_20}}=\frac{a_1 K_1 }{1 + K_1} -d_1 \end{array}$$ .$

Since $\frac{a_1 K_1 }{1 + K_1} -d_1>0 \Leftrightarrow \mu_1<K_1$, therefore, we can conclude that predator $y_1$ is persistent if $\mu_1<K_1$ and $\mu_2>K_2$.

2.If $\frac{K_2-1}{2}<\mu_2<K_2$, according to Theorem 2.5 in [33] and discussions above, we can conclude that predator $y_1$ is persistent if $\mu_1<K_1$ and the following equation is strictly positive:

$ $$\begin{array}{l} &\frac{dy_1}{y_1dt} \Big\vert_{E_{K_10\mu_2\nu_2}}\\ =& \left[ \frac{a_1 x_1 }{1 + x_1} -d_1+ \rho_1\left( \frac{a_1 x_1 y_2}{1 + x_1}- \frac{a_2 x_2 y_2}{1 + x_2}\right)\right] \Big\vert_{E_{K_10\mu_2\nu_2}}\\ =&\frac{a_1 K_1 }{1 + K_1} -d_1+ \rho_1\left( \frac{a_1 K_1 \nu_2}{1 + K_1}- \frac{a_2 \mu_2 \nu_2}{1 + \nu_2}\right)\\ =&\frac{a_1 K_1 }{1 + K_1} -d_1+ \rho_1 \nu_2\left( \frac{a_1 K_1}{1 + K_1}- \frac{a_2 \mu_2}{1 + \mu_2}\right)=\frac{a_1 K_1 }{1 + K_1} -d_1+ \rho_1 \nu_2\left( \frac{a_1 K_1}{1 + K_1}-d_2\right)>0 \end{array}$$ .$

According to the proof of Theorem 3.2, we can see that sufficient condition that $\frac{dy_1}{y_1dt} \Big\vert_{E_{K_10\mu_2\nu_2}}>0$ holds is the same as sufficient condition for the boundary equilibrium $E_{K_10\mu_2\nu_2}$ being unstable when $\mu_1<K_1$. Therefore, we can conclude that predator $y_1$ is persistent if one of the following inequalities hold

(a)$\mu_j<K_j, \mu_i>K_i$. Or

(b)$\frac{K_i-1}{2}<\mu_i<K_i \;and\; K_j>\max\Big\{\mu_j,\frac{d_i}{a_j-d_i}\Big\}$. Or

(c)$\frac{K_i-1}{2}<\mu_i<K_i,\,\,\mu_j<K_j<\frac{d_i}{a_j-d_i} \;and\;\rho_j<\frac{K_j(a_j-d_j)-d_j}{\nu_i\left[d_i-K_j(a_j-d_i)\right]}$.

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 $\frac{K_j-1}{2}<\mu_j<K_j, \mu_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 $\frac{K_j-1}{2}<\mu_j<K_j$, then Theorem 3.2 indicates that the omega limit set of Model (4) when $y_i=0$ is $E_{\mu_1\nu_1K_20}$ when $i=2,j=1$ while its omega limit set is $E_{K_10\mu_2\nu_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:

$ $$\begin{array}{l} &\frac{dy_1}{y_1dt} \Big\vert_{E_{K_10\mu_2\nu_2}}\\ =& \left[ \frac{a_1 x_1 }{1 + x_1} -d_1+ \rho_1\left( \frac{a_1 x_1 y_2}{1 + x_1}- \frac{a_2 x_2 y_2}{1 + x_2}\right)\right] \Big\vert_{E_{K_10\mu_2\nu_2}}\\ =&\frac{a_1 K_1 }{1 + K_1} -d_1+ \rho_1\left( \frac{a_1 K_1 \nu_2}{1 + K_1}- \frac{a_2 \mu_2 \nu_2}{1 + \nu_2}\right)\\ =&\frac{a_1 K_1 }{1 + K_1} -d_1+ \rho_1 \nu_2\left( \frac{a_1 K_1}{1 + K_1}- \frac{a_2 \mu_2}{1 + \mu_2}\right)=\frac{a_1 K_1 }{1 + K_1} -d_1+ \rho_1 \nu_2\left( \frac{a_1 K_1}{1 + K_1}-d_2\right)>0 \end{array}$$ .$

Since $\mu_i>K_i\Leftrightarrow \frac{a_1 K_1 }{1 + K_1} -d_1<0$, therefore, sufficient conditions for $\frac{dy_1}{y_1dt} \Big\vert_{E_{K_10\mu_2\nu_2}}>0$ is

$\frac{a_1 K_1}{1 + K_1}-d_2>0 \Leftrightarrow K_1>\frac{d_2}{a_1-d_2} \;and\; \rho_1>\frac{d_1-K_1(a_1-d_1)}{v_2\left[K_1(a_1-d_2)-d_2\right]}.$

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$,

$\frac{K_j-1}{2}<\mu_j<K_j, \mu_i>K_i,\,0<\frac{d_j}{a_i-d_j}<K_i<\mu_i $

$\;and\; \rho_i>\frac{d_i-K_i(a_i-d_i)}{\nu_j\left[K_i(a_i-d_j)-d_j\right]}.$

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

$\frac{K_i-1}{2}<\mu_i<K_i,\, \mu_j>\frac{K_j-1}{2} \;and\; K_j>\max\Big\{\mu_j,\frac{d_i}{a_j-d_i}\Big\}.$

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

$\frac{K_i-1}{2}<\mu_i<K_i,\,\frac{K_j-1}{2}<\mu_j<K_j<\frac{d_i}{a_j-d_i} \;and\;\rho_j<\frac{K_j(a_j-d_j)-d_j}{\nu_i\left[d_i-K_j(a_j-d_i)\right]}.$

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$,

$\frac{K_i-1}{2}<\mu_i<K_i,\,K_i>\max\Big\{\mu_i,\frac{d_j}{a_i-d_j}\Big\},\,$

$ \frac{K_j-1}{2}<\mu_j<K_j<\frac{d_i}{a_j-d_i} \;and\;\rho_j<\frac{K_j(a_j-d_j)-d_j}{\nu_i\left[d_i-K_j(a_j-d_i)\right]}. $

Based on the discussion above, we can conclude that the statement of Theorem 3.7 holds.

Proof of Theorem 3.8.

Proof.If $\mu_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 $\mu_i>K_i$ for both $i=1,2$.

The interior equilibrium $(x^*_1,y^*_1,x^*_2,y_2^*)$ is determined by the positive intersections of the nullclines (11) $x_1=F(x_2)=\frac{f_t(x_2)}{f_b(x_2)}$ and $x_2=G(x_1)=\frac{g_t(x_1)}{g_b(x_1)}$ where

$ $$\begin{align*} f_t(x_2)&=a_2\left[r_2\rho_1x_2\left(K_2-x_2\right)+K_2d_1\right],\,f_b(x_2)=r_2\rho_1x_2\left(K_2a_1-K_2a_2-a_1\right)\\ &-r_2\rho_1x_2^2(a_1-a_2)+K_2(a_1r_2\rho_1+a_1a_2-a_2d_1) \end{align*}$$ $

and

$ $$\begin{align*} g_t(x_1)&=a_1\left[r_1\rho_2x_1\left(K_1-x_1\right)+K_1d_2\right],\, g_b(x_1)=r_1\rho_2x_1\left(K_1a_2-K_1a_1-a_2\right)\\ &-r_1\rho_2x_1^2(a_2-a_1)+K_1(a_2r_1\rho_2+a_1a_2-a_1d_2). \end{align*}$$ $

Notice that the nullclines $x_1=F(x_2)=\frac{f_t(x_2)}{f_b(x_2)}$ and $x_2=G(x_1)=\frac{g_t(x_1)}{g_b(x_1)}$ has the following properties:

1.$F(0)=\frac{f_t(0)}{f_b(0)}=\frac{a_2K_2d_1}{K_2(a_1r_2\rho_1+a_1a_2-a_2d_1)}=\frac{a_2d_1}{a_1r_2\rho_1+a_1a_2-a_2d_1}$ and $F(K_2)=\frac{f_t(K_2)}{f_b(K_2)}=\frac{d_1}{a_1-d_1}=\mu_1$.

2.$f_t(x_2)=a_2\left[r_2\rho_1x_2\left(K_2-x_2\right)+K_2d_1\right]\geq a_2K_2 d_1>0forx_2\in [0, K_2]$ and

$f_b(x_2)\big\vert_{a_1=a_2=a}=a\left[r_2\rho_1(K_2-x_2)+K_2(a-d_1)\right].$

3.$G(0)=\frac{g_t(0)}{g_b(0)}=\frac{a_1K_1d_2}{K_1(a_2r_1\rho_2+a_1a_2-a_1d_2)}=\frac{a_1d_2}{a_2r_1\rho_2+a_1a_2-a_1d_2}$ and $G(K_1)=\frac{g_t(K_1)}{g_b(K_1)}=\frac{d_2}{a_2-d_2}=\mu_2$.

4.$g_t(x_1)=a_1\left[r_1\rho_2x_1\left(K_1-x_1\right)+K_1d_2\right]\geq a_1 K_1d_1>0\mbox{ for } x_1\in [0, K_1]$ and

$g_b(x_1)\big\vert_{a_1=a_2=a}=a\left[r_1\rho_2(K_1-x_1)+K_1(a-d_2)\right].$

According to Theorem 3.1, we know that population of prey $x_i$ for $i=1,2$ has the following properties:

$\limsup\limits_{t\rightarrow\infty} x_i(t)\leq K_i.$

Thus, we can restrict the function $F(x_2)$ on the domain of $[0, K_2]$ and $G(x_1)$ on the domain of $[0, K_1]$. Since $f_t(x_2)\geq a_2K_2 d_1>0forx_2\in [0, K_2]$ and $g_t(x_1)\geq a_1 K_1d_1>0\mbox{ for } x_1\in [0, K_1]$, thus, Model (4) has no interior equilibrium if $f_b(x_2)<0$for$x_2\in [0, K_2]$ or $g_b(x_1)<0$for$x_1\in [0, K_1]$.

Now we assume that $a_1>a_2$, then we have $f_b(x_2)<0$for$x_2\in [0, K_2]$ if

$r_2^2\rho_1^2\left(K_2a_1-K_2a_2-a_1\right)^2<4K_2r_2\rho_1(a_2-a_1)(a_1r_2\rho_1+a_1a_2-a_2d_1)$

$\Leftrightarrow \rho_1<\frac{4K_2a_2(a_1-a_2)(d_1-a_1)}{r_2\left(K_2a_1-K_2a_2+a_1\right)^2}$

while $f_b(x_2)>0$for$x_2\in [0, K_2]$ if $a_1>d_1$ since

$f_b(0)=K_2(a_1r_2\rho_1+a_2(a_1-d_1))>0 \;and\; f_b(K_2)=a_2K_2(a_1-d_1)>0.$

And $g_b(x_1)>0$ for $x_1\in [0, K_1]$ if

$r_1^2\rho_2^2\left(K_1a_2-K_1a_1-a_2\right)^2<4K_1r_1\rho_2(a_1-a_2)(a_2r_1\rho_2+a_1a_2-a_1d_2)$

$ \Leftrightarrow \rho_2<\frac{4K_1a_1(a_1-a_2)(a_2-d_2)}{r_1\left(K_1a_2-K_1a_1+a_2\right)^2}. $

Similar cases can be made for $a_1<a_2$, therefore we can conclude that Model (4) has no interior equilibrium if either

$ a_1>a_2, \,\rho_1<\frac{4K_2a_2(a_1-a_2)(d_1-a_1)}{r_2\left(K_2a_1-K_2a_2+a_1\right)^2}$

or

$a_2>a_1,\,\rho_2<\frac{4K_1a_1(a_2-a_1)(d_2-a_2)}{r_1\left(K_1a_2-K_1a_1+a_2\right)^2}$

hold.

Now we focus on sufficient conditions lead to both $f_b(x_2)>0$ for $x_2\in [0, K_2]$ and $g_b(x_1)>0$for$x_1\in [0, K_1]$. Assume that $a_1>a_2$ and $a_i>d_i$ for both $i=1,2$. Notice that $f_b(0)=K_2(a_1r_2\rho_1+a_2(a_1-d_1))> f_b(K_2)=a_2K_2(a_1-d_1)>0$ with

$f_b(x_2)=r_2\rho_1x_2\left(K_2a_1-K_2a_2-a_1\right)-r_2\rho_1x_2^2(a_1-a_2)+K_2(a_1r_2\rho_1+a_1a_2-a_2d_1).$

Therefore, we have $f_b(x_2)>f_b(K_2)$ for all $x_2\in [0, K_2]$. Since $g_b(x_1)$ is a degree 2 polynomial with the positive coefficient in the degree 2 and $g_b(0)=K_1(a_2r_1\rho_2+a_1(a_2-d_2))> g_b(K_1)=a_1K_1(a_2-d_2)>0$. Therefore, we have $g_b(x_1)>0$ for all $x_1\in [0, K_1]$ if the following conditions hold

$\rho_2<\frac{4K_1a_1(a_2-a_1)(d_2-a_2)}{r_1\left(K_1a_2-K_1a_1+a_2\right)^2}.$

The discussion so far also indicates that we have both $f_b(x_2)>0$for$x_2\in [0, K_2]$ and $g_b(x_1)>0$for$x_1\in [0, K_1]$ if the following inequalities hold for $i=1,j=2$ or $i=2,j=1$

$a_i>\max\{a_j,d_i\},\, \rho_j<\frac{4K_ia_i(a_j-a_i)(d_j-a_j)}{r_i\left(K_ia_j-K_ia_i+a_j\right)^2} .$

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

$x^c_1=\frac{K_1\left(r_1\rho_2+a_1-d_2-\sqrt{(a_1-d_2)(r_1\rho_2+a_1-d_2)}\right)}{r_1\rho_2} \in (0,K_1)$

and

$x^c_2=\frac{K_2\left(r_2\rho_1+a_2-d_1-\sqrt{(a_2-d_1)(r_2\rho_1+a_2-d_1)}\right)}{r_2\rho_1} \in (0,K_2).$

If $F(x_2^c)<K_1 \;and\; G(x_1^c)<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\times G$ maps $[0,K_2]\times[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\in [0,K_2], \, x_1\in [0,K_1]$. Now we focus on sufficient condition that leads to $F(x_2^c)<K_1 \;and\; G(x_1^c)<K_2$ when $a_1>a_2$. Since

$\max\limits_{0\leq x_2\leq K_2}\{F(x_2)\}=F(x_2^c)\leq\frac{\max\limits_{0\leq x_2\leq K_2}\{f_t(x_2)\}}{\min\limits_{0\leq x_2\leq K_2}\{f_b(x_2)\}}=\frac{f_b(K_2/2)}{f_b(K_2)}=\frac{K_2r_2\rho_1+4d_1}{4(a_1-d_1)}$

and

$\max\limits_{0\leq x_1\leq K_1}\{G(x_1)\}=G(x_1^c)\leq\frac{\max\limits_{0\leq x_1\leq K_1}\{g_t(x_1)\}}{\min\limits_{0\leq x_1\leq K_1}\{g_b(x_1)\}}=\frac{g_t(K_1/2)}{g_b\left(\frac{K_2a_1-K_2a_2-a_1}{2(a_1-a_2)}\right)}$

$=\frac{a_1K_1(r_1K_1\rho_2/4+d_2)}{K_1a_1(a_2-d_2)-\frac{r_1\rho_2(K_1a_1-K_1a_2-a_2)^2}{4(a_1-a_2)}},$

therefore, we have $F(x_2^c)<K_1 \;and\; G(x_1^c)<K_2$ when $a_1>a_2$ if the following inequalities hold

$\frac{K_2r_2\rho_1+4d_1}{4(a_1-d_1)}\leq K_1 \Leftrightarrow \rho_1\leq \frac{4(K_1a_1-K_1d_1-d_1)}{K_2r_2}.$

and

$ $$\begin{array}{lcl} K_2&\geq&\frac{a_1K_1(r_1K_1\rho_2/4+d_2)}{K_1a_1(a_2-d_2)-\frac{r_1\rho_2(K_1a_1-K_1a_2-a_2)^2}{4(a_1-a_2)}}\\ &\Leftrightarrow&\\ \rho_2&\leq&\frac{4K_1a_1(K_2a_2-K_2d_2-d_2)}{a_1r_1K_1^2+r_2K_2(K_1a_1-K_1a_2-a_2)^2}. \end{array}$$ $

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_i>\max\{a_j,d_1,d_2\}, a_j>\max\{d_1,d_2\}, \,\rho_i\leq \frac{4(K_ia_i-K_id_i-d_i)}{K_jr_j}$

and

$\rho_j<\min\Big\{\frac{4K_ia_i(a_j-a_i)(d_j-a_j)}{r_i\left(K_ia_j-K_ia_i+a_j\right)^2},\frac{4K_ia_i(K_ja_j-K_jd_j-d_j)}{a_ir_iK_i^2+r_jK_j(K_ia_i-K_ia_j-a_j)^2}\Big\}.$

In addition, since both $F(x_2)$ and $G(x_1)$ are unimode maps in their domain with unique local maximum, thus, we have

$x_1=F(x_2)\geq F(0)=\frac{a_2d_1}{a_1r_2\rho_1+a_1a_2-a_2d_1}$

$\;and\;x_2=G(x_1)\geq G(0)=\frac{a_1d_2}{a_2r_1\rho_2+a_1a_2-a_1d_2}.$

Therefore, we have $\frac{a_id_j}{a_jr_i\rho_j+a_ia_j-a_id_j}<x_j^*<K_j$ for both $i=1,j=2$ and $i=2,j=1$.

Now assume that $a_1=a_2=a$, then both $f_b$ and $g_b$ are reduced to linear decreasing functions, i.e.,

$f_b(x_2)\big\vert_{a_1=a_2=a}=a\left[r_2\rho_1(K_2-x_2)+K_2(a-d_1)\right] \;and\;g_b(x_1)\big\vert_{a_1=a_2=a}$

$ =a\left[r_1\rho_2(K_1-x_1)+K_1(a-d_2)\right] $

which indicates that $f_b(x_2)<0$ for $x_2\in [0, K_2]$ if $d_1>a+r_2\rho_1$ and $g_b(x_1)<0$for$x_1\in [0, K_1]$ if $d_2>a+r_1\rho_2$. Therefore, if $a_1=a_2=a$ and either $d_1>a+r_2\rho_1$ or $d_2>a+r_1\rho_2$ holds, then Model (4) has no interior equilibrium. On the other hand, both $f_b(x_2)>0$for$x_2\in [0, K_2]$ and $g_b(x_1)>0$for$x_1\in [0, K_1]$ if $a>\max\{d_1,d_2\}$. Then apply the discussions for the case $a_1\neq a_2$, then we can conclude that Model (4) has at least one interior equilibrium if

$a>\max\{d_1,d_2\},\,F(x^c_2)<K_1,\,\;and\; G(x^c_1)<K_2.$

Applying the similar arguments for the case $a_i>a_j$, we have

$\max\limits_{0\leq x_2\leq K_2}\{F(x_2)\}=F(x_2^c)\leq\frac{\max_{0\leq x_2\leq K_2}\{f_t(x_2)\}}{\min_{0\leq x_2\leq K_2}\{f_b(x_2)\}}=\frac{f_b(K_2/2)}{f_b(K_2)}=\frac{K_2r_2\rho_1+4d_1}{4(a_1-d_1)}$

and

$ \max\limits_{0\leq x_1\leq K_1}\{G(x_1)\}=G(x_1^c)\leq\frac{\max\limits_{0\leq x_1\leq K_1}\{g_t(x_1)\}}{\min\limits_{0\leq x_1\leq K_1}\{g_b(x_1)\}}=\frac{g_t(K_1/2)}{g_b\left(K_1\right)}=\frac{K_1r_1\rho_2+4d_2}{4(a_2-d_2)}. $

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\}, \rho_i<\frac{4(K_ia-K_id_i-d_i)}{K_jr_j}$

for both $i=1,j=2$ and $i=2,j=1$. In addition, $\frac{d_j}{r_i\rho_j+a-d_j}<x_j^*<K_j$ hold for both $i=1,j=2$ and $i=2,j=1$.

Proof of Theorem 3.9.

Proof.Suppose that $\min\{a_1,a_2\}>d$, then we have $\mu_i=\frac{d}{a_i-d}, \nu_i=\frac{r_i(K_i-\mu_i)(1+\mu_i)}{a_iK_i}$ where $(\mu_i,\nu_i)$ is the unique interior equilibrium of the single patch model (5) in the absence of the dispersal in predator for both $i=1,2$. Now recall from the null clines (9), we have

$ $$\begin{array}{lcl} x_i(x_j)&=&\frac{\rho_i q_j(x_j)p_j(x_j)+d_i}{a_i(1+\rho_i q_j(x_j))-(\rho_i q_j(x_j)p_j(x_j)+d_i)}=\frac{\rho_i q_j(x_j)p_j(x_j)+d_i}{\rho_i q_j(x_j)\left[a_i-p_j(x_j)\right]+a_i-d_i} \end{array}$$ $

which indicates that

$x_i(\mu_j)= \frac{\rho_i q_j(\mu_j)p_j(\mu_j)+d}{\rho_i q_j(\mu_j)\left[a_i-p_j(\mu_j)\right]+a_i-d}=\frac{d}{a_i-d}=\mu_i.$

This implies that $x_i=\mu_i$ for $i=1,2$ is the positive solutions of the null clines (9). Therefore, we can solve that $(\mu_1,\nu_1,\mu_2,\nu_2)$ is an interior solution of Model (4) if $\min\{a_1,a_2\}>d=d_1=d_2$. By substituting the equilibrium $(\mu_1,\nu_1,\mu_2,\nu_2)$ into the Jacobian matrix (14), we obtain its characteristic equation as follows:

$H(\lambda)=\lambda^4+(\alpha_1+\alpha_2)\lambda^3+[\alpha_1\alpha_2+d(\beta_1+\beta_2)]\lambda^2+d(\alpha_1\beta_2+\alpha_2\beta_1)\lambda+d^2(\beta_1\beta_2-\gamma_1\gamma_2=0$

where

$ $$\begin{array}{l} {\alpha _i}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; = - \frac{{{r_i}{\mu _i}({K_i}{a_i} - {K_i}d - {a_i} - d)}}{{{K_i}{a_i}}}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \Rightarrow [{\alpha _i} > 0 \Leftrightarrow {K_i}{a_1} - {K_i}d - {a_i} - d < 0 \Leftrightarrow \frac{{{K_i} - 1}}{2} < {\mu _i} < {K_i}]\\ {\beta _i}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; = \frac{{{\nu _i}({\nu _j}{\rho _i} + 1){{({a_i} - d)}^2}}}{{{a_i}}} > 0\\ {\gamma _i}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; = \frac{{{\rho _i}{\nu _i}{\nu _j}{{({a_j} - d)}^2}}}{{{a_j}}} > 0\\ {\beta _1}{\beta _2} - {\gamma _1}{\gamma _2} = \frac{{{\nu _1}{\nu _2}{{({a_1} - d)}^{({a_2} - d)}}^2({\nu _1}{\rho _2} + {\nu _2}{\rho _1} + 1)}}{{{a_1}{a_2}}} > 0. \end{array}$$ $

Then the real parts of the solutions of $H(\lambda)$ are all negative if $\alpha_1+\alpha_2>0$ while the solutions of $H(\lambda)$ has positive if $\alpha_1+\alpha_2<0$. Notice that the single patch $i$ (5) has global stability at $(\mu_i,\nu_i)$ if $\alpha_i>0\Leftrightarrow \frac{K_i-1}{2}<\mu_i<K_i.$ Therefore, the interior equilibrium is locally asymptotically stable of $\frac{K_i-1}{2}<\mu_i<K_i$ for both $i=1,2$ while it is unstable if

$\frac{\mu_1(K_1a_1-K_1d-a_1d)}{K_1a_1}+\frac{r\mu_2(K_2a_2-K_2d-a_2-d)}{K_2a_2}>0.$

Assume that $\alpha_1\alpha_2<0$and $\alpha_1+\alpha_2>0$. Then the real parts of all solutions of $H(\lambda)$ can be still negative if $\alpha_1\alpha_2+d(\beta_1+\beta_2)>0$ and $d(\alpha_1\beta_2+\alpha_2\beta_1)$. By algebraic calculations, we can conclude that if $\alpha_i<0$ and the dispersal of predator $y_i$ is large enough, then the interior equilibrium $(\mu_1,\nu_1,\mu_2,\nu_2)$ can still be locally asymptotically stable, where $\rho_i$ should satisfy the following condition:

$ \rho_i>\max\Big\{\frac{-\nu_j-r_j\mu_i\mu_j(K_ia_i-K_id-a_i-d)(K_ja_j-K_jd-a_j-d)}{(K_iK_ja_j\nu_jd\nu_i(a_i-d)^2)}, $

$ \frac{-\frac{\mu_i\nu_jK_j(\nu_i\rho_j+1)(a_j-d)^2(K_ia_i-K_id-a_i-d)}{r_j\mu_j\nu_iK_i(a_i-d)^2(K_ja_j-K_jd-a_j-d)}-1}{\nu_j}\Big\} . $

Now if $a_1=a_2=a, r_1=r_2=1, K_1=K_2=K, d_1=d_2=d$. The discussions above implies that Model (4) has the same stability at $E^i=(\mu,\nu,\mu,\nu)$ as the stability of the single patch model (5) at $(\mu, \nu)$ where $\mu=\frac{d}{a-d}$ and $\nu=\frac{(K-\mu)(1+\mu)}{aK}$. Therefore, $E^i$ is locally asymptotically stable if $\frac{K-1}{2}<\mu<K$ while it is unstable if $\mu>\frac{K-1}{2}$.

Now we should show that Model (4) has the unique $E^i=(\mu,\nu,\mu,\nu)$ whenever $a>d$. Notice that $F(x_2)$ and $G(x_1)$ have the following properties in the symmetric case (i.e., $a_1=a_2=a, r_1=r_2=1, K_1=K_2=K, d_1=d_2=d$):

1.$F(x_2)=\frac{\rho_1x_2\left(K-x_2\right)+Kd}{\rho_1(K-x_2)+aK(1-d)}$ with $F(0)=\frac{d}{\rho_1+a -d}$ and $F(\mu)=F(K)=\mu$. In addition, $F(x_2)>0$ for all $x_2\in [0,K]$ and $F(x_2)$ has a unique critical point $x^c_2$ for $x_2\in [0,K]$ where

$\mu<x^c_2=\frac{K\left(\rho_1+a-d-\sqrt{(a-d)(\rho_1+a-d)}\right)}{\rho_1} \in (0,K).$

2.$G(x_1)=\frac{\rho_2x_1\left(K-x_1\right)+Kd}{\rho_2(K-x_1)+aK(1-d)}$ with $G(0)=\frac{d}{\rho_2+a -d}$ and $G(\mu)=G(K)=\mu$. In addition, $G(x_1)>0$ for all $x_1\in [0,K]$ and $G(x_1)$ has a unique critical point $x^c_1$ for $x_1\in [0,K]$ where

$\mu<x^c_1=\frac{K\left(\rho_2+a-d-\sqrt{(a-d)(\rho_2+a-d)}\right)}{\rho_2} \in (0,K).$

The discussions above indicate that both $F(x_2)$ and $G(x_1)$ are unimode maps with a unique interception at $x_1=x_2=\mu$.

Proof of Theorem 4.1.

Proof.Proof of Item 1 can be obtained by adopting the proof provided in Theorem 3.1. We omit details.

The stability of $E_{0000},\,E_{K_1000},\,E_{00K_20},\,E_{K_10K_20}$ can be obtained from eigenvalues of the Jacobian matrix of Model (12) evaluated at these equilibria through simple algebraic calculations. We omit details. But we will return to the local stability of $E^b_i$ when we prove Item 4.

Item 3(a): If $x_i=0$, then we have $\lim_{t\rightarrow\infty} y_1(t)=\lim_{t\rightarrow\infty} y_2(t)=0$. This implies that the omega limit set of Model (13) is $y_1=y_2=0$. Therefore, prey $x_i$ persists by applying Theorem 2.5 of [33] since

$\frac{dx_i}{x_idt}\Big\vert_{x_i=0}=r_i>0.$

Item 3(b): Recall $p_i(x)=\frac{a_i x}{1+x}$ and $q_i(x)=\frac{r_i(K_i-x)(1+x)}{a_iK_i}$ for $i=1,2$. Then we construct the following Lyapunov function

$ $$\begin{aligned} V_{ij}(x_i,y_i,y_j) =(\rho_j +d_j)\int^{x_i}_{K_i} \frac{p_i(\xi)-p_i(K_i)}{p_i(\xi)}d\xi + (\rho_j +d_j)y_i+\rho_iy_j. \end{aligned}$$ $

If $ \hat{\mu_i}>K_i$, then we have

$\frac{dV_{ij}(x_i,y_i,y_j)}{dt} =(\rho_j+d_j)\left[(p_i(x_i)-p_i(K_i)\right]q_i(x_i)+y_i(p_i(K_i)-\hat{d}_i)<0$

since $p_i(K_i)-\hat{d}_i<0\Leftrightarrow \hat{\mu_i}>K_i.$ Therefore, Model (13) has global stability at $(K_i,0,0)$ if $\hat{\mu_i}>K_i$.

Item 3(c): We construct the following Lyapunov function

$ $$\begin{aligned} V_{ij}(x_i,y_i,y_j) =(\rho_j +d_j)\int^{x_i}_{\hat{\mu}_i} \frac{p_i(\xi)-p_i(\hat{\mu}_i)}{p_i(\xi)}d\xi + (\rho_j +d_j)\int^{y_i}_{\hat{\nu}_i} \frac{\eta_i-\hat{\nu}_i}{\eta_i}d\eta_i+\rho_i\int^{y_j}_{\hat{\nu}_j^i} \frac{\eta_j-\hat{\nu}_j^i}{\eta_j}d\eta_j. \end{aligned}$$ $

If $\frac{K_i-1}{2}< \hat{\mu_i}<K_i$, then we have

$\frac{dV_{ij}(x_i,y_i,y_j)}{dt} $

$ =(\rho_j+d_j)\left[(p_i(x_i)-p_i(\hat{\mu}_i)\right]\left[q_i(x_i)-\hat{\nu}_i\right]-\frac{\rho_i\hat{\nu}_i((\rho_j+d_j)y_j-\rho_jy_i)^2}{(\rho_j+d_j)y_iy_j}<0.$

Therefore, Model (13) has global stability at $(\hat{\mu_i},\hat{\nu_i},\hat{\nu}_j^i)$ if $\frac{K_i-1}{2}<\hat{\mu_i}<K_i.$

Item 4: If $x_j=0$, then Model (12) reduces to Model (13) who has global stability at $(K_i,0,0)$ if $\hat{\mu_i}>K_i$ while has $(\hat{\mu_i},\hat{\nu_i},\hat{\nu}_j^i)$ if $\frac{K_i-1}{2}<\hat{\mu_i}<K_i.$ Therefore, by applying Theorem 2.5 of [33], prey $x_j$ persists if

$\frac{dx_j}{x_jdt}\Big\vert_{x_j=0,y_j=0}=r_j>0 \mbox{ when } \hat{\mu_i}>K_i$

and

$\frac{dx_j}{x_jdt}\Big\vert_{x_j=0,y_j=\hat{\nu}_j^i}=r_j-a_j\hat{\nu}_j^i>0 \mbox{ when } \frac{K_i-1}{2}<\hat{\mu_i}<K_i.$

The persistence of both prey can be easily obtained from the persistence of one prey.

If $\frac{K_i-1}{2}<\hat{\mu_i}<K_i$ and $r_j-a_j\hat{\nu}_j^i<0$, then $\frac{dx_j}{x_jdt}\Big\vert_{x_j=0,y_j=\hat{\nu}_j^i}=r_j-a_j\hat{\nu}_j^i<0$. According to Theorem 2.18 by [34], we can conclude that the boundary equilibrium $(\hat{\mu_i},\hat{\nu_i},0,\hat{\nu}_j^i)$ is locally asymptotically stable. This proves the stability condition of $E^b_i$ of Item 2.

Item 5: We construct the following Lyapunov function

$ $$\begin{aligned} V(x_i,y_i,x_j,y_j) =&(\rho_i +d_i)\int^{x_j}_{\hat{\mu}_j} \frac{p_j(\xi_j)-p_j(\hat{\mu}_j)}{p_j(\xi_j)}d\xi_j + (\rho_i +d_i)\int^{y_j}_{\hat{\nu}_j} \frac{\eta_j-\hat{\nu}_j}{\eta_j}d\eta_j\\&+\rho_i x_i+\int^{y_i}_{\hat{\nu}_i^j} \frac{\eta_i-\hat{\nu}_i^j}{\eta_i}d\eta_i. \end{aligned}$$ $

Then we have

$ $$\begin{array}{l} &\frac{dV(x_i,y_i,x_j,y_j)}{dt} \\ =&(\rho_j+d_j)\left[(p_i(x_i)-p_i(\hat{\mu}_i)\right]\left[q_i(x_i)-\hat{\nu}_i\right]-\frac{\rho_i\hat{\nu}_i((\rho_j+d_j)y_j-\rho_jy_i)^2}{(\rho_j+d_j)y_iy_j}\\&+\rho_jp_i(x_i)\left[q_i(x_i)-\hat{\nu}_i^j\right]\\ <&(\rho_j+d_j)\left[(p_i(x_i)-p_i(\hat{\mu}_i)\right]\left[q_i(x_i)-\hat{\nu}_i\right]-\frac{\rho_i\hat{\nu}_i((\rho_j+d_j)y_j-\rho_jy_i)^2}{(\rho_j+d_j)y_iy_j}\\&+\rho_jp_i(x_i)\left[q_i(\frac{K_i-1}{2})-\hat{\nu}_i^j\right]\\ =&(\rho_j+d_j)\left[(p_i(x_i)-p_i(\hat{\mu}_i)\right]\left[q_i(x_i)-\hat{\nu}_i\right]-\frac{\rho_i\hat{\nu}_i((\rho_j+d_j)y_j-\rho_jy_i)^2}{(\rho_j+d_j)y_iy_j}\\&+\rho_jp_i(x_i)\left[\frac{r_i (K_i+1)^2}{4a_iK_i}-\hat{\nu}_i^j\right]\end{array}$$ .$

Therefore, if $\frac{K_j-1}{2}<\hat{\mu_j}<K_j \;and\; \frac{r_i (K_i+1)^2}{4a_iK_i}<\hat{\nu}_j^i$ hold, then we have $\frac{dV(x_i,y_i,x_j,y_j)}{dt}<0$. This implies that Model (13) has global stability at $(0,\hat{\nu}_i^j,\hat{\mu_j},\hat{\nu_j})$.

Item 6: Define $V(y_1,y_2)=\rho_2y_1+\rho_1y_2$, then we have

$\frac{dV}{dt}=\rho_2(p_1(x_1)-d_1)y_1+\rho_1(p_2(x_2)-d_2)y_2.$

Notice that $\limsup_{t\rightarrow\infty} x_i(t)\leq K_i$ for both $i=1,2$. Then if $\mu_i>K_i$ for both $i=1,2$, then we have $\max\{p_1(K_1)-d_1,p_2(K_2)-d_2\}<-\delta<0$. This implies that

$\frac{dV}{dt}=\rho_2(p_1(x_1)-d_1)y_1+\rho_1(p_2(x_2)-d_2)y_2<-\delta(\rho_2y_1+\rho_1y_2).$

Therefore, both predators go extinct if $\mu_i>K_i$ for both $i=1$ and $i=2$. Since both $\limsup_{t\rightarrow\infty} y_i(t)=0$ for both $i=1,2$. Then we have Model (12) reduced to the following uncoupled prey model

$x_i'=r_i x_i\left(1-\frac{x_i}{K_i}\right)$

which converges to $x_i=K_i$. Thus, Model (13) has global stability at $(K_1,0,K_2,0)$ when $\mu_i>K_i$ for both $i=1,2$.

On the other hand, if $0<\mu_i<K_i$ for both $i=1,2$, then we have $\max\{p_1(K_1)-d_1,p_2(K_2)-d_2\}>\delta>0$. This implies that

$\frac{dV}{dt}=\rho_2(p_1(x_1)-d_1)y_1+\rho_1(p_2(x_2)-d_2)y_2>\delta(\rho_2y_1+\rho_1y_2).$

Therefore, both predators persist if $0<\mu_i<K_i$ for both $i=1$ and $i=2$.

Item 7 can be obtained from Item 4 and Item 6.

Item 8 can be obtained from eigenvalues of the Jacobian matrix of Model (12) evaluated at the symmetric interior equilibrium $(\mu,\nu,\mu,\nu)$ through simple algebraic calculations. We omit details. The global stability of $(\mu,\nu,\mu,\nu)$ when $\frac{K-1}{2}<\mu<K$ can be obtained by constructing the following Lyapunov function

$ $$\begin{aligned} V(x_1,y_1,x_2,y_2) =&\rho_2\int^{x_1}_{\mu} \frac{p_1(\xi_1)-p_1(\mu)}{p_1(\xi_1)}d\xi_1 + \rho_2\int^{y_1}_{\nu} \frac{\eta_1-\nu}{\eta_1}d\eta_1\\&+\rho_1\int^{x_2}_{\mu} \frac{p_2(\xi_2)-p_2(\mu)}{p_2(\xi_2)}d\xi_2+\rho_1\int^{y_2}_{\nu} \frac{\eta_2-\nu}{\eta_2}d\eta_2\end{aligned}$$ $

which gives

$\frac{dV(x_1,y_1,x_2,y_2)}{dt} =\rho_2(p_1(\xi_1)-p_1(\mu))(q_1(x_1)-K)+\rho_1(p_2(\xi_2)-p_2(\mu))(q_2(x_2)-K)$

$+\frac{\rho_1\rho_2\nu(y_1-y_2)\left(\frac{1}{y_1}-\frac{1}{y_2}\right)}{y_1y_2}.$

---

Acknowledgments

This research is partially supported by NSF-DMS(1313312); NSF-IOS/DMS (1558127) and The James S. McDonnell Foundation 21st Century Science Initiative in Studying Complex Systems Scholar Award (UHC Scholar Award 220020472). We also would like to acknowledge the partial support from Natural Science Foundation of Jiangsu Province (BK20140927), Tianyuan Fund for Mathematics of NSFC (11426132, 11601226), and from Nanjing Tech University. The research of K.M is also partially supported by the Department of Education GAANN (P200A120192). S.K.S. is partially supported by the NBHM post-doctoral fellowship. All authors would like to thank Dr. Andrea Bruder for the discussions on the modeling dispersal strategies in the early stage of this manuscript.

---