Predation-induced dispersal toward fitness for predator invasion in predator–prey models

Abstract

In this paper, we consider a predator–prey model with nonuniform predator dispersal, called predation-induced dispersal (PID), which represents predator motility depending on the maximal predation rate and the predator death rate in a spatially heterogeneous region. We study the local stability of the semitrivial steady state when predators are absent for models with PID and linear dispersal. We then investigate the local/global bifurcation from the semitrivial steady state of these models. Finally, we compare the results of the model with PID to the results of the model with linear dispersal. We conclude that the nonuniform dispersal of predators obeying PID increases fitness for predator invasion when rare; thus, predators with PID can invade a region with an increased probability even in cases wherein predators dispersed linearly cannot invade a certain region. Based on the results, we provide an ecological interpretation with the simulations.

Keywords:

• Predation-induced diffusion
• evolution of dispersal
• local stability
• global bifurcation

AMS Classification::

• 35K57
• 92D25
• 92D40

1. Introduction

In ecology, species dispersal can occur by various factors such as species density in its habitat, its response to the environment, and interactions with other species causing migration to other habitats. Therefore, a species' response to its habitat environment and interactions with other species are important factors in modelling a more natural dispersal [45,50]. Many authors have studied population models with the evolution of dispersal, depending on other species, represented by reaction–diffusion equations: self and cross-diffusion models [34,39–41] and the prey-taxis model [7,20,25,37,51,53]. Additionally, considerable studies have been performed on models with starvation-driven diffusion that depend on resources [6,8,10,11,30–32]. Migration of species affords such as gaining better-quality resources, fleeing disadvantageous areas, and evading hostile species, which are crucial for species survival. In terms of evolution, the dispersion principle according to environmental conditions suggests how a species evolves to adapt to the environment of the species' habitat and whether this dispersal way provides an evolutionary advantage for the species' survival. Considerable studies have been performed on species dispersal from an evolutionary perspective by establishing discrete and continuous mathematical models [4,12,26,35,36,43–45].

In a heterogeneous habitat, rapid diffusion of species gives a disadvantage to its survival when its dispersal is spatially uniform [17]. In other words, the dispersal of biological organisms in a heterogeneous environment hinders their survival. Thus when establishing a biological model with a spatially heterogeneous habitat, uniform diffusion is inappropriate to represent the biological dispersal of a species. Due to this disadvantageous effect of diffusion to species, it is difficult to construct a proper model that describes the dispersal behaviour of species in habitats with spatial heterogeneity. Particularly in the predator–prey interactions, the effect of dispersal on the environmental conditions in the local habitat has not been well studied. Thus, it is worthwhile to study nonuniform dispersal depending on the location in the habitat to examine which dispersal is evolutionarily advantageous for biological predator–prey species.

Among the main reasons for the dispersal of species is to escape unfavorable regions with food shortages, harsh climate conditions, or pressure from predators for survival. Since species migration is unprofitable to them if the habitat conditions are uniform spatially and temporally, the habitat environment should be considered spatially and temporally heterogeneous. A low dispersal rate benefits species survival when the environment is spatially nonuniform [16,17,24,26,28,42]. On the other hand, a high dispersal rate favors species survival if the environment is changeable over time [16,24,26,42,52]. For a habitat with spatial and temporal heterogeneity, interaction with other species in a region determines the optimal dispersal rate for species survival [16,24]. In previous studies, it has been assumed that species have uniform and random movement, regardless of how the environmental conditions are given. Meanwhile, for the predator–prey model, various types of functional responses to predation have been examined extensively. Many models have a functional response in which ignoring an effect of predator density, i.e. the function that describes the density of prey consumed by its predator depends only on the prey. However, considerable explicit biological and physiological evidence [1,3,13,18] has shown that in various situations, when predators have to search for food in heterogeneous situations, they exhibit a ratio-dependent functional response in which the per capita predator growth rate should be a function of the prey-to-predator abundance ratio [2,4,21–23,27,33,46,48].

In this work, it is assumed that a species tends to stay in a current location when such a location is favorable for survival; thus, the species will make its motility low, as will the dispersal rate. In contrast, if the species is located in an unfavorable region, its motility will be high to escape such a region. In other words, the species will move with rapid dispersal. This movement strategy indicates that even though a species usually does not have prior information on a habitat condition, such motility change generally induces the species to move toward a better environment.

To describe such a phenomenon, we introduce the predator dispersal strategy, predation-induced dispersal (PID), which represents the change in motility of a predator depending on the predation rate and predator death rate in a region. More precisely, if predators are located in a low-predation region wherein prey can easily hide and escape from predators or in a dangerous place where the predator death rate is high, they have to migrate to other regions. In other words, when the predation is low or the predator death rate is high, they disperse faster and more actively and tend to avoid such locations. In contrast, if the predation is high or the predator death rate is low, the predator dispersal is small because the predators are satisfied with their hunting and tend to remain in that region.

In this paper, we suggest a dispersal strategy that cooperates with fitness for predator invasion in a heterogeneous habitat. To achieve our goal, we examine the predator–prey model with ratio-dependent functional responses in which the predator's movement follows the rule of PID to understand how PID affects the fitness for a species invasion. First, we investigate the local stability of the semitrivial steady state when predators are absent for models with PID and linear dispersal. We then study the local/global bifurcation from the semitrivial steady state of these models. Finally, we compare the results of the model with PID to those with linear dispersal. We conclude that the nonuniform dispersal of predators obeying PID increases fitness for predator invasion when rare; thus, predators with PID can invade a region with an increased probability even in cases wherein predators dispersed linearly cannot invade a certain region.

The outline of this article is as follows. In Section 2, we set up the mathematical model in which predators have PID. The stability of the semitrivial steady state of the models with linear dispersal and PID is investigated with ecological interpretations in Section 3. Section 4 is devoted to obtaining the result of the coexistence state of the PID model by using the linearized eigenvalue and local/global bifurcation theory. In particular, the results are compared with those of a model based on linear dispersal. In Section 5, we summarize the obtained results and discuss the biological interpretations with numerical simulations.

2. Mathematical models with PID

In this section, we present a predator–prey model with PID in a spatially heterogeneous environment under no flux boundary conditions. First, we introduce the PID strategy mathematically. Let β and c be the predation rate and predator death rate, respectively; $\text{β}/m$ is the maximal predation rate, where m is the half-capturing saturation parameter. Denote a constant i as the measure of satisfaction in the environment defined by $i=\frac{mc}{\text{β}}.$It is noteworthy that i increases when the maximal predation is low, and i decreases when the death rate is low. That means that i indicates whether a given environment is good for predator survival. If a predator's habitat is spatially heterogeneous where the predation rate $\text{β}\left(x\right)$ depends on a location x, satisfaction measure i can be extended to function $i\left(x\right)$ with location variable x as $i\left(x\right)=\frac{mc}{\text{β}\left(x\right)}.$Then, we consider a mobility function $p\left(i\right(x\left)\right)$, which is an increasing function with respect to i. The increasing property of $p\left(i\right(x\left)\right)$ represents that predators tend to quickly leave the location when they are in an area where $i\left(x\right)$ is large and stay in the location when they are in an area where $i\left(x\right)$ is small. At the individual level, the species with the mobility trait $p\left(i\right(x\left)\right)$ decides its moving rate only depending on local habitat conditions, not having a purpose in where to go. Thus, such a decision does not guarantee locating a better place because it does not depend on the arrival point. The dispersal of predator populations with such traits can be represented by $\mathrm{\Delta }\left(p\right(i\left(x\right)\left)v\right)$. When predators are satisfied with the environment, they reduce their mobility to stay in the region; conversely, when they are in a region where there is not sufficient food to survive, they increase their mobility to avoid such unfavorable locations.

For example, if we consider $p\left(i\right(x\left)\right):=i(x{)}^k={\left(\frac{mc}{\text{β}\left(x\right)}\right)}^k$for some positive constant k, then, $\mathrm{\Delta }\left(p\right(i\left(x\right)\left)v\right)=\mathrm{\nabla }\cdot p\left(i\right(x\left)\right)(\mathrm{\nabla }v-kv\mathrm{\nabla }(\ln \text{β}\left(x\right)),$which is diffusion with a logarithmic form of advection that represents predators moving toward a high-predation habitat [29].

With the preceding motivations, in this study, we investigate the following predator–prey interaction systems with PID and ratio-dependent functional responses in a spatially heterogeneous environment: (1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) where $\mathrm{\Omega }\subset {\mathbb{R}}^N$ is a bounded domain and the zero-flux condition is given on smooth boundary $\mathrm{\partial }\mathrm{\Omega }$, where $\overrightarrow{n}$ is the outward unit normal vector on the boundary. The functions $u(x,t)$ and $v(x,t)$ represent the prey and predator population densities, respectively, at location x and time t. The predator's reaction term indicates that the predator is fed by only one prey and not by any other resource. Thus, the model describes the interaction by which the predator disperses, $p\left(i\right)$, which obeys the PID rule defined below; however, the prey diffuses linearly at a constant rate. Moreover, $K\left(x\right)$, $\frac{\alpha \left(x\right)}{m}$, and $\frac{\text{β}\left(x\right)}{m}$ are positive functions on $x\in \mathrm{\Omega }$ and represent the carrying capacity, maximal capturing rate by the predator, and maximal predation rate, respectively, in a spatially heterogeneous environment; c,r and m are positive constants that represent the death rate of the predator, the intrinsic growth rate of the prey, and the half-capturing saturation constant, respectively.

In this work, we consider PID under the following assumptions: (G) $\begin{array}{rl}& p\left(i\right)>0\text{is an increasing and bounded smooth function for}i\ge 0.\\ & p\left(i\right)+i(1-i)p^{\prime }\left(i\right)>0\text{ for }i>1.\end{array}$(G) Assumption (S) will be used for local/global bifurcation of the model. We denote $\underset{i>0}{max}p\left(i\right)=h$ and $\underset{i>0}{min}p\left(i\right)=\ell$.

Here is an example of $p\left(i\right)$ satisfying (G) and (S). If we consider a step function (Figure 1(a)), as in the simple case of $p^0\left(i\right)$ defined by $p^0\left(i\right)=\{\begin{array}{ll}\ell & \text{for}0<i\le 1,\\ h& \text{for}1<i<\mathrm{\infty },\end{array}$dispersal function $p^0$ can be approximated by a smooth motility function defined by convolution, (2) $p^{ϵ}:=p^0\ast {\eta }^{ϵ},$(2) where $ϵ>0$ is small and ${\eta }^{ϵ}$ is a smooth symmetric mollifier with $\mathrm{supp}{\eta }^{ϵ}\subset (0,ϵ)$ and $\int {\eta }^{ϵ}\left(x\right)\mathrm{d}x=1$. Then, $p^{ϵ}\left(i\right)$ is a smooth function, and $p^{ϵ}$ converges to $p^0$ as ε approaches 0. Additionally, $p^{ϵ}$ satisfies (??) and (??) for a sufficiently small $ϵ>0$.

Figure 1. The shapes of $p\left(i\right)$. (a) A step function case $p\left(i\right)=p^0\left(i\right)$. (b) $p\left(i\right)=\nu \delta \left(i\right)$ with bounds ℓ and h. $p\left(i\right)$ is not necessarily symmetric with $p\left(1\right)=\nu$.

Now, we introduce a mobility function, $p\left(i\right)$, with a shape that satisfies (G) and (S). For the given mobility function $p\left(i\right)$, we determine $\delta \left(i\right)$ such that: $p\left(i\right(x\left)\right)=p\left(1\right)\delta \left(i\right(x\left)\right),$where $\delta \left(i\right)$ satisfies (G), (S), and $\delta \left(1\right)=1$. We call $\delta \left(i\right)$ the shape of $p\left(i\right)$.

If $p\left(i\right)\equiv \nu$ for a constant ν, predator dispersal is not affected by the surrounding environment. Since $\frac{\text{β}\left(x\right)}{m}-c=0$ for i = 1, the equality implies that the net growth rate is given by: $v\left(\text{β}\right(x)\frac{u}{mu+v}-c)\approx 0,$if the prey population is sufficiently large. In this case, the predator need not be concerned about the environment at its location, so it exhibits natural dispersal, ν, independent of the environment at i = 1. Thus, we set $p\left(1\right)=\nu$ to compare the results obtained from the conditions for the PID model to those with a constant dispersal. Therefore, PID is defined as $p\left(i\right)=\nu \delta \left(i\right)$ (Figure 2(b)). Now, we will consider high- and low- predation regions ${\mathrm{\Omega }}_H$, and ${\mathrm{\Omega }}_L$, respectively, defined as: ${\mathrm{\Omega }}_H:=\{x\in \mathrm{\Omega }:\frac{\text{β}\left(x\right)}{m}>c\};{\mathrm{\Omega }}_L:=\{x\in \mathrm{\Omega }:\frac{\text{β}\left(x\right)}{m}<c\}.$The criterion for dividing ${\mathrm{\Omega }}_H$ and ${\mathrm{\Omega }}_L$ is that the maximal predation rate is equal to the death rate, $\frac{\text{β}\left(x\right)}{m}=c$. For expansion, the domain types are defined as follows: Ω is called a high (low) predation domain if ${\int }_{\mathrm{\Omega }}(\frac{\text{β}\left(x\right)}{m}-c)>(<)0.$Note that if $\mathrm{\Omega }={\mathrm{\Omega }}_L$, Ω is an unlivable domain.

Figure 2. Stability region of $({\theta }_{\mu },0)$ with respect to m and ν for linear dispersal.

The spatial heterogeneity of the environment is represented by a nonconstant function $\text{β}\left(x\right)$ and such spatial heterogeneity is crucial to the stability analysis in this study.

We also consider a predator–prey model with uniform diffusion ($p\left(i\right)=\nu$) in system (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ), where ν is a positive constant: (3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) We then analyze the local stability of the semitrivial steady state of (3(3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) ) and determine the condition wherein (3(3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) ) has a coexistence steady state. The effect of regional heterogeneity and PID on the stability of the semitrivial steady states of model (3(3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) ) is investigated. Furthermore, we compare the stability conditions of the model incorporating PID with those of a model with uniform diffusion. In addition, the coexistence steady state of the models is examined for two different dispersal strategies: PID and linear dispersal.

3. Stability of the semitrivial steady state

We consider a predator–prey interaction system (3(3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) ) with linear dispersal represented by a constant diffusion and system (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ) with PID. To investigate invasion by predators, we analyze the stability of the semitrivial state $({\theta }_{\mu },0)$ and compare the conditions for these two models to examine the effect of PID on fitness for species invasion when rare.

First, we consider the following scalar equation: (4) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(4) Let $f(x,u)=r(1-\frac{u}{K\left(x\right)})$. Then, $f_u=-r\frac{1}{K\left(x\right)}<0$, and $uf(x,0)=ru$. Since the principal eigenvalue ${\lambda }_1(\mu \mathrm{\Delta }+r)$ is positive, the single Equation (4(4) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(4) ) has a unique positive solution denoted by ${\theta }_{\mu }$ (see [5]). Therefore, when predators are absent, system (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ) (also model (3(3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) )) has a unique semitrivial steady state.

Since initial data $u(x,0)$ and $v(x,0)$ are nonnegative and not identically zero, $u(x,t)>0$ and $v(x,t)>0$ for all $x\in \overline{\mathrm{\Omega }}$ and all t>0 by the maximum principle [47]. Moreover, $\left(u\right(x,t),v(x,t\left)\right)$ are the classical solutions of systems (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ) and (3(3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) ) and exist for all times t>0. A positive solution $\left(u\right(x,t),v(x,t\left)\right)$ of system (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ) (or (3(3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) )), which is $u(x,t)>0$, $v(x,t)>0$ is called the coexistence of system (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ) (or (3(3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) )). In addition, we say that a positive steady state $\left(u\right(x),v(x\left)\right)$ of system (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ) (or (3(3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) )) is the coexistence state when both components are positive; it is a semitrivial state when one component is positive and the other is zero. Note that in systems (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ) and (3(3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) ), there is only one semitrivial steady state $({\theta }_{\mu },0)$.

It is also noteworthy that the prey in system (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ) or (3(3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) ) always survives if a certain condition is satisfied.

Theorem 3.1

Let $\left(u\right(x,t),v(x,t\left)\right)$ be a solution of (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ) or (3(3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) ). If $r-\alpha \left(x\right)>0$ for all $x\in \mathrm{\Omega }$, then $u(x,t)>0$ and $\underset{t\to \mathrm{\infty }}{lim inf}u(x,t)>0$.

Proof.

From the first equation of (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ), we have (5) $u_t=\mu \mathrm{\Delta }u+\frac{ru}{K\left(x\right)}\left(K\right(x)-u-\frac{\alpha \left(x\right)K\left(x\right)v}{r(mu+v)})\ge \mu \mathrm{\Delta }u+\frac{ru}{K\left(x\right)}\left(K\right(x)(1-\frac{\alpha \left(x\right)}{r})-u).$(5) Let $\tilde{u}(x,t)$ be a solution of the equation: (6) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+\frac{ru}{K\left(x\right)}\left(K\right(x)(1-\frac{\alpha \left(x\right)}{r})-u)}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(6) From [15], (6(6) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+\frac{ru}{K\left(x\right)}\left(K\right(x)(1-\frac{\alpha \left(x\right)}{r})-u)}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(6) ) has a unique positive steady-state solution $\tilde{u}\left(x\right)$ which is globally asymptotically stable. By the comparison principle and (5(5) $u_t=\mu \mathrm{\Delta }u+\frac{ru}{K\left(x\right)}\left(K\right(x)-u-\frac{\alpha \left(x\right)K\left(x\right)v}{r(mu+v)})\ge \mu \mathrm{\Delta }u+\frac{ru}{K\left(x\right)}\left(K\right(x)(1-\frac{\alpha \left(x\right)}{r})-u).$(5) ), $u(x,t)\ge \tilde{u}(x,t)$ for $(x,t)\in \mathrm{\Omega }\times (0,\mathrm{\infty })$. Additionally, we have $\underset{t\to \mathrm{\infty }}{lim inf}u(x,t)\ge \underset{t\to \mathrm{\infty }}{lim inf}\tilde{u}(x,t)=\tilde{u}\left(x\right)>0.$

3.1. Model with linear dispersal

In this subsection, we present the results of the linear stability for system (3(3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) ). Since (3(3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) ) is a special case wherein $p\left(i\right)=\nu$ in (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ), we omit the proofs here. We provide proofs for the results of the PID system (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ) in Section 3.2.

In this subsection, we state the results of the stability of the semitrivial steady state. Linearizing system (3(3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) ) at $({\theta }_{\mu },0)$, we obtain the following eigenvalue problem: (7) $\{\begin{array}{ll}{\displaystyle \lambda \varphi =\mu \mathrm{\Delta }\varphi +r(1-\frac{2{\theta }_{\mu }}{K\left(x\right)})\varphi -\frac{\alpha \left(x\right)}{m}\psi }& \\ {\displaystyle \lambda \psi =\nu \mathrm{\Delta }\psi +(\frac{\text{β}\left(x\right)}{m}-c)\psi }& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }\varphi \cdot \overrightarrow{n}=\mathrm{\nabla }\psi \cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(7) For the linear stability of $({\theta }_{\mu },0)$, we have the following:

Lemma 3.2

If ${\lambda }_1(\nu \mathrm{\Delta }+\frac{\text{β}\left(x\right)}{m}-c)>0$, then $({\theta }_{\mu },0)$ is linearly unstable, and if ${\lambda }_1(\nu \mathrm{\Delta }+\frac{\text{β}\left(x\right)}{m}-c)<0$, then $({\theta }_{\mu },0)$ is linearly asymptotically stable.

Denote the average of function $\text{β}\left(x\right)$ by $\overline{\text{β}}$: $\overline{\text{β}}:=\frac{1}{|\mathrm{\Omega }|}{\int }_{\mathrm{\Omega }}\text{β}\left(x\right)\mathrm{d}x.$From Lemma 3.2, we have the stability of $({\theta }_{\mu },0)$.

Theorem 3.3

Let $({\theta }_{\mu },0)$ be the semitrivial steady-state solution of (3(3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) ).

(i)	Suppose that $mc<\overline{\text{β}}$. Then $({\theta }_{\mu },0)$ is linearly unstable for any $\nu >0$.
(ii)	Let M be the positive constant satisfying $\underset{x\in \mathrm{\Omega }}{max}\text{β}\left(x\right)=Mc$ for a given $\text{β},c$. Then, for a given $m\in (\overline{\text{β}}/c,M)$, there exists ${\nu }^{\ast }$ depending on $\mathrm{\Omega },m,\text{β}$, and c such that if $\nu >{\nu }^{\ast }$, then $({\theta }_{\mu },0)$ is linearly asymptotically stable, and if $\nu <{\nu }^{\ast }$, then $({\theta }_{\mu },0)$ is linearly unstable.
(iii)	For a given ν, there exists $m_{\ast }\in (\overline{\text{β}}/c,M)$ depending on $\mathrm{\Omega },\nu ,\text{β}$ and c such that if $m<m_{\ast }$, then $({\theta }_{\mu },0)$ is linearly unstable, and if $m>m_{\ast }$, then $({\theta }_{\mu },0)$ is linearly asymptotically stable.

According to the above theorem, if either the predator death rate or the half-capturing saturation constant is sufficiently small, predators can invade a region when rare. For a certain range of the half-capturing saturation constant, the threshold value relates to the dispersal rate of the predator and determines the stability of the semitrivial steady state where the predator is absent. In addition, for a given predator dispersal rate, the threshold value relates to the half-capturing saturation constant that determines the stability of the semitrivial steady state where there is no predator. Figure 2 represents such phenomena.

From Theorems 3.1 and 3.3, it is noteworthy that the instability of the semitrivial steady state of (3(3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) ) and positivity of $u(x,t)$ imply the coexistence steady states of system (3(3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) ). This can be shown using the fixed-point index theory as in [10,29]. Instead, we will provide the coexistence results by using the local/global bifurcations in the next section.

Corollary 3.4

Assume that $r-\alpha \left(x\right)>0$ for $x\in \mathrm{\Omega }$.

(i)	If $mc<\overline{\text{β}}$, then (3(3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) ) has the coexistence for any $\nu >0$.
(ii)	Let M be a positive constant satisfying $\underset{x\in \mathrm{\Omega }}{max}\text{β}\left(x\right)=Mc$ for a given $\text{β},c$. Then, for a given $m\in (\overline{\text{β}}/c,M)$, there exists ${\nu }^{\ast }$ depending on $\mathrm{\Omega },m,\text{β}$ and c such that if $\nu >{\nu }^{\ast }$, then (3(3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) ) has the coexistence.
(iii)	For a given ν, there exists $m_{\ast }\in (\overline{\text{β}}/c,M)$ such that if $m<m_{\ast }$, then (3(3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) ) has the coexistence.

3.2. Model with PID

Recall that PID is a dispersal in which a predator's movement is affected by the predation rate in the circumstances of the location in the habitat. More precisely, if a predator favors a certain region since its predation rate is high (this implies that there are good enough prey in the region to satisfy the predator), it tends to stay there. Thus, its motility is low, and it disperses slowly. However, if a predator feels there is a lack of prey, it will avoid that region, and its motility will become high. To represent this phenomenon, PID is defined in Section 1 as follows: (8) $p\left(i\right)=\nu \delta \left(i\right),$(8) where $\delta \left(i\right)$ is the shape of mobility function $p\left(i\right)$, satisfying (G), (S), and $\delta \left(1\right)=1$. Here, ν represents the predator dispersal rate independent of the given environmental elements $\text{β},m$ and c. The setting (8(8) $p\left(i\right)=\nu \delta \left(i\right),$(8) ) will be used to compare PID with linear dispersal.

In this subsection, we investigate the stability of the semitrivial steady-state solution $({\theta }_{\mu },0)$ of (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ).

Letting $V:=p\left(i\right)v$, (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ) is equivalent to (9) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle V_t=p\left(i\right)[\mathrm{\Delta }V+\frac{V}{p\left(i\right)}\left(\text{β}\right(x)\frac{u}{mu+v}-c)]}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }V\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,V(x,0)=V_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(9) where $V_0=p\left(i\right)v_0$.

Denote $g(x,u,V)=ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}$and $h(x,u,V)=\frac{V}{p\left(i\right)}\left(\text{β}\right(x)\frac{u}{mu+v}-c).$By simple calculation, we obtain $\left(\begin{array}{cc}{g}_u& g_V\\ h_u& h_V\end{array}\right)=\left(\begin{array}{cc}r(1-{\displaystyle \frac{2u}{K\left(x\right)}})-{\displaystyle \frac{\alpha \left(x\right)v^2}{(mu+v{)}^2}}& -{\displaystyle \frac{1}{p\left(i\right)}}{\displaystyle \frac{\alpha \left(x\right)m{u}^2}{(mu+v{)}^2}}\\ {\displaystyle \frac{\text{β}\left(x\right)v^2}{(mu+v{)}^2}}& {\displaystyle \frac{1}{p\left(i\right)}}({\displaystyle \frac{\text{β}\left(x\right)m{u}^2}{(mu+v{)}^2}}-c)\end{array}\right),$and the linearized eigenvalue problem at $({\theta }_{\mu },0)$: (10) $\{\begin{array}{ll}{\displaystyle \lambda \varphi =\mu \mathrm{\Delta }\varphi +r(1-\frac{2{\theta }_{\mu }}{K\left(x\right)})\varphi -\frac{\alpha \left(x\right)}{p\left(i\right)m}\psi }& \\ {\displaystyle \lambda \psi =p\left(i\right)[\mathrm{\Delta }\psi +\frac{1}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c)\psi ]}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }\varphi \cdot \overrightarrow{n}=\mathrm{\nabla }\psi \cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(10) For the linear stability of $({\theta }_{\mu },0)$, we introduce the lemma, which can be proven similarly to [10,31]:

Lemma 3.5

If ${\lambda }_1(\mathrm{\Delta }+\frac{1}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c\left)\right)>0$, then $({\theta }_{\mu },0)$ is linearly unstable, and if ${\lambda }_1(\mathrm{\Delta }+\frac{1}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c\left)\right)<0$, then $({\theta }_{\mu },0)$ is linearly asymptotically stable.

Now, we define (11) $R:=\frac{{\int }_{{\mathrm{\Omega }}_L}[c-\text{β}(x\left)/m\right]}{{\int }_{{\mathrm{\Omega }}_H}\left[\text{β}\right(x)/m-c]}$(11) for a given $\text{β}\left(x\right)$, c, and m. This ratio R indicates the degree of predation in the predator's living environment by considering the difference between the maximal predation and the predator death rate. For example, if R<1, then the region that represents a strong hunting case is more occupied in the entire domain Ω than a region wherein it is difficult for the predator to hunt prey.

Denote the average of the function $\text{β}\left(x\right)$ by $\overline{\text{β}}$: $\overline{\text{β}}:=\frac{1}{|\mathrm{\Omega }|}{\int }_{\mathrm{\Omega }}\text{β}\left(x\right)\mathrm{d}x.$

Lemma 3.6

Let $\text{β}\left(x\right),m,c$ and $p\left(i\right)$ be given.

(i)	If $\overline{\text{β}}-mc>0$, then ${\int }_{\mathrm{\Omega }}\frac{1}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c)\mathrm{d}x>0$,
(ii)	If $\overline{\text{β}}-mc<0$, there exists $\rho =\rho (p,m,\text{β},c)\in (1,\frac{h}{\ell }]$ such that if $\rho >R$, then ${\int }_{\mathrm{\Omega }}\frac{1}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c)\mathrm{d}x>0$, and if $\rho <R$, then ${\int }_{\mathrm{\Omega }}\frac{1}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c)\mathrm{d}x<0$.

Proof.

(i) follows from the monotonicity of $p\left(i\right)$.

(ii) First, we decompose the integral as ${\int }_{\mathrm{\Omega }}\frac{1}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c)\mathrm{d}x={\int }_{{\mathrm{\Omega }}_H}\frac{1}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c)\mathrm{d}x+{\int }_{{\mathrm{\Omega }}_L}\frac{1}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c)\mathrm{d}x.$Since $p\left(i\right)$ is an increasing function of i, we have ${\int }_{{\mathrm{\Omega }}_H}\frac{1}{p\left(1\right)}(\frac{\text{β}\left(x\right)}{m}-c)\mathrm{d}x\le {\int }_{{\mathrm{\Omega }}_H}\frac{1}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c)\mathrm{d}x\le {\int }_{{\mathrm{\Omega }}_H}\frac{1}{p\left(0\right)}(\frac{\text{β}\left(x\right)}{m}-c)\mathrm{d}x.$Additionally, the continuity of $p\left(i\right)$ implies that there exists ${\xi }_1\in (0,1)$ such that ${\int }_{{\mathrm{\Omega }}_H}\frac{1}{p\left({\xi }_1\right)}(\frac{\text{β}\left(x\right)}{m}-c)\mathrm{d}x={\int }_{{\mathrm{\Omega }}_H}\frac{1}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c)\mathrm{d}x.$Similarly, we can find ${\xi }_2\in (1,\underset{x\in \mathrm{\Omega }}{max}i(x\left)\right)$ satisfying ${\int }_{{\mathrm{\Omega }}_L}\frac{1}{p\left({\xi }_2\right)}(\frac{\text{β}\left(x\right)}{m}-c)\mathrm{d}x={\int }_{{\mathrm{\Omega }}_L}\frac{1}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c)\mathrm{d}x.$Therefore, ${\int }_{\mathrm{\Omega }}\frac{1}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c)\mathrm{d}x={\int }_{{\mathrm{\Omega }}_H}\frac{1}{p\left({\xi }_1\right)}(\frac{\text{β}\left(x\right)}{m}-c)\mathrm{d}x+{\int }_{{\mathrm{\Omega }}_L}\frac{1}{p\left({\xi }_2\right)}(\frac{\text{β}\left(x\right)}{m}-c)\mathrm{d}x.$Hence, by choosing $\rho =\frac{p\left({\xi }_2\right)}{p\left({\xi }_1\right)}>1$, we get our assertion.

Thus, we obtain the stability of $({\theta }_{\mu },0)$.

Theorem 3.7

Let $({\theta }_{\mu },0)$ be the unique semitrivial steady-state solution of (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ), and ρ be the number in Lemma 3.6(ii).

1. If $mc<\overline{\text{β}}$, then, $({\theta }_{\mu },0)$ is linearly unstable for any $p\left(i\right)>0$.

2. If $\rho >R$, then, $({\theta }_{\mu },0)$ is linearly unstable for any $p\left(i\right)>0$.

3. Suppose that $\rho <R$, and that M and ${\nu }^{\ast }$ are given constants in Theorem 3.3(ii). Then, for a given m<M, there exists $p^{\ast }>{\nu }^{\ast }$ depending on $\mathrm{\Omega },m,\text{β},\delta$, and c such that if $\nu >p^{\ast }$, then $({\theta }_{\mu },0)$ is linearly asymptotically stable, and if $\nu <p^{\ast }$, then $({\theta }_{\mu },0)$ is linearly unstable.

4. Let $m_{\ast }$ be in Theorem 3.3(iii) for a given ν. Then, for a given p, there exists $m^{\ast }\in (m_{\ast },M)$ depending on $\mathrm{\Omega },p,\text{β}$, and c such that if $m<m^{\ast }$, then $({\theta }_{\mu },0)$ is linearly unstable, and if $m>m^{\ast }$, then $({\theta }_{\mu },0)$ is linearly asymptotically stable.

Proof.

By Lemma 3.5, it suffices to check the sign of ${\lambda }_1(\mathrm{\Delta }+\frac{1}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c))={\lambda }_1(\mathrm{\Delta }+\frac{1}{\nu \delta \left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c\left)\right)$. For simplicity, we denote ${\lambda }_1(\mathrm{\Delta }+\frac{1}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c))$ by ${\tilde{\lambda }}_1$. By Theorem A.2 and Lemma 3.6, we obtain (i) and (ii).

(iii) Since ${\lambda }_1(\nu \mathrm{\Delta }+\frac{1}{\delta \left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c\left)\right)$ and ${\tilde{\lambda }}_1$ have the same sign, the existence of such $p^{\ast }$ follows from Theorem A.2. Now, to show that $p^{\ast }>{\nu }^{\ast }$, it suffices to verify that ${\tilde{\lambda }}_1$ is positive for $\nu ={\nu }^{\ast }$. Let ϕ be a principal eigenfunction corresponding to ${\lambda }_1(\mathrm{\Delta }+\frac{1}{{\nu }^{\ast }}(\frac{\text{β}\left(x\right)}{m}-c\left)\right)=0$ with $||\varphi |{|}_2=1$. Then $\begin{array}{rl}{\tilde{\lambda }}_1& \ge {\int }_{\mathrm{\Omega }}-|\mathrm{\nabla }\varphi {|}^2+\frac{1}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c){\varphi }^2\mathrm{d}x\\ & ={\int }_{\mathrm{\Omega }}-|\mathrm{\nabla }\varphi {|}^2+\frac{1}{{\nu }^{\ast }}(\frac{\text{β}\left(x\right)}{m}-c){\varphi }^2\mathrm{d}x+{\int }_{\mathrm{\Omega }}(\frac{1}{p\left(i\right)}-\frac{1}{{\nu }^{\ast }})(\frac{\text{β}\left(x\right)}{m}-c){\varphi }^2\mathrm{d}x\\ & ={\int }_{\mathrm{\Omega }}(\frac{1}{\nu \delta \left(i\right)}-\frac{1}{{\nu }^{\ast }})(\frac{\text{β}\left(x\right)}{m}-c){\varphi }^2\mathrm{d}x\\ & >{\int }_{{\mathrm{\Omega }}_L}(\frac{1}{{\nu }^{\ast }}-\frac{1}{\nu })(c-\frac{\text{β}\left(x\right)}{m}){\varphi }^2\mathrm{d}x+{\int }_{{\mathrm{\Omega }}_H}(\frac{1}{\nu }-\frac{1}{{\nu }^{\ast }})(\frac{\text{β}\left(x\right)}{m}-c){\varphi }^2\mathrm{d}x=0.\end{array}$(iv) Let ${\psi }_1$ be a principal eigenfunction corresponding to ${\tilde{\lambda }}_1$. For given $\text{β}\left(x\right)$ and c, we choose M such that $\underset{x\in \mathrm{\Omega }}{max}(\frac{\text{β}\left(x\right)}{m}-c)=0$. Then, if m>M, $\frac{\text{β}\left(x\right)}{m}-c<0$ for all $x\in \mathrm{\Omega }$. Thus, we have ${\tilde{\lambda }}_1=\frac{{\int }_{\mathrm{\Omega }}-|\mathrm{\nabla }{\psi }_1{|}^2+\frac{1}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c){\psi }_1^2\mathrm{d}x}{{\int }_{\mathrm{\Omega }}{\psi }_1^2\mathrm{d}x}<0.$For $m>\overline{\text{β}}/c$, ${\lambda }_1(\mathrm{\Delta }+\frac{1}{\nu }(\frac{\text{β}\left(x\right)}{m}-c\left)\right)$ is positive. Observing that $\frac{\text{β}\left(x\right)}{m}-c<\frac{1}{\delta \left(i\right(x\left)\right)}(\frac{\text{β}\left(x\right)}{m}-c)$, we have ${\lambda }_1(\mathrm{\Delta }+\frac{1}{\nu }(\frac{\text{β}\left(x\right)}{m}-c))<{\lambda }_1(\mathrm{\Delta }+\frac{1}{\nu \delta \left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c))$by Lemma A.1, and thus ${\tilde{\lambda }}_1>0$ for $m>\overline{\text{β}}/c$. Since ${\tilde{\lambda }}_1$ is a strictly decreasing function to m, we can find $m^{\ast }\in (\overline{\text{β}}/c,M)$ such that ${\tilde{\lambda }}_1=0$ for $m=m^{\ast }$. Now, we only need to show that $m^{\ast }>m_{\ast }$, where $m_{\ast }$ satisfies ${\lambda }_1(\mathrm{\Delta }+\frac{1}{\nu }(\frac{\text{β}\left(x\right)}{m_{\ast }}-c\left)\right)=0$. Since for $m=m_{\ast }$, we have $0={\lambda }_1(\mathrm{\Delta }+\frac{1}{\nu }(\frac{\text{β}\left(x\right)}{m_{\ast }}-c))<{\lambda }_1(\mathrm{\Delta }+\frac{1}{\nu \delta \left(i\right)}(\frac{\text{β}\left(x\right)}{m_{\ast }}-c)),$$m^{\ast }>m_{\ast }$ follows from the monotonicity of ${\tilde{\lambda }}_1$ with respect to m.

In Theorem 3.7, as in the model with linear dispersal, if either the predator death rate or the half-capturing saturation rate is low, predators can invade a region. However, for the given half-capturing rate, if the number R is relatively large, then another threshold value exists other than that occurring in the model with linear diffusion (Theorem 3.3) related to predator dispersal such that the stability of the semitrivial steady state $({\theta }_{\mu },0)$ changes. Moreover, for the given PID, another threshold value other than that appearing in the model with linear diffusion (Theorem 3.3) relates to the half-capturing rate such that the stability of the semitrivial steady state $({\theta }_{\mu },0)$ changes. Figure 3 represents the contents of Theorem 3.7. It shows that if m is less than $\overline{\text{β}}/c$, then $({\theta }_{\mu },0)$ is unstable. However, curve ${\lambda }_1(m,\nu )=0$ for the linear diffusion case is shifted to the right of curve ${\lambda }_1(m,\nu )=0$ for PID, which means that the area for the instability of $({\theta }_{\mu },0)$ is enlarged in the $m\nu$-plane. Thus there is a range of half-capturing saturation rates m such that the coexistence of system (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ) exists even though (3(3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) ) cannot have coexistence. This implies that PID affects the evolution toward fitness for predator invasion in predator–prey interactions. Similarly, we can interpret the results with respect to ν (Theorem 3.7(iii)). The coexistence of system (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ) can be shown from Theorem 3.7 as Corollary 3.4.

Figure 3. Comparison of stability region of $({\theta }_{\mu },0)$ with respect to m and ν between linear dispersal and PID.

Corollary 3.8

Assume that $r-\alpha \left(x\right)>0$ for all $x\in \mathrm{\Omega }$ and that ρ is the number in Lemma 3.6(ii).

(i)	Suppose that $\rho >R$. Then, there is a coexistence of (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ) for any i>0.
(ii)	Suppose that $\rho <R$, and M and ${\nu }^{\ast }$ are given constants in Theorem 3.3. Then, for given m<M, there exists $p^{\ast }>{\nu }^{\ast }$ depending on $\mathrm{\Omega },m,\text{β},\delta$, and c such that if $\nu <p^{\ast }$, then there is a coexistence of (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ).
(iii)	Let $m_{\ast }$ be in Theorem 3.3 for a given ν. Then, for a given p, there exists $m^{\ast }\in (m_{\ast },M)$ depending on $\mathrm{\Omega },p,\text{β}$ and c such that if $m<m^{\ast }$, then there is a coexistence solution of (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ).

Note that the number $\rho \in (1,\frac{h}{\ell })$ is not specified in the preceding results.

3.3. Particular case of $p\left(i\right)$

Let us consider the particular form $p=p^{ϵ}$ defined in (2(2) $p^{ϵ}:=p^0\ast {\eta }^{ϵ},$(2) ) with $\ell =\nu$ where ν is the predator diffusion rate in system (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ). Note that $p^{ϵ}$ is smooth and satisfies (G) and (S). Thus, we can read Lemma 3.6 for this particular PID as follows:

Lemma 3.9

Let $\text{β},m,c$ and $p^{ϵ}$ be given.

(i)

Suppose $\overline{\text{β}}-mc>0$. Then, ${\int }_{\mathrm{\Omega }}\frac{1}{p^{ϵ}\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c)\mathrm{d}x>0$;

(ii)

Suppose $\overline{\text{β}}-mc<0$. If $\frac{h}{\ell }>R$, then ${\int }_{\mathrm{\Omega }}\frac{1}{p^{ϵ}\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c)\mathrm{d}x>0$, and if $\frac{h}{\ell }<R$, then ${\int }_{\mathrm{\Omega }}\frac{1}{p^{ϵ}\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c)\mathrm{d}x<0$ for a sufficiently small $ϵ>0$, where R is defined in (11(11) $R:=\frac{{\int }_{{\mathrm{\Omega }}_L}[c-\text{β}(x\left)/m\right]}{{\int }_{{\mathrm{\Omega }}_H}\left[\text{β}\right(x)/m-c]}$(11) ).

Furthermore, we can rewrite the stability result Theorem 3.7 of system (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ) for PID defined as (2(2) $p^{ϵ}:=p^0\ast {\eta }^{ϵ},$(2) ) (Figure 4).

Figure 4. Stability region of $({\theta }_{\mu },0)$ with respect to $\frac{h}{\ell }$ and ℓ.

Theorem 3.10

Let $({\theta }_{\mu },0)$ be a unique semitrivial steady-state solution of (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ), and ε be small such that Lemma 3.9(ii) holds.

(i)	Suppose that $\frac{h}{\ell }>R$. Then, $({\theta }_{\mu },0)$ is linearly unstable for any $\ell >0$.
(ii)	Suppose that $\frac{h}{\ell }<R$, and M and ${\nu }^{\ast }$ are given constants in Theorem 3.3(ii). Then, for a given m<M, there exists ${\ell }^{\ast }>{\nu }^{\ast }$ depending on $\mathrm{\Omega },m,\text{β},\ell ,h,ϵ$, and c such that if $\ell >{\ell }^{\ast }$, then $({\theta }_{\mu },0)$ is linearly asymptotically stable, and if $\ell <{\ell }^{\ast }$, then $({\theta }_{\mu },0)$ is linearly unstable.
(iii)	Let $m_{\ast }$ be in Theorem 3.3(iii) for a given ν. Then, for a given ℓ and h, there exists $m^{\ast }\in (m_{\ast },M)$ depending on $\mathrm{\Omega },\ell ,h,\text{β}$, and c such that if $m<m^{\ast }$, then $({\theta }_{\mu },0)$ is linearly unstable, and if $m>m^{\ast }$, then $({\theta }_{\mu },0)$ is linearly asymptotically stable.

The above theorem indicates that the dispersal sensitivity, $\frac{h}{l}$, is high, and the instability of the semitrivial state follows; thus, the predators can invade a region when rare, regardless of any PID. However, when the dispersal sensitivity, $\frac{h}{\ell }$, for the spatial environment is relatively small, another threshold value other than that occurring in the model with linear diffusion (Theorem 3.3) determines whether the predators can invade a region (see Figure 4). In this case, the predators obeying the PID can invade a region when rare, even if they cannot invade when dispersed linearly. Therefore there is a high probability of invading a region compared with the case with linear dispersal. The dispersal sensitivity, $\frac{h}{l}$, is crucial to determining the invasion of predators in a region. This particular case of $p\left(i\right)$ interprets the results more evidently than with the general shape of $p\left(i\right)$.

4. Coexistence states

4.1. Global/local bifurcation by linear dispersal

We state the instability (coexistence) region for $m,\nu$. Denote ${\lambda }_1(m,\nu )={\lambda }_1(\nu \mathrm{\Delta }+\frac{\text{β}\left(x\right)}{m}-c)$for a given β, c. From Lemma A.1, we have

Theorem 4.1

For $(m,\nu )\in (\overline{\text{β}}/c,M)\times {\mathbb{R}}^{\mathbb{+}}$,

(i)	${\lambda }_1(m,\nu )$ is decreasing with respect to m and ν.
(ii)	$\underset{\nu \to 0}{lim}{\lambda }_1(m,\nu )=\underset{x_{\in }\overline{\mathrm{\Omega }}}{max}(\frac{\text{β}}{m}-c)$.
(iii)	$\underset{\nu \to \mathrm{\infty }}{lim}{\lambda }_1(m,\nu )=\frac{1}{|\mathrm{\Omega }|}{\int }_{\mathrm{\Omega }}\left(\frac{\text{β}}{m}-c\right)$. Furthermore, $\frac{d\nu }{\mathrm{d}m}<0$.

Note that combining the above theorems, we obtain a continuous curve (Figure 2) : $\left\{\right(m,\nu )\in (\overline{\text{β}}/c,M)\times {\mathbb{R}}^{\mathbb{+}}:{\lambda }_1(m,\nu )=0\}.$Now, we obtain the global/local bifurcation and coexistence state of (3(3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) ): (12) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(12) First, we find the a priori bound of (12(12) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(12) ).

Lemma 4.2

Assume that $(u,v)$ is a positive solution of (12(12) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(12) ). Then, $||u|{|}_{\mathrm{\infty }}\le \underset{\overline{\mathrm{\Omega }}}{max}K\left(x\right)\mathrm{a}\mathrm{n}\mathrm{d}||v|{|}_{\mathrm{\infty }}\le \frac{\underset{\overline{\mathrm{\Omega }}}{max}\text{β}\left(x\right)\underset{\overline{\mathrm{\Omega }}}{max}K\left(x\right)}{c}.$

Assume that $\mu ,\nu$, $K\left(x\right)$, $\alpha \left(x\right)$, $\text{β}\left(x\right)$ and c are given. For all values of m, we have the branch of trivial solutions ${\mathrm{\Gamma }}^0=\left\{\right(m,0,0):m\in {\mathbb{R}}^{\mathbb{+}}\}$ and of semitrivial steady states ${\mathrm{\Gamma }}^1=\left\{\right(m,{\theta }_{\mu },0):m\in {\mathbb{R}}^{\mathbb{+}}\}$ to elliptic system (12(12) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(12) ). By applying the local bifurcation theorem of Crandall and Rabinowitz [14], we obtain a branch of positive solutions that bifurcates from ${\mathrm{\Gamma }}^1$. For p>n, we define Banach spaces X and Y as $X:=W_N^{2,p}\left(\mathrm{\Omega }\right)\times W_N^{2,p}\left(\mathrm{\Omega }\right),Y:=L^p\left(\mathrm{\Omega }\right)\times L^p\left(\mathrm{\Omega }\right),$where $W_N^{2,p}\left(\mathrm{\Omega }\right)=\{u\in W^{2,p}(\mathrm{\Omega }):\mathrm{\nabla }u\cdot \overrightarrow{n}=0\text{on}\mathrm{\partial }\mathrm{\Omega }\}$. Then, by the Sobolev embedding theorem, $X\subset C_0^1\left(\mathrm{\Omega }\right)\times C_0^1\left(\mathrm{\Omega }\right)$, where $C_0^1\left(\mathrm{\Omega }\right)$ is a set of continuously differentiable functions with compact support in Ω.

From Theorems A.2 and 3.3, there exists $m_{\ast }$ such that ${\lambda }_1(\nu \mathrm{\Delta }+\frac{\text{β}\left(x\right)}{m_{\ast }}-c)=0$. Thus, there is the positive function ${\psi }^{\ast }$ which is a solution of the equation: $\{\begin{array}{ll}{\displaystyle \nu \mathrm{\Delta }{\psi }^{\ast }+(\frac{\text{β}\left(x\right)}{m_{\ast }}-c){\psi }^{\ast }=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }{\psi }^{\ast }\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$where ${\int }_{\mathrm{\Omega }}{{\psi }^{\ast }}^2\mathrm{d}x=1$.

Choosing $m\in (0,\mathrm{\infty })$ as a bifurcation parameter, we obtain the local bifurcation property as follows.

Lemma 4.3

Assume that $\mu ,\nu ,K\left(x\right),\alpha \left(x\right),\text{β}\left(x\right)$, and c are given. Then, a branch of positive solutions to elliptic Equation (12(12) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(12) ) bifurcates from the semitrivial steady state solution curve ${\mathrm{\Gamma }}^1$ if and only if $m=m_{\ast }$. In other words, all positive solutions of the elliptic equation near $(m_{\ast },{\theta }_{\mu },0)\in \mathbb{R}\times C_0^1\left(\mathrm{\Omega }\right)\times C_0^1\left(\mathrm{\Omega }\right)$ can be represented by ${\mathrm{\Gamma }}^2=\left\{\right(m\left(s\right),{\theta }_{\mu }+s({\varphi }^{\ast }+u(s\left)\right),s({\psi }^{\ast }+v(s\left)\right)):0<s\le \delta \}$where ${\varphi }^{\ast }=(\mu \mathrm{\Delta }+r-\frac{2r{\theta }_{\mu }}{K\left(x\right)}{)}^{-1}(\frac{\alpha \left(x\right)}{m}{\psi }^{\ast })$. Here, $\left(m\right(s),u(s),v(s\left)\right)$ is a smooth function with respect to s and satisfies $\left(m\right(0),u(0),v(0\left)\right)=(m_{\ast },0,0)$ and ${\int }_{\mathrm{\Omega }}v\left(s\right){\psi }^{\ast }\mathrm{d}x=0$.

Before we investigate the global nature of the solution curve in Lemma 4.3, we give the bound of the values m for the existence of positive solutions to (12(12) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(12) ).

Theorem 4.4

If $(m,u,v)$ is a nontrivial, nonnegative solution of (12(12) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(12) ), then $m<m_{\ast }$.

The global bifurcation result for parameter m is obtained as follows:

Theorem 4.5

Assume that $\mu ,\nu ,K\left(x\right),\alpha \left(x\right),\text{β}\left(x\right)$, and c are given. If $m\in (0,m_{\ast })$, then (12(12) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(12) ) has at least one positive solution. The set of positive solutions of (12(12) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(12) ) with respect to bifurcation parameter m forms a continuum ${\mathrm{\Gamma }}^{\ast }\subset \mathbb{R}\times X$ satisfying the following properties:

(i)	${\mathrm{\Gamma }}^{\ast }$ bifurcates from ${\mathrm{\Gamma }}^1$ at $m=m_{\ast }$,
(ii)	$\{m:(m,u,v)\in {\mathrm{\Gamma }}^{\ast }\}=(0,m_{\ast })$.

4.2. Global/local bifurcation by PID

We investigate the global/local bifurcation and coexistence of the time-independent system of (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ): (13) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{V}{p\left(i\right)}\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }V\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$(13) where $V=p\left(i\right)v$ and $p\left(i\right)=\nu \delta \left(i\right)$. We regard $m\in (0,\mathrm{\infty })$ as a bifurcation parameter. First, we obtain the a priori bound of (13(13) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{V}{p\left(i\right)}\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }V\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$(13) ).

Lemma 4.6

Assume that $(u,V)$ is a positive solution of (13(13) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{V}{p\left(i\right)}\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }V\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$(13) ). Then $||u|{|}_{\mathrm{\infty }}\le \underset{\overline{\mathrm{\Omega }}}{max}K\left(x\right)\mathrm{a}\mathrm{n}\mathrm{d}||V|{|}_{\mathrm{\infty }}\le \frac{h}{c}\underset{\overline{\mathrm{\Omega }}}{max}\text{β}\left(x\right)\underset{\overline{\mathrm{\Omega }}}{max}K\left(x\right).$

Proof.

Let $x_0,x_1\in \overline{\mathrm{\Omega }}$ be points satisfying $u\left(x_0\right)=||u|{|}_{\mathrm{\infty }}$ and $V\left(x_1\right)=||V|{|}_{\mathrm{\infty }}$. By the maximum principle, $1-\frac{u\left(x_0\right)}{K\left(x_0\right)}-\alpha \left(x_0\right)\frac{v\left(x_0\right)}{mu\left(x_0\right)+v\left(x_0\right)}\ge 0$ and $\text{β}\left(x_1\right)\frac{u\left(x_1\right)}{mu\left(x_1\right)+v\left(x_1\right)}-c\ge 0$. From this, we have $\begin{array}{rl}1-\frac{u\left(x_0\right)}{K\left(x_0\right)}-\alpha \left(x_0\right)\frac{v\left(x_0\right)}{mu\left(x_0\right)+v\left(x_0\right)}& \ge 0\\ ⟺1-\frac{u\left(x_0\right)}{K\left(x_0\right)}& \ge 0.\end{array}$Thus, $u\left(x_0\right)\le K\left(x_0\right)\le \underset{\overline{\mathrm{\Omega }}}{max}K$. Similarly, $\begin{array}{rl}& \text{β}\left(x_1\right)\frac{u\left(x_1\right)}{mu\left(x_1\right)+v\left(x_1\right)}-c\ge 0\\ & ⟺c\left(mu\right(x_1)+v(x_1\left)\right)\le \text{β}\left(x_1\right)u\left(x_1\right)\\ & ⟺cv\left(x_1\right)\le \text{β}\left(x_1\right)u\left(x_1\right).\end{array}$Therefore, $V\left(x_1\right)=p\left(s\right(x_1\left)\right)v\left(x_1\right)\le \frac{h}{c}\text{β}\left(x_1\right)u\left(x_1\right)\le \frac{h}{c}\underset{\overline{\mathrm{\Omega }}}{max}\text{β}\underset{\overline{\mathrm{\Omega }}}{max}K$.

We assume that $\mu ,\nu$, $K\left(x\right)$, $\alpha \left(x\right)$, $\text{β}\left(x\right)$, and c are given. For all values of m, we have the branch of trivial solutions ${\mathrm{\Gamma }}^0=\left\{\right(m,0,0):m\in {\mathbb{R}}^{\mathbb{+}}\}$ and of semitrivial steady states ${\mathrm{\Gamma }}^1=\left\{\right(m,{\theta }_{\mu },0):m\in {\mathbb{R}}^{\mathbb{+}}\}$ of elliptic system (13(13) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{V}{p\left(i\right)}\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }V\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$(13) ). Using a similar argument as in Section 4.1, we can obtain the branch of positive solutions that bifurcates from ${\mathrm{\Gamma }}^1$.

From the proof of Theorem 3.7, ${\lambda }_1(\mathrm{\Delta }+\frac{1}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c))=0$ for $m=m^{\ast }$, and there is a positive function ${\psi }^{\ast }$ satisfying $\{\begin{array}{ll}{\displaystyle \mathrm{\Delta }{\psi }^{\ast }+\frac{1}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m^{\ast }}-c){\psi }^{\ast }=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }{\psi }^{\ast }\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$where ${\int }_{\mathrm{\Omega }}{{\psi }^{\ast }}^2\mathrm{d}x=1$.

We have the following local bifurcation property:

Lemma 4.7

Assume that $\mu ,p,K\left(x\right),\alpha \left(x\right),\text{β}\left(x\right)$, and c are given. Then, the branch of positive solutions of elliptic Equation (13(13) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{V}{p\left(i\right)}\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }V\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$(13) ) bifurcates from the semitrivial steady state solution curve ${\mathrm{\Gamma }}^1$ if and only if $m=m^{\ast }$. In other words, all positive solutions of the elliptic equation near $(m^{\ast },{\theta }_{\mu },0)\in {\mathbb{R}}^{+}\times C_0^1\left(\mathrm{\Omega }\right)\times C_0^1\left(\mathrm{\Omega }\right)$ can be expressed by ${\mathrm{\Gamma }}^2=\left\{\right(m\left(s\right),{\theta }_{\mu }+s({\varphi }^{\ast }+u(s\left)\right),s({\psi }^{\ast }+V(s\left)\right)):0<s\le \delta \}$where ${\varphi }^{\ast }=(\mu \mathrm{\Delta }+[r-\frac{2r{\theta }_{\mu }}{K\left(x\right)}\left]{)}^{-1}\right(\frac{\alpha \left(x\right)}{mp\left(i\right)}{\psi }^{\ast })$, $\left(m\right(s),u(s),V(s\left)\right)$ is a smooth function with respect to s and satisfies $\left(m\right(0),u(0),V(0\left)\right)=(m^{\ast },0,0)$ and ${\int }_{\mathrm{\Omega }}V\left(s\right){\psi }^{\ast }\mathrm{d}x=0$.

Proof.

Let $U=u-{\theta }_{\mu }$. Define the operator $F:{\mathbb{R}}^{+}\times C_0^1\left(\mathrm{\Omega }\right)\times C_0^1\left(\mathrm{\Omega }\right)\to L^p\times L^p$ by $F(m,U,V)=\left(\begin{array}{c}\mu \mathrm{\Delta }(U+{\theta }_{\mu })+(U+{\theta }_{\mu })[r-\frac{r(U+{\theta }_{\mu })}{K\left(x\right)}-\frac{\alpha \left(x\right)V}{mp\left(i\right)(U+{\theta }_{\mu })+V}]\\ \mathrm{\Delta }V+\frac{1}{p\left(i\right)}V[\frac{\text{β}\left(x\right)p\left(i\right)(U+{\theta }_{\mu })}{mp\left(i\right)(U+{\theta }_{\mu })+V}-c]\end{array}\right).$Since $(u,V)=({\theta }_{\mu },0)$ is the semitrivial steady state of elliptic Equation (13(13) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{V}{p\left(i\right)}\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }V\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$(13) ), we have $F(m,0,0)=0$. By simple calculations, we obtain $F_{(U,V)}(m,0,0)[\varphi ,\psi ]=\left(\begin{array}{c}\mu \mathrm{\Delta }\varphi +[r-\frac{2r{\theta }_{\mu }}{K\left(x\right)}]\varphi -\frac{\alpha \left(x\right)}{p\left(i\right)m}\psi \\ \mathrm{\Delta }\psi +\frac{1}{p\left(i\right)}[\frac{\text{β}\left(x\right)}{m}-c]\psi \end{array}\right).$Then, $\mathrm{Ker}{F}_{(U,V)}(m^{\ast },0,0)=\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\{{\varphi }^{\ast },{\psi }^{\ast }\}$, where ${\varphi }^{\ast }=(\mathrm{\Delta }+[r-\frac{2r{\theta }_{\mu }}{K\left(x\right)}\left]{)}^{-1}\right(\frac{\alpha \left(x\right)}{p\left(i\right)m}{\psi }^{\ast })<0$ in Ω. If $(\tilde{\varphi },\tilde{\psi })\in \text{Range}{F}_{(U,V)}(m^{\ast },0,0)$, then there exists $(\varphi ,\psi )\in X$ such that (14) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }\varphi +[r-\frac{2r{\theta }_{\mu }}{K\left(x\right)}]\varphi -\frac{\alpha \left(x\right)}{p\left(i\right)m}\psi =\tilde{\varphi }}& \\ {\displaystyle \nu \mathrm{\Delta }\psi +\frac{1}{p\left(i\right)}[\frac{\text{β}\left(x\right)}{m}-c]\psi =\tilde{\psi }}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }\varphi \cdot \overrightarrow{n}=\mathrm{\nabla }\psi \cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(14) By the Fredholm alternative theorem, the second equation of (14(14) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }\varphi +[r-\frac{2r{\theta }_{\mu }}{K\left(x\right)}]\varphi -\frac{\alpha \left(x\right)}{p\left(i\right)m}\psi =\tilde{\varphi }}& \\ {\displaystyle \nu \mathrm{\Delta }\psi +\frac{1}{p\left(i\right)}[\frac{\text{β}\left(x\right)}{m}-c]\psi =\tilde{\psi }}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }\varphi \cdot \overrightarrow{n}=\mathrm{\nabla }\psi \cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(14) ) is solvable if and only if ${\int }_{\mathrm{\Omega }}\tilde{\psi }{\psi }^{\ast }\mathrm{d}x=0$. For such a solution ψ, the first equation has a unique solution ϕ since $\mu \mathrm{\Delta }+r-\frac{2r{\theta }_{\mu }}{K\left(x\right)}$ is invertible. Thus, by the properties of the compact operator, $\mathrm{codim}(\mathrm{Range}{F}_{(U,V)}(m^{\ast },0,0\left)\right)=1.$Furthermore, we have $F_{(U,V),m}(m^{\ast },0,0)[{\varphi }^{\ast },{\psi }^{\ast }]=\left(\begin{array}{c}-\alpha \left(x\right)\frac{p\left(i\right)+i{p}^{\prime }\left(i\right)}{p^2\left(i\right)}{\psi }^{\ast }\\ -\frac{\text{β}\left(x\right)\left(p\right(i)+i(1-i\left)p^{\prime }\right(i\left)\right)}{m^2{p}^2\left(i\right)}{\psi }^{\ast }\end{array}\right).$Since ${\int }_{\mathrm{\Omega }}\frac{\text{β}\left(x\right)\left(p\right(i)+i(1-i\left)p^{\prime }\right(i\left)\right)}{m^2{p}^2\left(i\right)}({\psi }^{\ast }{)}^2>0$ by (S), we have $F_{(U,V),m}(m^{\ast },0,0)[\varphi ,\psi ]\notin \text{Range}{F}_{(U,V)}(m^{\ast },0,0).$Consequently, we can apply the local bifurcation theorem in [14] to F at $(m^{\ast },0,0)$. Note that there are no other bifurcation points except $m=m^{\ast }$ by virtue of the Krein–Rutman theorem.

Next, we give the bound of values, m, for the existence of positive solutions of (13(13) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{V}{p\left(i\right)}\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }V\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$(13) ).

Theorem 4.8

If $(m,u,V)$ is a nontrivial, nonnegative solution of (13(13) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{V}{p\left(i\right)}\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }V\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$(13) ), then $m<m^{\ast }$.

Proof.

Suppose that $m\ge m^{\ast }$. From (13(13) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{V}{p\left(i\right)}\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }V\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$(13) ) of the V-equation with $V=p\left(i\right)v$, $\begin{array}{rl}0& =\mathrm{\Delta }V+\frac{V}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c+\frac{\text{β}\left(x\right)u}{mu+v}-\frac{\text{β}\left(x\right)}{m})=\mathrm{\Delta }V\\ & +\frac{V}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c)-\frac{\text{β}\left(x\right)v^2}{m(mu+v)}.\end{array}$Thus, we obtain $\mathrm{\Delta }V+\frac{V}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c)=\frac{\text{β}\left(x\right)v^2}{m(mu+v)}.$Multiplying V by both sides and integrating on Ω, we have ${\int }_{\mathrm{\Omega }}-|\mathrm{\nabla }V{|}^2+V^2\frac{1}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c)\mathrm{d}x={\int }_{\mathrm{\Omega }}\frac{\text{β}\left(x\right)v^{2}V}{m(mu+v)}\mathrm{d}x>0.$On the other hand, from the result in Theorem 3.7, ${\lambda }_1(\mathrm{\Delta }+\frac{1}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c\left)\right)\le 0$ when $m\ge m^{\ast }$. Therefore, ${\int }_{\mathrm{\Omega }}-|\mathrm{\nabla }V{|}^2+V^2\frac{1}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c)\mathrm{d}x<{\lambda }_1(\mathrm{\Delta }+\frac{1}{p\left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c)){\int }_{\mathrm{\Omega }}{V}^2\mathrm{d}x\le 0,$which is a contradiction. Hence, $m<m^{\ast }$.

The global bifurcation structure is obtained as follows:

Theorem 4.9

Let $\mu ,p,K\left(x\right),\alpha \left(x\right),\text{β}\left(x\right)$, and c be given. Assume that $r-\alpha \left(x\right)>0$ for all $x\in \mathrm{\Omega }$. Then, (13(13) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{V}{p\left(i\right)}\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }V\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$(13) ) has at least one positive solution if and only if $m<m^{\ast }$. The set of positive solutions of (13(13) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{V}{p\left(i\right)}\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }V\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$(13) ) for bifurcation parameter m forms continuum ${\mathrm{\Gamma }}^{\ast }\subset \mathbb{R}\times X$ satisfying the following properties:

(i)	${\mathrm{\Gamma }}^{\ast }$ bifurcates from ${\mathrm{\Gamma }}^1$ at $m=m^{\ast }$,
(ii)	$\{m:(m,u,V)\in {\mathrm{\Gamma }}^{\ast }\}=(0,m^{\ast })$.

Proof.

For the local bifurcation branch ${\mathrm{\Gamma }}^2$ from Lemma 4.7, we denote ${\mathrm{\Gamma }}^{\ast }$ the maximum connected set of ${\mathrm{\Gamma }}^2\subset {\mathrm{\Gamma }}^{\ast }$. First, we introduce the positive cone: $P:=\left\{\right(u,V):u>0,V>0\mathrm{i}\mathrm{n}\mathrm{\Omega }\}.$We claim that ${\mathrm{\Gamma }}^{\ast }\setminus \left\{\right(m^{\ast },{\theta }_{\mu },0\left)\right\}\subset {\mathbb{R}}^{\mathbb{+}}\times P$. Suppose that $(\hat{m},\hat{u},\hat{V})\in {\mathrm{\Gamma }}^{\ast }$ satisfies $(\hat{u},\hat{V})\in \mathrm{\partial }P$. Then, it follows that $\hat{u}\ge 0,\hat{V}\ge 0$ for $x\in \mathrm{\Omega }$, and so $\hat{u}\left(x_0\right)\hat{V}\left(x_0\right)=0\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{s}\mathrm{o}\mathrm{m}\mathrm{e}{x}_0\in \mathrm{\Omega }.$By the strong maximum principle, $\hat{u}\equiv 0$ or $\hat{V}\equiv 0$.

We suppose that $\hat{u}\equiv \hat{V}\equiv 0$. Then, $(\hat{m},\hat{u},\hat{V})$ lies on the trivial branch of solutions ${\mathrm{\Gamma }}^0$. Let $\left\{\right(m_n,u_n,V_n\left)\right\}$ be the sequence of positive solutions of (13(13) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{V}{p\left(i\right)}\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }V\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$(13) ) with $m_n$ in ${\mathrm{\Gamma }}^{\ast }\cap \left(\right(0,m^{\ast })\times P)$. By Harnack's inequality, there exists $C_1>0$ such that $C_1\underset{\overline{\mathrm{\Omega }}}{max}K\le u_n\le \underset{\overline{\mathrm{\Omega }}}{max}K$. By dividing the first equation of (13(13) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{V}{p\left(i\right)}\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }V\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$(13) ) by $||u_n|{|}_{\mathrm{\infty }}$ and integrating over Ω, we have (15) $\begin{array}{rl}0& ={\int }_{\mathrm{\Omega }}\frac{r{u}_n}{||u_n|{|}_{\mathrm{\infty }}}(1-\frac{u_n}{K\left(x\right)})-\alpha \left(x\right)\frac{u_n{v}_n}{||u_n|{|}_{\mathrm{\infty }}(m_n{u}_n+v_n)}\mathrm{d}x\end{array}$(15) (16) $\begin{array}{rl}& \ge {\int }_{\mathrm{\Omega }}\frac{r{u}_n}{||u_n|{|}_{\mathrm{\infty }}}(1-\frac{u_n}{K\left(x\right)})-\alpha \left(x\right)\frac{u_n{v}_n}{C_1\underset{\overline{\mathrm{\Omega }}}{max}K(m_n{u}_n+v_n)}\mathrm{d}x\end{array}$(16) where $v_n=\frac{V_n}{p\left(i\right)}$. Note that $\frac{u_n}{||u_n|{|}_{\mathrm{\infty }}}>0$ for sufficiently large n. Therefore, for large n, $\frac{r{u}_n}{||u_n|{|}_{\mathrm{\infty }}}(1-\frac{u_n}{K\left(x\right)})-\alpha \left(x\right)\frac{u_n{v}_n}{C_1\underset{\overline{\mathrm{\Omega }}}{max}K(m_n{u}_n+v_n)}>0,$which the second integration (15(15) $\begin{array}{rl}0& ={\int }_{\mathrm{\Omega }}\frac{r{u}_n}{||u_n|{|}_{\mathrm{\infty }}}(1-\frac{u_n}{K\left(x\right)})-\alpha \left(x\right)\frac{u_n{v}_n}{||u_n|{|}_{\mathrm{\infty }}(m_n{u}_n+v_n)}\mathrm{d}x\end{array}$(15) ) is positive, and so it is a contradiction. Thus, there is no sequence ${\mathrm{\Gamma }}^{\ast }\cap \left(\right(0,m^{\ast })\times P)$ converging to $(\hat{m},\hat{u},\hat{V})$. Therefore, either $\hat{u}$ or $\hat{V}$ should be nonzero. Suppose that $\hat{u}$ is zero. Then, $\mathrm{\Delta }V-\frac{cV}{p\left(i\right)}=0$. However, since this equation cannot have a nonzero solution, $\hat{u}$ is not identically 0. Suppose that $\hat{V}$ is zero. Since $m^{\ast }$ is the only bifurcation point, $\hat{m}=m^{\ast }$, which is impossible. Therefore, ${\mathrm{\Gamma }}^{\ast }\setminus \left\{\right(m^{\ast },{\theta }_{\mu },0\left)\right\}\subset {\mathbb{R}}^{\mathbb{+}}\times P$.

Next, by global bifurcation theory [38,49], one of the following nonexcluding properties holds:

1. ${\mathrm{\Gamma }}^{\ast }$ is unbounded in ${\mathbb{R}}^{+}\times X$;

2. ${\mathrm{\Gamma }}^{\ast }$ meets the semitrivial steady state set $\left\{\right(m,{\theta }_{\mu },0):m\in {\mathbb{R}}^{\mathbb{+}}\}$ at another bifurcation point, except $(m^{\ast },{\theta }_{\mu },0)$;

3. ${\mathrm{\Gamma }}^{\ast }$ connects to another component on the boundary of ${\mathbb{R}}^{+}\times P$.

From the a priori bound of solutions of (13(13) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{V}{p\left(i\right)}\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }V\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$(13) ) (Lemmas 4.6) and 4.8, (i) cannot hold. Additionally, (ii) cannot occur since the positive solutions bifurcate from the semitrivial steady state solution curve ${\mathrm{\Gamma }}^1$ only at $m=m^{\ast }$ by Lemma 4.7. Hence, alternative (iii) must hold. Then, the positive solution $(m,u,V)$ of (13(13) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{V}{p\left(i\right)}\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }V\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$(13) ) on the boundary of ${\mathbb{R}}^{+}\times P$ has the form $(m,{\theta }_{\mu },0)$, $(m,0,0)$ and $(0,u,V)$. Since $m^{\ast }$ is a unique bifurcation point and $(0,0)$ is unstable for any m under the assumption $r-\alpha \left(x\right)>0$ for all $x\in \mathrm{\Omega }$, $(0,u,V)$ is the only possible form of the positive solution. Therefore, $\{m:(m,u,V)\in {\mathrm{\Gamma }}^{\ast }\}=(0,m^{\ast })$.

Figure 5 is the bifurcation diagram representing the curve of the density of predators from the semitrivial state $({\theta }_{\mu },0)$ for a system with random dispersal and PID of predators, as half-capturing saturation constant m is a bifurcation parameter. The bifurcation curves begin at threshold values $m_{\ast }$ and $m^{\ast }$ with $m_{\ast }<m^{\ast }$, which occur in Theorems 3.3(iii) and 3.7(iv), respectively. Thus, there is a range of m such that $m_{\ast }<m<m^{\ast }$, which indicates that if the predators have motility obeying PID, they can coexist with the prey even though there is no coexistence when the predators diffuse linearly.

Figure 5. Bifurcation diagram with parameter m.

Furthermore, we can take ν as a bifurcation parameter for systems (1(1) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\mathrm{\Delta }\left(p\right(i\left)v\right)+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }\left(p\right(i\left)v\right)\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega },\end{array}$(1) ) and (3(3) $\{\begin{array}{ll}{\displaystyle u_t=\mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}}& \\ {\displaystyle v_t=\nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv}& \text{in}\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }\times (0,\mathrm{\infty }),\\ u(x,0)=u_0\left(x\right)\ge 0,v(x,0)=v_0\left(x\right)\ge 0& \text{in}\mathrm{\Omega }.\end{array}$(3) ).

The bifurcation diagram in Figure 6 is provided for a bifurcation parameter ν. The curves start at two distinct threshold values, ${\nu }^{\ast }$ and $p^{\ast }$, from Theorems 3.3(ii) and 3.7(iii), respectively. Thus, as in the case where m serves as a bifurcation parameter, there are values ν such that ${\nu }^{\ast }<\nu <p^{\ast }$. This result again shows that when the predators move following PID, the coexistence state exists even though the predators cannot survive with the prey when the former diffuses linearly. Since the results are similar to the results obtained when we take m as a parameter, we present them in Appendix 2.

Figure 6. Bifurcation diagram with parameter ν.

5. Conclusion

In ecology, the motility of species in a certain region depends on the circumstances of the local habitat. In this research, we studied the evolution of dispersal toward fitness for the invasion of predators, called PID as a dispersal strategy by examining the effect of PID on the stability of semitrivial steady states. First, we considered a predator–prey model with linear dispersal. If the predator death rate or the half-capturing saturation rate is low for any dispersal rate of the predator, the predator can invade a region when rare. Moreover, when the half-capturing saturation constant is in a particular range, the threshold value relates to the predator dispersal rate, which determines the stability of the prey where the predator is absent. Furthermore, for a given predator dispersal rate, the threshold value is related to the half-capturing saturation constant that changes the prey's stability where the predator is absent. In addition, the instability of the semitrivial steady state implies the coexistence of the time-dependent system from the positive density of the prey because there is only one semitrivial steady state, $({\theta }_{\mu },0)$ (Figure 2). Whereas, when the predator follows PID, if the predator death rate or the predation degree is sufficiently low for any PID of the predator, it can invade a region when rare. In addition, given the half-capturing rate, when the predation degree is relatively high, another threshold value other than that occurring in the model with linear dispersal relates to the dispersal rate of the predator such that the stability of the semitrivial steady state $({\theta }_{\mu },0)$ changes. Furthermore, for a given PID, another threshold value other than that occurring in the model with linear dispersal relates to the half-capturing rate such that the stability of the semitrivial steady state $({\theta }_{\mu },0)$ is changeable (Figure 3). Thus, the predator has a higher probability of invading a region than in the case with linear dispersal (Figure 3).

In particular, when PID is defined as the shape of a step function with bounds h and ℓ, we obtain more applicable conditions for stability. More precisely, when the dispersal sensitivity $\frac{h}{\ell }$ is high, the predator can invade a region, regardless of any PID (Figure 4). That is, the high sensitivity of predator dispersal always makes predators survive even in a low predation domain (Figure 7). In contrast, if the dispersal sensitivity $\frac{h}{\ell }$ is relatively low, another threshold value other than that occurring in the model with linear dispersal determines the invasion of the predator in a region (Figure 4). In this case, the predator with PID can invade when rare in a region if it cannot invade when dispersed linearly (Figures 8 and 9). Thus, dispersal sensitivity is crucial in determining the invasion of predators in a region for the PID of step-function PID. We can conclude that the evolution of dispersal toward fitness for predator invasion to a given region for predators progresses effectively when their movement obeys PID rather than linear dispersal.

Figure 7. Effect of sensitivity $\frac{h}{\ell }$ when the predator follows PID : m = 1.15, $\ell =0.52$. (a) Stabilities of $({\theta }_{\mu },0)$ when $\frac{h}{\ell }=1.4$. (b) Instability of $({\theta }_{\mu },0)$ when $\frac{h}{\ell }=5.0$.

Figure 8. Effect of PID : m = 0.78. (a) Stabilities of $({\theta }_{\mu },0)$ when the predator diffusion rate $\nu =0.01$. (b) Instability of $({\theta }_{\mu },0)$ when the predator follows PID with $\ell =0.01$ and $\frac{h}{\ell }=1.4$.

Figure 9. Effect of PID on predator diffusion : m = 1.15. (a) Stabilities of $({\theta }_{\mu },0)$ when the predator diffusion rate $\nu =0.2$. (b) Instability of $({\theta }_{\mu },0)$ when the predator follows PID with $\ell =0.2$, $\frac{h}{\ell }=1.4$.

Finally, we obtained the coexistence results of the time-independent system obtained using a bifurcation branch, starting a unique semitrivial steady state of models with a linear diffusion and PID, also supporting that the PID of predators can increase fitness for invasion of the predators when rare.

In this paper, we adopted the ratio-dependent functional response to describe the predator–prey interaction. If the prey-dependent functional response, such as Holling type II, is considered, we can see a similar effect of PID on predator survival, but it would be more complicated, and some different results from those of this paper, depending on the prey's conditions, will also be obtained. Finally, we note that the uniqueness and stability of the coexistence state is an important issue in the mathematical ecology model, especially for predator–prey models in a spatially heterogeneous environment. We leave this challenging issue open and hope to be solved in future works.

Acknowledgments

The authors greatly appreciate the valuable comments and suggestions of anonymous reviewers of this paper and the editor, which significantly improved the original version.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT & Future Planning (NRF-2020R1F1A1A01054935, and NRF-2022R1F1A1063068) and Ministry of Education (NRF-2021R1A6A3A01086879).

Appendices

Appendix 1.

Eigenvalue problem

We give some properties of the following eigenvalue problem that will be mainly used to obtain the stability results: (A1) $\{\begin{array}{ll}\lambda \varphi =\nu \mathrm{\Delta }\varphi +g\left(x\right)\varphi & \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }\varphi \cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(A1) We call the largest eigenvalue of (A1(A1) $\{\begin{array}{ll}\lambda \varphi =\nu \mathrm{\Delta }\varphi +g\left(x\right)\varphi & \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }\varphi \cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(A1) ) the principal eigenvalue, denoted by ${\lambda }_1(\nu \mathrm{\Delta }+g(x\left)\right)$, and defined as ${\lambda }_1(\nu \mathrm{\Delta }+g(x\left)\right)=\underset{\varphi \in H^1\left(\mathrm{\Omega }\right)\setminus \left\{0\right\}}{sup}\frac{{\int }_{\mathrm{\Omega }}-\nu |\mathrm{\nabla }\varphi {|}^2+g\left(x\right){\varphi }^2\mathrm{d}x}{{\int }_{\mathrm{\Omega }}{\varphi }^2\mathrm{d}x}.$Then, there exists the positive principal eigenfunction, ${\varphi }_1$, corresponding to ${\lambda }_1(\nu \mathrm{\Delta }+g(x\left)\right)$.

We introduce the following lemmas about the well-known property of the principal eigenvalue and some conditions for determining the sign of the principal eigenvalue ${\lambda }_1(\nu \mathrm{\Delta }+g(x\left)\right)$. We omit the proofs here; see [9,19].

Lemma A.1

(i)	${\lambda }_1(\nu \mathrm{\Delta }+g(x\left)\right)$ is monotonically decreasing with respect to ν.
(ii)	${\lambda }_1(\nu \mathrm{\Delta }+\tilde{g}(x\left)\right)\ge {\lambda }_1(\nu \mathrm{\Delta }+g(x\left)\right)$ for $\tilde{g}\left(x\right)\ge g\left(x\right)$.
(iii)	$\underset{\nu \to 0}{lim}{\lambda }_1(\nu \mathrm{\Delta }+g(x\left)\right)=\underset{x\in \overline{\mathrm{\Omega }}}{max}g\left(x\right)$.
(iv)	$\underset{\nu \to \mathrm{\infty }}{lim}{\lambda }_1(\nu \mathrm{\Delta }+g(x\left)\right)=\frac{1}{|\mathrm{\Omega }|}{\int }_{\mathrm{\Omega }}g\left(x\right)\mathrm{d}x$.

Theorem A.2

(i)	If $g\left(x\right)<0$, then ${\lambda }_1(\nu \mathrm{\Delta }+g(x\left)\right)<0$.
(ii)	If ${\int }_{\mathrm{\Omega }}g\left(x\right)\mathrm{d}x>0$, then ${\lambda }_1(\nu \mathrm{\Delta }+g(x\left)\right)>0$.
(iii)	If ${\int }_{\mathrm{\Omega }}g\left(x\right)\mathrm{d}x<0$ and $g\left(x\right)$ is positive somewhere in Ω, then there exists ${\nu }^{\ast }>0$ such that ${\lambda }_1(\nu \mathrm{\Delta }+g(x\left)\right)\ge 0$ if and only if $\nu \le {\nu }^{\ast }$. The equality holds when $\nu ={\nu }^{\ast }$.

Appendix 2.

Global bifurcation by parameter ν

A.1. Model with linear diffusion

We choose $\nu \in (0,\mathrm{\infty })$ as a bifurcation parameter. Then, we obtain a bifurcation curve from the semitrivial steady state of elliptic Equation (12(12) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(12) ) by considering ν as a bifurcation parameter. From Theorem 3.3, there are no bifurcation phenomena for $m\ge M$ and $m\le \overline{\text{β}}/c$. Thus, we assume that $m\in (\overline{\text{β}}/c,M)$. Then, we have the corresponding results to the above case in which m acts as a bifurcation parameter when m is replaced by ν.

From Theorems A.2 and 3.3, ${\lambda }_1({\nu }^{\ast }\mathrm{\Delta }+\frac{\text{β}\left(x\right)}{m}-c)=0$ for $\nu ={\nu }^{\ast }$: ${\psi }_{\ast }>0$ is the solution of the equation $\{\begin{array}{ll}{\displaystyle \mathrm{\Delta }{\psi }_{\ast }+\frac{1}{{\nu }^{\ast }}(\frac{\text{β}\left(x\right)}{m}-c){\psi }_{\ast }=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }{\psi }_{\ast }\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$where ${\int }_{\mathrm{\Omega }}{{\psi }_{\ast }}^2\mathrm{d}x=1$.

The local bifurcation property holds as follows.

Lemma A.3

Assume that $\mu ,K\left(x\right),\alpha \left(x\right),\text{β}\left(x\right)$, and c are given, and $m\in (\overline{\text{β}}/c,M)$. Then, a positive solution of elliptic Equation (12(12) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(12) ) bifurcates from the semitrivial steady state solution curve ${\mathrm{\Gamma }}_1$ if and only if $\nu ={\nu }^{\ast }$. In other words, all positive solutions of (12(12) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(12) ) near $({\nu }^{\ast },{\theta }_{\mu },0)\in {\mathbb{R}}^{+}\times C_0^1\left(\mathrm{\Omega }\right)\times C_0^1\left(\mathrm{\Omega }\right)$ can be expressed by ${\mathrm{\Gamma }}_2=\left\{\right(\nu \left(s\right),{\theta }_{\mu }+s({\varphi }_{\ast }+u(s\left)\right),s({\psi }_{\ast }+v(s\left)\right)):0<s\le \delta \},$where ${\varphi }_{\ast }=(\mu \mathrm{\Delta }+r-\frac{2r{\theta }_{\mu }}{K\left(x\right)}{)}^{-1}(\frac{\alpha \left(x\right)}{m}{\psi }_{\ast })$. Here, $\left(\nu \right(s),u(s),v(s\left)\right)$ is a smooth function with respect to s and satisfies $\left(\nu \right(0),u(0),v(0\left)\right)=({\nu }^{\ast },0,0)$ and ${\int }_{\mathrm{\Omega }}v\left(s\right){\psi }_{\ast }\mathrm{d}x=0$.

We state the global bifurcation result for system (12(12) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(12) ) with the bound of $\nu <{\nu }^{\ast }$ such that (12(12) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(12) ) has a coexistence steady state.

Theorem A.4

Let $\mu ,K\left(x\right),\alpha \left(x\right),\text{β}\left(x\right),c$, and $m\in (\overline{\text{β}}/c,M)$ be given. Assume that $r-\alpha \left(x\right)>0$ for all $x\in \mathrm{\Omega }$. Then (12(12) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(12) ) has at least one positive solution if and only if $\nu <{\nu }^{\ast }$. The set of positive solutions of (12(12) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \nu \mathrm{\Delta }v+\text{β}\left(x\right)\frac{uv}{mu+v}-cv=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(12) ) with bifurcation parameter ν forms a continuum ${\mathrm{\Gamma }}_{\ast }\subset \mathbb{R}\times X$ satisfying the following:

(i)	${\mathrm{\Gamma }}_{\ast }$ bifurcates from ${\mathrm{\Gamma }}_1$ at $\nu ={\nu }^{\ast }$,
(ii)	$\{\nu :(\nu ,u,v)\in {\mathrm{\Gamma }}_{\ast }\}=(0,{\nu }^{\ast })$.

A.2. Model with PID

Next, we investigate the bifurcation curve on the semitrivial steady state set of elliptic Equation (13(13) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{V}{p\left(i\right)}\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }V\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$(13) ) by considering ν as a parameter, where δ is a shape of p such that $p=\nu \delta$. From Theorem 3.7(iv), there are no bifurcation phenomena for $m\ge M$ or $m<m_{\ast }$. Thus, we assume that $m\in (m_{\ast },M)$. Consider the equivalent form of (13(13) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{V}{p\left(i\right)}\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }V\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$(13) ) (A2) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{1}{\nu \delta }V\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(A2) Regarding ν as a bifurcation parameter, we have the branch of trivial solutions, ${\mathrm{\Gamma }}_0=\left\{\right(\nu ,0,0):\nu \in {\mathbb{R}}^{+}\}$ and the branch of semitrivial steady states, ${\mathrm{\Gamma }}_1=\left\{\right(\nu ,{\theta }_{\mu },0):\nu \in {\mathbb{R}}^{+}\}$, of elliptic system (A2(A2) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{1}{\nu \delta }V\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(A2) ). Then, as before, we obtain a branch of positive solutions that bifurcates from ${\mathrm{\Gamma }}_1$. For p>n, we define Banach spaces X and Y as $X:=W_N^{2,p}\left(\mathrm{\Omega }\right)\times W_N^{2,p}\left(\mathrm{\Omega }\right),Y:=L^p\left(\mathrm{\Omega }\right)\times L^p\left(\mathrm{\Omega }\right),$where $W_N^{2,p}\left(\mathrm{\Omega }\right)=\{u\in W^{2,p}(\mathrm{\Omega }):\mathrm{\nabla }u\cdot \overrightarrow{n}=0\text{on}\mathrm{\partial }\mathrm{\Omega }\}$. By the Sobolev embedding theorem, $X\subset C_0^1\left(\mathrm{\Omega }\right)\times C_0^1\left(\mathrm{\Omega }\right)$.

From the proof of Theorem 3.7, ${\lambda }_1(\mathrm{\Delta }+\frac{1}{p^{\ast }\delta \left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c\left)\right)=0$. There is a positive function ${\psi }_{\ast }$ is a solution of the equation $\{\begin{array}{ll}{\displaystyle \mathrm{\Delta }{\psi }_{\ast }+\frac{1}{p^{\ast }\delta \left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c){\psi }_{\ast }=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }{\psi }_{\ast }\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$where ${\int }_{\mathrm{\Omega }}{{\psi }_{\ast }}^2\mathrm{d}x=1$.

The following local bifurcation property holds.

Lemma A.5

Assume that μ, K, α, β, and c are given, and $m\in (m_{\ast },M)$. The branch of positive solutions of elliptic Equation (A2(A2) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{1}{\nu \delta }V\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(A2) ) bifurcates from the semitrivial steady state solution curve ${\mathrm{\Gamma }}_1$ if and only if $\nu =p^{\ast }$. In other words, all positive solutions of (A2(A2) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{1}{\nu \delta }V\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(A2) ) near $(p^{\ast },{\theta }_{\mu },0)\in \mathbb{R}\times C_0^1\left(\mathrm{\Omega }\right)\times C_0^1\left(\mathrm{\Omega }\right)$ can be represented by ${\mathrm{\Gamma }}_2=\left\{\right(\nu \left(s\right),{\theta }_{\mu }+s({\varphi }_{\ast }+u(s\left)\right),s({\psi }_{\ast }+V(s\left)\right)):0<s\le \delta \}$where ${\varphi }_{\ast }=(\mu \mathrm{\Delta }+[r-\frac{2r{\theta }_{\mu }}{K\left(x\right)}\left]{)}^{-1}\right(\frac{\alpha \left(x\right)}{mp\left(i\right)}{\psi }_{\ast })$. Here, $\left(\nu \right(s),u(s),V(s\left)\right)$ is a smooth function with respect to s and satisfies $\left(\nu \right(0),u(0),V(0\left)\right)=(p^{\ast },0,0)$ and ${\int }_{\mathrm{\Omega }}V\left(s\right){\psi }_{\ast }\mathrm{d}x=0$.

There is a bound of values ν for the existence of positive solutions of (A2(A2) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{1}{\nu \delta }V\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(A2) ).

Theorem A.6

If $(\nu ,u,V)$ is a nontrivial, nonnegative solution of (A2(A2) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{1}{\nu \delta }V\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(A2) ), then $\nu <p^{\ast }$.

Proof.

Suppose that $\nu \ge p^{\ast }$. Then, from the second Equation (A2(A2) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{1}{\nu \delta }V\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }v\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(A2) ), we have $\mathrm{\Delta }V+\frac{1}{\nu \delta \left(i\right)}V(\frac{\text{β}\left(x\right)}{m}-c)=\frac{\text{β}\left(x\right)v^2}{m(mu+v)}.$Multiply v by both sides and integration, then ${\int }_{\mathrm{\Omega }}-|\mathrm{\nabla }V{|}^2+\frac{1}{\nu \delta \left(i\right)}{V}^2(\frac{\text{β}\left(x\right)}{m}-c)\mathrm{d}x={\int }_{\mathrm{\Omega }}\frac{\text{β}\left(x\right)v^{2}V}{m(mu+v)}\mathrm{d}x>0.$From the result in Theorem 3.7, ${\lambda }_1(\mathrm{\Delta }+\frac{1}{\nu \delta \left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c\left)\right)\le 0$ when $\nu \ge p^{\ast }$. Therefore, ${\int }_{\mathrm{\Omega }}-|\mathrm{\nabla }V{|}^2+\frac{1}{\nu \delta \left(i\right)}{V}^2(\frac{\text{β}\left(x\right)}{m}-c)\mathrm{d}x<{\lambda }_1(\mathrm{\Delta }+\frac{1}{\nu \delta \left(i\right)}(\frac{\text{β}\left(x\right)}{m}-c)){\int }_{\mathrm{\Omega }}{V}^2\mathrm{d}x\le 0,$and so we have a contradiction. Hence, $\nu <p^{\ast }$.

Finally, we obtain the following global bifurcation result of (13(13) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{V}{p\left(i\right)}\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }V\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$(13) ) for parameter ν as a similar argument of Theorem 4.9.

Theorem A.7

Let $\mu ,K,\alpha ,\text{β},c$, and $m\in (m_{\ast },M)$ be given. Assume that $r-\alpha \left(x\right)>0$ for all $x\in \mathrm{\Omega }$. Then (13(13) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{V}{p\left(i\right)}\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }V\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$(13) ) has at least one positive solution if and only if $\nu \in (0,p^{\ast })$. The set of positive solutions of (13(13) $\{\begin{array}{ll}{\displaystyle \mu \mathrm{\Delta }u+ru(1-\frac{u}{K\left(x\right)})-\alpha \left(x\right)\frac{uv}{mu+v}=0}& \\ {\displaystyle \mathrm{\Delta }V+\frac{V}{p\left(i\right)}\left(\text{β}\right(x)\frac{u}{mu+v}-c)=0}& \text{in}\mathrm{\Omega },\\ \mathrm{\nabla }u\cdot \overrightarrow{n}=\mathrm{\nabla }V\cdot \overrightarrow{n}=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$(13) ) with bifurcation parameter ν forms continuum ${\mathrm{\Gamma }}_{\ast }\subset \mathbb{R}\times X$ satisfying the following properties:

(i)	${\mathrm{\Gamma }}_{\ast }$ bifurcates from ${\mathrm{\Gamma }}_1$ at $\nu =p^{\ast }$,
(ii)	$\{\nu :(\nu ,u,v)\in {\mathrm{\Gamma }}_{\ast }\}=(0,p^{\ast })$.