Lanzhong Zhang
College of Mathematics and Statistics, Northwest Normal University, Lanzhou, China.
DOI: 10.4236/jamp.2026.145092   PDF    HTML   XML   13 Downloads   101 Views   Citations

Abstract

This paper studies a class of predator-prey models characterized by predator attraction and fear effects. By employing the upper and lower solution method and the regularity theory of elliptic equations, we prove the uniform boundedness of positive solutions. Additionally, using the prey growth rate as a bifurcation parameter, we apply the Crandall-Rabinowitz local bifurcation theory to investigate the local bifurcation structure of positive solutions branching off from two semi-trivial solutions and analyze the local stability of the bifurcated solutions.

Keywords

Holling-Tanner Predation Model, Predator-Taxis, Fear Effect

Share and Cite:

1. Introduction

In ecosystems, the coexistence and evolution of species are influenced by multiple factors, including population migration, interspecific interactions, and environmental heterogeneity. In recent years, mathematical studies of predator-prey models have garnered increasing attention, where diffusion, chemotaxis, and fear effects are considered key mechanisms affecting the spatial distribution and dynamic behavior of populations.

In modeling predator-prey interactions, chemotaxis describes the directed movement response of organisms to external stimuli (such as food or predation risk). Kareiva and Odell [1] were among the first to propose a prey-taxis model to describe the aggregation behavior of predators towards areas of higher prey density. Subsequently, numerous studies have further explored predator-prey models incorporating prey-taxis [2]-[6], with reference [6] investigating prey-taxis in the Holling-Tanner predator-prey model. This study revealed the significant impact of the taxis coefficient as a bifurcation parameter on the existence of non-normal steady-state solutions and elucidated the crucial role of prey-taxis in the formation of spatial patterns.

In addition to the active search of predators for prey, prey also employ various anti-predation strategies to evade predation risk. Most species possess the ability to perceive predation risk and adjust their behavior accordingly [7]. Among these, predator-taxis is an important anti-predation behavior, where prey migrate towards areas of lower predator density in response to perceived predation risk. Previous research has proposed and analyzed predator-prey models that incorporate predator-taxis [8]-[12].

Moreover, fear effects are also significant factors influencing population dynamics. Wang et al. [13] were the first to introduce fear effects into a predator-prey model with a Holling-II functional response, finding that higher levels of fear could stabilize predator-prey dynamics by eliminating periodic solutions, while lower levels of fear might lead to multiple limit cycles through subcritical Hopf bifurcation, resulting in bistability. Subsequent studies [14]-[17] have further explored the impact of fear effects in various predator-prey models.

Inspired by the aforementioned studies, this paper considers a more complex predator-prey model. Unlike the research in reference [6], which focuses on predator foraging behavior (i.e., prey-taxis), this paper introduces both predator-taxis and fear effects into the Holling-Tanner predator-prey model. Additionally, the model incorporates environmental factors of spatial heterogeneity, where the intra-specific competition coefficient of predators, $m\left(x\right)$ , depends on spatial location, thereby reflecting the complexity of real ecological environments. Building upon these motivations, we consider the full time-dependent model, which is a parabolic system describing the evolution of prey density $u\left(x,t\right)$ and predator density $v\left(x,t\right)$

$\{\begin{array}{ll}{u}_t=d_1\Delta u+\alpha \nabla \cdot \left(u\nabla v\right)+\frac{r_{0}u}{1+k_0\alpha v}-du-a{u}^2-\frac{puv}{1+qu},\hfill & x\in \Omega ,t>0,\hfill \\ v_t=d_2\Delta v+sv-\frac{m\left(x\right)v^2}{r+u}+\frac{cpuv}{1+qu},\hfill & x\in \Omega ,t>0,\hfill \\ \frac{\partial u}{\partial \mu }=0,\frac{\partial v}{\partial \mu }=0,\hfill & x\in \partial \Omega ,t>0,\hfill \\ u\left(x,0\right)=u_0\left(x\right)>0,v\left(x,0\right)=v_0\left(x\right)>0,\hfill & x\in \Omega .\hfill \end{array}$ (1)

This paper is devoted to the analysis of its steady-state problem, which yields the stationary solutions that dictate the long-term spatial patterns. Therefore, we investigate the following elliptic system

$\{\begin{array}{ll}-d_1\Delta u-\alpha \nabla \cdot \left(u\nabla v\right)=\frac{r_{0}u}{1+k_0\alpha v}-du-a{u}^2-\frac{puv}{1+qu},\hfill & x\in \Omega ,\hfill \\ -d_2\Delta v=sv-\frac{m\left(x\right)v^2}{r+u}+\frac{cpuv}{1+qu},\hfill & x\in \Omega ,\hfill \\ \frac{\partial u}{\partial \mu }=0,\frac{\partial v}{\partial \mu }=0,\hfill & x\in \partial \Omega .\hfill \end{array}$ (2)

Let Ω be a bounded domain in ${ℝ}^n$ with a smooth boundary $\partial \Omega$ , and let $\mu$ denote the unit outward normal vector on the boundary. The variables $u$ and $v$ represent the population densities of the prey and predator, respectively. The positive constants $d_1$ and $d_2$ are their random diffusion coefficients. The advection term $-\alpha \nabla \cdot \left(u\nabla v\right)$ models predator-taxis, describing the directed movement of the prey away from regions of high predator density. The coefficient $\alpha >0$ represents the prey’s sensitivity to predation risk. The same risk perception triggers a fear effect, which reduces the prey’s reproduction rate; this is captured in the term $\frac{r_{0}u}{1+k_0\alpha v}$ , where $k_0>0$ scales the impact of fear on the birth rate $r_0$ . The parameter $d$ signifies the natural mortality rate of the prey, and $a$ is the intra-specific competition coefficient among prey. The expression $\frac{pu}{1+qu}$ denotes the Holling-II functional response of the predator. The term $s$ represents

the predator’s intrinsic growth rate, and $-\frac{m\left(x\right)v^2}{r+u}$ indicates intra-specific competition among predators, where $m\left(x\right)\ge 0$ is a spatially dependent coefficient and $r>0$ is a half-saturation constant for the prey-dependent carrying capacity. The term $\frac{cpuv}{1+qu}$ represents biomass conversion, with $c>0$ being the conversion efficiency. The homogeneous Neumann boundary conditions signify that the system is isolated at the boundary.

The novelty of this work with respect to the existing literature is twofold. While studies like [6] incorporate prey-taxis into the Holling-Tanner framework, and [8] -[12] analyze predator-taxis models, our work uniquely couples predator-taxis and fear effects through the same sensitivity coefficient $\alpha$ , which appears simultaneously in the cross-diffusion term $-\alpha \nabla \cdot \left(u\nabla v\right)$ and the fear-modified growth rate $\frac{r_{0}u}{1+k_0\alpha v}$ . This dual coupling creates a strongly nonlinear mathematical structure absent in all the aforementioned works. More importantly, unlike the constant-coefficient competition terms in [6] [8]-[12], we introduce a spatially degenerate intra-specific competition coefficient $m\left(x\right)$ for the predator. This degeneracy, characterized by $m\left(x\right)\equiv 0$ in a subregion ${\Omega }_0$ , models a refuge scenario where predator competition completely vanishes locally. The vanishing of $m\left(x\right)$ fundamentally alters the mathematical framework: it does not merely perturb a uniform state but gives rise to a new class of semi-trivial solutions and bifurcation structures, posing significant analytical challenges in establishing a priori bounds—a task that is central to our present analysis.

For the spatially dependent parameter $m\left(x\right)$ , we assume it exhibits degeneracy within the region Ω, specifically that $m\left(x\right)$ satisfies

$m\left(x\right)\equiv 0,x\in {\overline{\Omega }}_0,m\left(x\right)>0,x\in \overline{\Omega }\backslash {\overline{\Omega }}_0$(3)

where ${\overline{\Omega }}_0$ is a smooth closed subdomain within Ω. From a biological perspective, condition (3) indicates that ${\Omega }_0$ is a suitable subregion for the survival of the predator population, as the intra-specific competition pressure among predators completely vanishes in the subregion ${\Omega }_0$ .

The strong coupling characteristics of the model (2) primarily stem from two key factors: first, the system includes the cross-diffusion term $\alpha \nabla \cdot \left(u\nabla v\right)$ ; and second, the taxis coefficient $\alpha$ appears simultaneously in the reaction term on the right-hand side of the equations. This dual association renders the model (2) a strongly coupled partial differential equation system, posing significant challenges for subsequent theoretical analysis.

The main objective of this paper is to investigate the uniform boundedness, local bifurcation structure, and local asymptotic stability of the positive steady-state solutions of the aforementioned model. Specifically, we first establish a priori estimates for the positive solutions, from which we prove the uniform boundedness of the positive steady-state solutions. Next, taking the growth rate of the prey population $r_0$ as a bifurcation parameter, we utilize the Crandall-Rabinowitz local bifurcation theorem to analyze the existence of positive solution branches bifurcating from the semi-trivial solution curve, and further discuss the linear stability of the branching solutions.

The primary contribution of this paper is a rigorous theoretical analysis, establishing uniform boundedness, existence via bifurcation, and local asymptotic stability for the positive steady-state solutions. We emphasize that the analytical thresholds derived herein—particularly the bifurcation parameter $r_0$ and the stability criteria—are not merely abstract mathematical conditions; ecologically, they delineate critical regimes for species coexistence. Specifically, when the prey growth rate $r_0$ crosses the bifurcation point, the system transitions from a predator-free state to a coexistence state, thereby predicting the precise conditions under which both species can persist in the heterogeneous landscape. The local stability or instability of the bifurcating solution translates directly into ecological forecasts: a stable positive steady state corresponds to the formation of a stationary spatial pattern of coexistence, whereas instability signals the onset of cyclical population dynamics or localized extinctions, especially near the boundary of the predator refuge ${\Omega }_0$ where the degeneracy of $m\left(x\right)$ takes effect. We acknowledge that this analytical work serves as a necessary theoretical foundation for subsequent numerical investigations. Without the a priori estimates and stability thresholds established here, any computational exploration of pattern formation would lack a solid mathematical basis, particularly given the degeneracy of the competition coefficient $m\left(x\right)$ in the refuge region. The numerical visualization of these spatial patterns and the ecological implications of parameter variations will be addressed in our future work.

2. Existence and a Priori Estimates of Positive Solutions

First, we define several symbols that will be used in the subsequent sections of this paper. Let $T$ be a differential operator, $\varphi$ be a Hölder continuous function on Ω, and $d$ be a positive constant. Denote ${\lambda }_1^D\left(\pm d,T+\varphi ,\Omega \right)$ and ${\lambda }_1^N\left(\pm d,T+\varphi ,\Omega \right)$ as the first eigenvalues of the operator $\mp d\text{Δ}+T+\varphi$ in the domain Ω under Dirichlet and Neumann boundary conditions, respectively. When $T=\varphi =0$ , we omit the operator $T$ and the function $\varphi$ , denoting ${\lambda }_1^D\left(\pm d,\Omega \right)={\lambda }_1^D\left(\pm d,0,\Omega \right)$ and ${\lambda }_1^N\left(\pm d,\Omega \right)={\lambda }_1^N\left(\pm d,0,\Omega \right)$ ; when only $T=0$ , we omit the operator $T$ , denoting ${\lambda }_1^D\left(\pm d,\varphi ,\Omega \right)={\lambda }_1^D\left(\pm d,0+\varphi ,\Omega \right)$ and ${\lambda }_1^N\left(\pm d,\varphi ,\Omega \right)={\lambda }_1^N\left(\pm d,0+\varphi ,\Omega \right)$ .

Furthermore, the norms in the spaces $L^p\left(\Omega \right)$ ($p\ge 1$ ) and $C\left(\overline{\Omega }\right)$ are defined as follows

${\Vert \varphi \Vert }_p={\left({\displaystyle {\int }_{\Omega }{\left|\varphi \left(x\right)\right|}^p\text{d}x}\right)}^{1/p},{\Vert \varphi \Vert }_{\infty }=\underset{\overline{\Omega }}{\max }\left|\varphi \left(x\right)\right|.$

According to Theorem 2.2 (i) in [18], when $s<{\lambda }_1^D\left(d_2,{\Omega }_0\right)$ , the boundary value problem

$\{\begin{array}{ll}-d_2\Delta v=sv-\frac{m\left(x\right)v^2}{r},\hfill & x\in \Omega ,\hfill \\ {\partial }_{n}v=0,\hfill & x\in \partial \Omega \hfill \end{array}$ (4)

has a unique positive solution $v_s^{d_2}\left(x\right)$ ; when $s\ge {\lambda }_1^D\left(d_2,{\Omega }_0\right)$ , there are no positive solutions to this problem, where the function $m\left(x\right)$ satisfies condition (3). This indicates that when $0<s<{\lambda }_1^D\left(d_2,{\Omega }_0\right)$ , the model (2) has a semi-trivial solution $\left(0,v_s^{d_2}\right)$ . Additionally, when $r_0>d$ , the model (2) also has another semi-trivial solution $\left(\frac{r_0-d}{a},0\right)$ , and for any parameter values, the model always has the trivial solution $\left(0,0\right)$ .

Next, we will prove that

$s<{\lambda }_1^D\left(d_2,{\Omega }_0\right)$ (5)

is a necessary condition for the existence of positive solutions of the model (2), as stated in the following lemma.

Lemma 2.1 The model (2) has positive solutions only if (5) holds.

Proof. We use proof by contradiction. Assume that when $s\ge {\lambda }_1^D\left(d_2,{\Omega }_0\right)$ , the model (2) has a positive solution $\left(u,v\right)$ , then $v$ must be an upper solution of the boundary value problem (4). Let $\phi$ be the eigenfunction corresponding to the eigenvalue ${\lambda }_1^D\left(d_2,{\Omega }_0\right)$ . Extending $\phi$ to zero in the entire domain Ω, for sufficiently small $ϵ>0$ , $ϵ\phi$ serves as a lower solution for the boundary value problem (4). By the method of upper and lower solutions for elliptic equations, there must exist a non-negative solution between the lower solution $ϵ\phi$ and the upper solution $v$ , which is strictly positive in the interior of the domain Ω. This indicates that the boundary value problem (4) has a positive solution, contradicting the conclusion of Theorem 2.2 (i) in [18] (which states that when $s\ge {\lambda }_1^D\left(d_2,{\Omega }_0\right)$ , the boundary value problem (4) has no positive solutions). Therefore, the assumption is false, and the original conclusion is proved.

From Lemma 2.1, we conclude that condition (5) is necessary for the existence of positive solutions of the model (2). Hence, in the subsequent content of this paper, we will always assume that condition (5) holds.

Theorem 2.2 Let $a,M,d$ be given positive constants satisfying the condition $s<{\lambda }_1^D\left(d_2,{\Omega }_0\right)$ . Then, when the parameter $r_0$ satisfies $d<r_0<aM+d$ , there exists a positive constant $W=W\left(M\right)$ that depends only on $M$ , such that any positive solution $\left(u,v\right)$ of the model (2) satisfies ${\Vert u\Vert }_{\infty }+{\Vert v\Vert }_{\infty }\le W$ .

Proof. We employ a proof by contradiction. Assume that the conclusion does not hold. Then there exists a sequence ${\left\{r_0^n\right\}}_{n=1}^{\infty }$ such that $d<r_0^n<aM+d$ , for which there exist positive solutions $\left(u_n,v_n\right)$ of the model (2) satisfying

${\Vert u_n\Vert }_{\infty }+{\Vert v_n\Vert }_{\infty }\to \infty ,n\to \infty .$ (6)

From the first equation of the model (2), $u_n$ satisfies

$-d_1\Delta u_n-\alpha \nabla \cdot \left(u_n\nabla v_n\right)=\frac{r_0^n{u}_n}{1+k_0\alpha v_n}-d{u}_n-a{u}_n^2-\frac{p{u}_n{v}_n}{1+q{u}_n},x\in \Omega .$

Noting that $\frac{r_0^n{u}_n}{1+k_0\alpha v_n}\le r_0^n{u}_n$ and $\frac{p{u}_n{v}_n}{1+q{u}_n}\ge 0$ , we obtain

$-d_1\text{Δ}{u}_n-\alpha \nabla \cdot \left(u_n\nabla v_n\right)\le r_0^n{u}_n-d{u}_n-a{u}_n^2.$

Let $U_n$ be the positive solution of the ordinary differential equation $r_0^{n}U-dU-a{U}^2=0$ , which is given by $U_n=\frac{r_0^n-d}{a}$ . By the maximum principle for elliptic equations, we find that

$u_n\left(x\right)\le U_n=\frac{r_0^n-d}{a}\le \frac{aM+d-d}{a}=M,\forall x\in \Omega .$

Thus, ${\Vert u_n\Vert }_{\infty }\le M$ , which indicates that $\left\{u_n\right\}$ is uniformly bounded. Combining this with equation (6), we conclude that ${\Vert v_n\Vert }_{\infty }\to \infty$ as $n\to \infty$ .

Define ${\widehat{v}}_n=\frac{v_n}{{\Vert v_n\Vert }_{\infty }}$ . Dividing both sides of the second equation of the model (2) by ${\Vert v_n\Vert }_{\infty }$ , we find that ${\widehat{v}}_n$ satisfies

$-d_2\Delta {\widehat{v}}_n={\widehat{v}}_n\left(s-\frac{m\left(x\right){\widehat{v}}_n{\Vert v_n\Vert }_{\infty }}{r+u_n}+\frac{cp{u}_n}{1+q{u}_n}\right),x\in \Omega ,{\partial }_{\mu }{\widehat{v}}_n=0,x\in \partial \Omega .$(7)

Since $m\left(x\right)\ge 0$ , ${\widehat{v}}_n\ge 0$ , ${\Vert v_n\Vert }_{\infty }>0$ , and $r+u_n>0$ , it follows that

$-\frac{m\left(x\right){\widehat{v}}_n{\Vert v_n\Vert }_{\infty }}{r+u_n}\le 0.$

Moreover, since $\frac{u_n}{1+q{u}_n}<\frac{1}{q}$ , we have

$-d_2\text{Δ}{\widehat{v}}_n<\left(s+\frac{cp}{q}\right){\widehat{v}}_n.$

Multiplying both sides of this inequality by ${\widehat{v}}_n$ and integrating over Ω using integration by parts yields

$d_2{\displaystyle {\int }_{\Omega }{\left|\nabla {\widehat{v}}_n\right|}^2\text{d}x}<\left(s+\frac{cp}{q}\right){\displaystyle {\int }_{\Omega }{\widehat{v}}_n^2\text{d}x}.$

Adding ${\displaystyle {\int }_{\Omega }{\widehat{v}}_n^2\text{d}x}$ to both sides, and noting that ${\Vert {\widehat{v}}_n\Vert }_{\infty }=1$ implies ${\displaystyle {\int }_{\Omega }{\widehat{v}}_n^2\text{d}x}\le \left|\Omega \right|$, we obtain

$d_2{\displaystyle {\int }_{\Omega }{\left|\nabla {\widehat{v}}_n\right|}^2\text{d}x}+{\displaystyle {\int }_{\Omega }{\widehat{v}}_n^2\text{d}x}\le \left(s+\frac{cp}{q}+1\right)\left|\Omega \right|.$

This shows that the sequence ${\left\{{\widehat{v}}_n\right\}}_{n=1}^{\infty }$ is uniformly bounded in $H^1\left(\Omega \right)$ . Therefore, there exists a subsequence (still denoted as $\left\{{\widehat{v}}_n\right\}$ ) and a function $\widehat{v}\in H^1\left(\Omega \right)$ such that

${\widehat{v}}_n\rightharpoonup \widehat{v}\text{weakly}\text{in}{H}^1\left(\Omega \right),{\widehat{v}}_n\to \widehat{v}\text{strongly}\text{in}{L}^2\left(\Omega \right).$

From ${\Vert {\widehat{v}}_n\Vert }_{\infty }=1$ and the strong convergence, using the argument in Proposition 4.1 from [18], we can deduce that $\widehat{v}\cancel{\equiv }0$ holds in $\overline{\Omega }$ , and it is evident that $0\le \widehat{v}\le 1$ .

Next, we will prove that $\widehat{v}=0$ holds almost everywhere in $\Omega \backslash {\overline{\Omega }}_0$ . Multiplying both sides of equation (7) by ${\widehat{v}}_n$ and integrating over Ω gives us

${\Vert v_n\Vert }_{\infty }{\displaystyle {\int }_{\Omega }\frac{m\left(x\right)}{r+u_n}{\widehat{v}}_n^3\text{d}x}={\displaystyle {\int }_{\Omega }\left(s+\frac{cp{u}_n}{1+q{u}_n}\right){\widehat{v}}_n^2\text{d}x}-d_2{\displaystyle {\int }_{\Omega }{\left|\nabla {\widehat{v}}_n\right|}^2\text{d}x}.$(8)

Since $\left\{{\widehat{v}}_n\right\}$ is bounded in $H^1\left(\Omega \right)$ and ${\Vert u_n\Vert }_{\infty }\le M$ , the right-hand side of Equation (8) is uniformly bounded with respect to $n$ . Given that ${\Vert v_n\Vert }_{\infty }\to \infty$ , it follows that the left-hand side must tend to zero

$\underset{n\to \infty }{\text{lim}}{\displaystyle {\int }_{\Omega }\frac{m\left(x\right)}{r+u_n}{\widehat{v}}_n^3\text{d}x}=0.$

Due to the boundedness of $\left\{u_n\right\}$ , there exists a subsequence such that $u_n\to u^{\text{*}}$ weakly in $L^2\left(\Omega \right)$ . Combining this with the strong convergence of ${\widehat{v}}_n$ , we take the limit to obtain

${\displaystyle {\int }_{\Omega }\frac{m\left(x\right)}{r+u^{\text{*}}}{\widehat{v}}^3\text{d}x}=0.$

Since $m\left(x\right)>0$ in $\Omega \backslash {\overline{\Omega }}_0$ and $r+u^{\text{*}}\ge r>0$ , we conclude that $\widehat{v}=0$ holds almost everywhere in $\Omega \backslash {\overline{\Omega }}_0$ .

Now consider the region ${\Omega }_0$ . Since $m\left(x\right)\equiv 0$ on ${\overline{\Omega }}_0$ , Equation (7) simplifies in ${\Omega }_0$ to

$-d_2\Delta {\widehat{v}}_n=\left(s+\frac{cp{u}_n}{1+q{u}_n}\right){\widehat{v}}_n,x\in {\Omega }_0.$

Multiplying this equation by a test function $\phi \in C_0^{\infty }\left({\Omega }_0\right)$ and integrating over ${\Omega }_0$ , we obtain

$d_2{\displaystyle {\int }_{{\Omega }_0}\nabla {\widehat{v}}_n\cdot \nabla \phi \text{d}x}={\displaystyle {\int }_{{\Omega }_0}\left(s+\frac{cp{u}_n}{1+q{u}_n}\right){\widehat{v}}_n\phi \text{d}x}.$

Taking the limit as $n\to \infty$ and using the weak convergence of ${\widehat{v}}_n$ along with the boundedness of $u_n$ , we find that multiplying both sides by the test function $\phi \in C_0^{\infty }\left({\Omega }_0\right)$ yields

$d_2{\displaystyle {\int }_{{\Omega }_0}\nabla \widehat{v}\cdot \nabla \phi \text{d}x}={\displaystyle {\int }_{{\Omega }_0}\left(s+\frac{cp{u}^{\text{*}}}{1+q{u}^{\text{*}}}\right)\widehat{v}\phi \text{d}x}.$

This indicates that ${\widehat{v}|}_{{\Omega }_0}\in H_0^1\left({\Omega }_0\right)$ is a non-negative weak solution of the equation

$-d_2\Delta \widehat{v}=\left(s+\frac{cp{u}^{*}}{1+q{u}^{*}}\right)\widehat{v},x\in {\Omega }_0.$

Since $\widehat{v}\cancel{\equiv }0$ and $\widehat{v}=0$ in $\Omega \backslash {\overline{\Omega }}_0$ , it must be that ${\widehat{v}|}_{{\Omega }_0}\cancel{\equiv }0$ . By the Krein–Rutman theorem, $s$ must be the principal eigenvalue of the operator $-d_2\text{Δ}-\frac{cp{u}^{\text{*}}}{1+q{u}^{\text{*}}}$ on ${\Omega }_0$ with homogeneous Dirichlet boundary conditions, i.e.,

$s={\lambda }_1^D\left(d_2,-\frac{cp{u}^{*}}{1+q{u}^{*}},{\Omega }_0\right).$

Since $u^{\text{*}}\ge 0$ in ${\Omega }_0$ , we have $-\frac{cp{u}^{\text{*}}}{1+q{u}^{\text{*}}}\le 0$ . By the monotonicity of eigenvalues with respect to potential functions, we obtain

${\lambda }_1^D\left(d_2,-\frac{cp{u}^{*}}{1+q{u}^{*}},{\Omega }_0\right)\ge {\lambda }_1^D\left(d_2,0,{\Omega }_0\right)={\lambda }_1^D\left(d_2,{\Omega }_0\right).$

Thus, $s\ge {\lambda }_1^D\left(d_2,{\Omega }_0\right)$ , which contradicts the assumption of the theorem that $s<{\lambda }_1^D\left(d_2,{\Omega }_0\right)$ . This contradiction indicates that the assumption made in the proof by contradiction is false, thereby proving the original statement.

3. Local Bifurcation Structure and Stability of Local Bifurcation Solutions

In this subsection, we take $r_0$ as the bifurcation parameter and, using the local bifurcation theory proposed by Crandall and Rabinowitz [19], we study the local bifurcation structure of positive solutions of model (2) starting from the semi-trivial solutions $\left(\frac{r_0-d}{a},0\right)$ and $\left(0,v_s^{d_2}\right)$ , as well as the stability of these local bifurcation solutions.

Now, we investigate the local bifurcation phenomena of model (2). First, we introduce the function

$\tau \left(r_0\right)={\lambda }_1^N\left(-d_1,\alpha \nabla v_s^{d_2}\nabla +\alpha \Delta v_s^{d_2}+\frac{r_0}{1+k_0\alpha v_s^{d_2}}-p{v}_s^{d_2},\Omega \right),$

where $0<s<{\lambda }_1^D\left(d_2,{\Omega }_0\right)$ and $r_0\in \left[d,+\infty \right)$ . Since $\frac{1}{1+k_0\alpha v_s^{d_2}}>0$ , we know that $\tau \left(r_0\right)$ is continuous and strictly increasing with respect to $r_0$ . Using a reasoning similar to Theorem 2.1 in [20] and combining standard regularity theory for elliptic equations, Sobolev embedding theorem, and the boundary smoothness assumption $\partial \Omega \in C^{2+ϵ}$ $\left(0<ϵ<1\right)$ , we can show that there exists a constant $M^{\circ }$ such that any positive solution $\left(u,v\right)\in C^2\left(\Omega \right)\times C^2\left(\Omega \right)$ of model (2) satisfies

$\underset{\overline{\Omega }}{\max }\left\{{\Vert \nabla u\Vert }_{C^1},{\Vert \nabla v\Vert }_{C^1},{\Vert \Delta v\Vert }_{C^1}\right\}\le M^{\circ }.$

Furthermore, assuming there exists a parameter $d$ such that $\tau \left(d\right)<d$ , regardless of whether $\tau \left(d\right)<0$ , $\tau \left(d\right)=0$ , or $\tau \left(d\right)>0$ , there exists a unique $r_0^{\text{**}}\in \left[d,+\infty \right)$ such that

$\tau \left(r_0^{**}\right)={\lambda }_1^N\left(-d_1,\alpha \nabla v_s^{d_2}\nabla +\alpha \Delta v_s^{d_2}+\frac{r_0}{1+k_0\alpha v_s^{d_2}}-p{v}_s^{d_2},\Omega \right)=d.$ (9)

Let $p>N$ , define the function spaces $X=W_n^{2,p}\left(\Omega \right)\times W_n^{2,p}\left(\Omega \right)$ and $Y=L^p\left(\Omega \right)\times L^p\left(\Omega \right)$ , where

$W_n^{2,p}\left(\Omega \right)=\left\{u\in W^{2,p}\left(\Omega \right):{\partial }_{n}u=0,x\in \partial \Omega \right\}.$

Model (2) has two semi-trivial solution curves ${\Theta }_u$ and ${\Theta }_v$ , which are given by

${\Theta }_u=\left\{\left(r_0,u,v\right):\left(r_0,\left(r_0-d\right)/a,0\right):r_0>d\right\},{\Theta }_v=\left\{\left(r_0,u,v\right):\left(r_0,0,v_s^{d_2}\right):r_0>0\right\}.$

Next, we separately study the local bifurcation structure of positive solutions of model (2) from the semi-trivial solution curves ${\Theta }_u$ and ${\Theta }_v$ .

Theorem 3.1 Let $r_0^{*}=d-\frac{sa}{cp+sq}$ , and $r_0^{**}$ be given by (9). The following conclusions hold

• If $-\frac{cp}{q}<s<0$ , then the sufficient and necessary condition for the model (2) to produce a positive solution branch from the semi-trivial solution branch ${\Theta }_u$ is $r_0=r_0^{\text{*}}$ . Furthermore, all positive solutions of model (2) in the neighborhood of the point $\left(r_0^{\text{*}},\frac{r_0^{\text{*}}-d}{a},0\right)\in ℝ\times X$ lie on a smooth curve ${\theta }_1$ , expressed as

${\theta }_1=\left\{\left(r_0\left(\theta \right),u\left(\theta \right),v\left(\theta \right)\right)=\left(r_0^{*}+\theta {r^{\prime }}_0\left(0\right)+o\left(\left|\theta \right|\right),\frac{r_0^{*}-d}{a}+\theta {\varphi }_{*}+o\left(\left|\theta \right|\right),\theta {\psi }_{*}+o\left(\left|\theta \right|\right)\right):0<\theta <{\theta }_1\right\},$

where ${\theta }_1>0$ is a constant, ${\varphi }_{\text{*}}$ and ${\psi }_{\text{*}}$ are determined functions in the space $W_n^{2,p}\left(\Omega \right)$ , and ${r^{\prime }}_0\left(0\right)$ is defined as

${r^{\prime }}_0\left(0\right)=\frac{{\displaystyle {\int }_{\Omega }\left(\frac{m\left(x\right)}{r+\frac{r_0^{\text{*}}-d}{a}}{\psi }_{\text{*}}^3-\frac{cp}{{\left(1+q\frac{r_0^{\text{*}}-d}{a}\right)}^2}{\varphi }_{\text{*}}{\psi }_{\text{*}}^2\right)\text{d}x}}{{\displaystyle {\int }_{\Omega }\frac{cpa}{{\left(a+q\left(r_0^{\text{*}}-d\right)\right)}^2}{\psi }_{\text{*}}^2\text{d}x}}.$

Additionally, when $p<\tilde{p}$ , the solution branch ${\theta }_1$ at the point $\left(r_0^{\text{*}},\frac{r_0^{\text{*}}-d}{a},0\right)$ undergoes a supercritical bifurcation; when $p>\tilde{p}$ , this bifurcation is subcritical, where the critical parameter $\tilde{p}$ satisfies

$\tilde{p}=\frac{{\displaystyle {\int }_{\Omega }\frac{m\left(x\right){\left(1+q\frac{r_0^{\text{*}}-d}{a}\right)}^2}{r+\frac{r_0^{\text{*}}-d}{a}}{\psi }_{\text{*}}^3\text{d}x}}{{\displaystyle {\int }_{\Omega }c{\varphi }_{\text{*}}{\psi }_{\text{*}}^2\text{d}x}}.$

• If $0<s<{\lambda }_1^D\left(d_2,{\Omega }_0\right)$ , then the sufficient and necessary condition for the model (2) to produce a positive solution branch from the semi-trivial solution branch ${\Theta }_v$ is $r_0=r_0^{\text{**}}$ . Furthermore, all positive solutions of model (2) in the neighborhood of the point $\left(r_0^{\text{**}},0,v_s^{d_2}\right)\in ℝ\times X$ lie on a smooth curve ${\theta }_2$ , expressed as

${\theta }_2=\left\{\left(r_0\left(\theta \right),u\left(\theta \right),v\left(\theta \right)\right)=\left(r_0^{\text{**}}+\theta {{\widehat{r}}^{\prime }}_0\left(0\right)+o\left(\left|\theta \right|\right),\theta {\varphi }_{\text{**}}+o\left(\left|\theta \right|\right),v_s^{d_2}+\theta {\psi }_{\text{**}}+o\left(\left|\theta \right|\right)\right):0<\theta <{\theta }_2\right\},$

where ${\theta }_2>0$ is a constant, ${\varphi }_{\text{**}}$ and ${\psi }_{\text{**}}$ are determined functions in the space $W_n^{2,p}\left(\Omega \right)$ , and ${{\widehat{r}}^{\prime }}_0\left(0\right)$ is defined as

$\begin{array}{c}{{\widehat{r}}^{\prime }}_0\left(0\right)=({\displaystyle {\int }_{\Omega }(\left(a-pq{v}_s^{d_2}\right){\varphi }_{**}^3-\left(\alpha \nabla {\varphi }_{**}\nabla {\psi }_{**}+\alpha {\varphi }_{**}\Delta {\psi }_{**}\right){\varphi }_{**}\begin{array}{c}\\ \\ \end{array}}\\ +\left(\frac{r_0^{**}{k}_0\alpha }{{\left(1+k_0\alpha v_s^{d_2}\right)}^2}+p\right){\varphi }_{**}^2{\psi }_{**})\text{d}x)/{\displaystyle {\int }_{\Omega }\frac{{\varphi }_{**}^2}{1+k_0\alpha v_s^{d_2}}\text{d}x}.\end{array}$

Additionally, when $a>\tilde{a}$ , the solution branch ${\theta }_2$ at the point $\left(r_0^{\text{**}},0,v_s^{d_2}\right)$ undergoes a supercritical bifurcation; when $a<\tilde{a}$ , this bifurcation is subcritical, where the critical parameter $\tilde{a}$ satisfies

$\begin{array}{c}\tilde{a}={\displaystyle {\int }_{\Omega }(pq{v}_s^{d_2}{\varphi }_{**}^3+\left(\alpha \nabla {\varphi }_{**}\nabla {\psi }_{**}+\alpha {\varphi }_{**}\Delta {\psi }_{**}\right){\varphi }_{**}\begin{array}{c}\\ \\ \end{array}}\\ -\left(p+\frac{r_0^{**}{k}_0\alpha }{{\left(1+k_0\alpha v_s^{d_2}\right)}^2}\right){\varphi }_{**}^2{\psi }_{**})\text{d}x/{\displaystyle {\int }_{\Omega }{\varphi }_{**}^3\text{d}x}.\end{array}$

Proof. The proof of conclusion (i) is completely similar to that of (ii), so we will only provide a detailed proof of conclusion (ii) here. Let $w=v-v_s^{d_2}$ , and define the nonlinear operator $L_1:ℝ\times X\to Y$ as

$L_1\left(r_0,u,w\right)={\left(f_1,f_2\right)}^{\top },$

where the component expressions are given by

$\begin{array}{l}{f}_1=d_1\Delta u+\alpha \nabla u\nabla v+\alpha u\Delta v+\frac{r_{0}u}{1+k_0\alpha v}-du-a{u}^2-\frac{puv}{1+qu},\\ f_2=d_2\Delta v+sv-\frac{m\left(x\right)v^2}{r+u}+\frac{cpuv}{1+qu}.\end{array}$

Clearly, the sufficient and necessary condition for the equation $L_1\left(r_0,u,w\right)=0$ to hold is that $\left(u,w+v_s^{d_2}\right)$ is a solution of model (2). Therefore, to find the positive solution branch of (2) near $\left(r_0,u,v\right)=\left(r_0,0,v_s^{d_2}\right)$ , we only need to find the positive solution branch of $L_1\left(r_0,u,w\right)=0$ near $\left(r_0,u,w\right)=\left(r_0,0,0\right)$ .

Let $L_{1\left(u,w\right)}\left(r_0,u,w\right)$ denote the Fréchet derivative of the operator $L_1$ with respect to the variables $\left(u,w\right)$ , and direct calculation gives

$L_{1\left(u,w\right)}\left(r_0,u,w\right){\left(\varphi ,\psi \right)}^{\top }=\left(\begin{array}{c}{B}_1\left(\varphi ,\psi \right)+B_2\left(\varphi ,\psi \right)\\ d_2\text{Δ}+B_3\left(\varphi ,\psi \right)\end{array}\right)$

where the coefficient operators $B_1,B_2,B_3$ are defined as

$\begin{array}{l}{B}_1=\frac{r_0}{1+k_0\alpha \left(w+v_s^{d_2}\right)}\varphi -d\varphi -2au\varphi -\frac{p\left(w+v_s^{d_2}\right)}{{\left(1+qu\right)}^2}\varphi -\frac{r_0{k}_0\alpha u}{{\left(1+k_0\alpha \left(w+v_s^{d_2}\right)\right)}^2}\psi -\frac{pu}{1+qu}\psi ,\\ B_2=d_1\Delta \varphi +\alpha \nabla \left(w+v_s^{d_2}\right)\nabla \varphi +\alpha \Delta \left(w+v_s^{d_2}\right)\varphi +\alpha \nabla u\nabla \psi +\alpha u\Delta \psi ,\\ B_3=\frac{m\left(x\right){\left(w+v_s^{d_2}\right)}^2}{{\left(r+u\right)}^2}\varphi +\frac{cp\left(w+v_s^{d_2}\right)}{{\left(1+qu\right)}^2}\varphi +s\psi -\frac{2m\left(x\right)\left(w+v_s^{d_2}\right)}{r+u}\psi +\frac{cpu}{1+qu}\psi ,\end{array}$

with $\varphi ,\psi \in W_n^{2,p}\left(\Omega \right)$ .

Based on the Krein-Rutman theorem, the existence of a positive solution $\left(\varphi ,\psi \right)$ satisfying $\varphi >0$ for the equation $L_{1\left(u,w\right)}\left(r_0,0,0\right){\left(\varphi ,\psi \right)}^{\top }={\left(0,0\right)}^{\top }$ is equivalent to the condition $r_0=r_0^{\text{**}}$ . Furthermore, it can be verified that the point $\left(r_0^{\text{**}},0,v_s^{d_2}\right)$ is the unique bifurcation point from the semi-trivial solution branch ${\Theta }_v$ of model (2) producing a positive solution branch.

Let $\left({\varphi }_{**},{\psi }_{**}\right)\in \ker L_{1\left(u,w\right)}\left(r_0^{**},0,0\right)$ , then $\left({\varphi }_{\text{**}},{\psi }_{\text{**}}\right)$ satisfies

$\{\begin{array}{ll}{d}_1\Delta {\varphi }_{**}+\alpha \nabla v_s^{d_2}\nabla {\varphi }_{**}+\alpha {\varphi }_{**}\Delta v_s^{d_2}+\frac{r_0{\varphi }_{**}}{1+k_0\alpha v_s^{d_2}}-d{\varphi }_{**}-p{v}_s^{d_2}{\varphi }_{**}=0,\hfill & x\in \Omega ,\hfill \\ d_2\Delta {\psi }_{**}+s{\psi }_{**}-\frac{2m\left(x\right)v_s^{d_2}{\psi }_{**}}{r}+\frac{m\left(x\right){\left(v_s^{d_2}\right)}^2}{r^2}{\varphi }_{**}+cp{v}_s^{d_2}{\varphi }_{**}=0,\hfill & x\in \Omega ,\hfill \\ {\partial }_n{\varphi }_{**}={\partial }_n{\psi }_{**}=0,\hfill & x\in \partial \Omega .\hfill \end{array}$ (10)

First, we analyze the first equation of (10). Based on this equation, we construct the following Neumann boundary value eigenvalue problem

$\{\begin{array}{ll}{d}_1\Delta {\varphi }_{**}+\alpha \nabla v_s^{d_2}\nabla {\varphi }_{**}+\alpha {\varphi }_{**}\Delta v_s^{d_2}+\frac{r_0{\varphi }_{**}}{1+k_0\alpha v_s^{d_2}}-d{\varphi }_{**}-p{v}_s^{d_2}{\varphi }_{**}=0,\hfill & x\in \Omega ,\hfill \\ {\partial }_n{\varphi }_{**}=0,\hfill & x\in \partial \Omega ,\hfill \end{array}$

According to the previously defined notation, the principal eigenvalue of the above eigenvalue problem can be denoted as

$\xi \left(r_0\right)={\lambda }_1^N\left(-d_1,\alpha \nabla v_s^{d_2}\nabla +\alpha \Delta v_s^{d_2}+\frac{r_0}{1+k_0\alpha v_s^{d_2}}-p{v}_s^{d_2},\Omega \right)=d,$

where the parameters satisfy $0<s<{\lambda }_1^D\left(d_2,{\Omega }_0\right)$ and $r_0\in \left(d,+\infty \right)$ . From the monotonicity of the principal eigenvalue of the elliptic operator, we know that the principal eigenvalue strictly increases with the potential function. Since the coefficient of $r_0$ , $\frac{1}{1+k_0\alpha v_s^{d_2}}>0$ , it follows that $\xi \left(r_0\right)$ is strictly increasing with respect to $r_0$ , and as $r_0\to +\infty$ , $\xi \left(r_0\right)\to +\infty$ . Combining $\xi \left(r_0^{\text{**}}\right)=d$ , we conclude that the first equation of (10) has a unique positive solution ${\varphi }_{**}\in W_n^{2,p}\left(\Omega \right)$ if and only if $r_0=r_0^{\text{**}}$ . Next, we analyze the second equation of (10). From equation (4), we know that

$\left(d_2\Delta +s-\frac{m\left(x\right)v_s^{d_2}}{r}\right)v_s^{d_2}=0\cdot v_s^{d_2},$

which indicates that ${\lambda }_1^N\left(d_2\Delta +s-\frac{m\left(x\right)v_s^{d_2}}{r},\Omega \right)=0$ . Again using the monotonicity of the principal eigenvalue, since the potential function satisfies $-\frac{2m\left(x\right)v_s^{d_2}}{r}<-\frac{m\left(x\right)v_s^{d_2}}{r}$ , we have

${\lambda }_1^N\left(d_2\Delta +s-\frac{2m\left(x\right)v_s^{d_2}}{r},\Omega \right)<{\lambda }_1^N\left(d_2\Delta +s-\frac{m\left(x\right)v_s^{d_2}}{r},\Omega \right)=0.$

For the elliptic operator Δ, its principal eigenvalue is the unique maximum eigenvalue among all eigenvalues. Therefore, all eigenvalues of the operator $d_2\Delta +s-\frac{2m\left(x\right)v_s^{d_2}}{r}$ are less than 0, indicating that this operator is strictly negative, and consequently, its inverse ${\left(d_2\Delta +s-\frac{2m\left(x\right)v_s^{d_2}}{r}\right)}^{-1}$ exists and is strictly negative. Thus, the second equation of (10) can be equivalently transformed into

${\psi }_{**}=-{\left(d_2\Delta +s-\frac{2m\left(x\right)v_s^{d_2}}{r}\right)}^{-1}\left(\frac{m\left(x\right){\left(v_s^{d_2}\right)}^2}{r^2}+cp{v}_s^{d_2}\right){\varphi }_{**}>0.$

Since ${\varphi }_{\text{**}}>0$ and $\frac{m\left(x\right){\left(v_s^{d_2}\right)}^2}{r^2}+cp{v}_s^{d_2}>0$ , combined with the negativity of the inverse operator, we conclude that ${\psi }_{\text{**}}>0$ . In summary, we can prove that $\text{ker}\left(L_{1\left(u,w\right)}\left(r_0^{\text{**}},0,0\right)\right)=\text{Span}\left\{\left({\varphi }_{\text{**}},{\psi }_{\text{**}}\right)\right\}$ , indicating that the dimension of this kernel space is

$\text{dim}\left(\text{ker}\left(L_{1\left(u,w\right)}\left(r_0^{\text{**}},0,0\right)\right)\right)=1.$

Next, we verify that $\text{codim}\left(\text{Range}\left(L_{1\left(u,w\right)}\left(r_0^{\text{**}},0,0\right)\right)\right)=1$ . Consider the adjoint operator $L_{1\left(u,w\right)}^{\text{**}}\left(r_0^{\text{**}},0,0\right)$ . According to functional analysis theory, the kernel space of the adjoint operator is orthogonal to the range of the original operator, i.e., $\text{Range}{\left(L_{1\left(u,w\right)}\left(r_0^{\text{**}},0,0\right)\right)}^{\perp }=\text{ker}\left(L_{1\left(u,w\right)}^{\text{**}}\left(r_0^{\text{**}},0,0\right)\right)$ . The matrix form of the operator $L_{1\left(u,w\right)}^{\text{**}}\left(r_0^{\text{**}},0,0\right)$ can be expressed as

$L_{1\left(u,w\right)}^{\text{**}}\left(r_0^{\text{**}},0,0\right)=\left(\begin{array}{cc}{A}_1& A_2\\ 0& A_3\end{array}\right)$

where

$\begin{array}{l}{A}_1=d_1\Delta +\alpha \nabla v_s^{d_2}\nabla +\alpha \Delta v_s^{d_2}+\frac{r_0}{1+k_0\alpha v_s^{d_2}}-dp{v}_s^{d_2},\\ A_2=\frac{m\left(x\right){\left(v_s^{d_2}\right)}^2}{r^2}+cp{v}_s^{d_2},\\ A_3=d_2\Delta +s-\frac{2m\left(x\right)v_s^{d_2}}{r},\end{array}$

To determine the structure of $\text{ker}\left(L_{1\left(u,w\right)}^{\text{**}}\left(r_0^{\text{**}},0,0\right)\right)$ , let ${\left({{\varphi }^{\prime }}_{\text{**}},{{\psi }^{\prime }}_{\text{**}}\right)}^{\top }\in \text{ker}\left(L_{1\left(u,w\right)}^{\text{**}}\left(r_0^{\text{**}},0,0\right)\right)$ , then it satisfies

$L_{1\left(u,w\right)}^{\text{**}}\left(r_0^{\text{**}},0,0\right){\left({{\varphi }^{\prime }}_{\text{**}},{{\psi }^{\prime }}_{\text{**}}\right)}^{\top }={\left(0,0\right)}^{\top },$

Expanding gives the following Neumann boundary value problem

$\{\begin{array}{ll}{d}_1\Delta {{\varphi }^{\prime }}_{**}+\alpha \nabla v_s^{d_2}\nabla {{\varphi }^{\prime }}_{**}+\alpha {{\varphi }^{\prime }}_{**}\Delta v_s^{d_2}+\frac{r_0{{\varphi }^{\prime }}_{**}}{1+k_0\alpha v_s^{d_2}}-d{{\varphi }^{\prime }}_{**}-p{v}_s^{d_2}{{\varphi }^{\prime }}_{**}\hfill & \hfill \\ +\frac{m\left(x\right){\left(v_s^{d_2}\right)}^2}{r^2}{{\varphi }^{\prime }}_{**}+cp{v}_s^{d_2}{{\varphi }^{\prime }}_{**}=0,\hfill & x\in \Omega ,\hfill \\ d_2\Delta {{\psi }^{\prime }}_{**}+s{{\psi }^{\prime }}_{**}-\frac{2m\left(x\right)v_s^{d_2}{{\psi }^{\prime }}_{**}}{r}=0,\hfill & x\in \Omega ,\hfill \\ {\partial }_n{{\varphi }^{\prime }}_{**}={\partial }_n{{\psi }^{\prime }}_{**}=0,\hfill & x\in \partial \Omega .\hfill \end{array}$ (11)

First, we analyze the second equation of (11). As previously mentioned, the operator $d_2\Delta +s-\frac{2m\left(x\right)v_s^{d_2}}{r}$ is a strictly negative operator, and its inverse exists and is also strictly negative, hence this equation only has the trivial solution ${{\psi }^{\prime }}_{\text{**}}\equiv 0$ . Substituting ${{\psi }^{\prime }}_{\text{**}}\equiv 0$ into the first Equation of (11), we simplify to obtain

$d_1\Delta {{\varphi }^{\prime }}_{**}+\alpha \nabla v_s^{d_2}\nabla {{\varphi }^{\prime }}_{**}+\alpha {{\varphi }^{\prime }}_{**}\Delta v_s^{d_2}+\frac{r_0{{\varphi }^{\prime }}_{**}}{1+k_0\alpha v_s^{d_2}}-d{{\varphi }^{\prime }}_{**}-p{v}_s^{d_2}{{\varphi }^{\prime }}_{**}=0,$

This can be equivalently rewritten in the form of an eigenvalue equation

$\left(d_1\Delta +\alpha \nabla v_s^{d_2}\nabla +\alpha \Delta v_s^{d_2}+\frac{r_0}{1+k_0\alpha v_s^{d_2}}-d-p{v}_s^{d_2}\right)\cdot {{\varphi }^{\prime }}_{**}=0\cdot {{\varphi }^{\prime }}_{**}.$

Here, 0 is the unique eigenvalue of this eigenvalue equation, and the corresponding eigenvector space is one-dimensional. According to the Fredholm alternative theorem, the range of the operator is orthogonal to the kernel space of its adjoint operator

$\text{Range}\left(L_{1\left(u,w\right)}^{\text{**}}\left(r_0^{\text{**}},0,0\right)\right)={\left[\text{ker}\left(L_{1\left(u,w\right)}^{\text{**}}\left(r_0^{\text{**}},0,0\right)\right)\right]}^{\perp }.$

Combining the dimension of the kernel space $\text{dim}\left(\text{ker}\left(L_{1\left(u,w\right)}^{\text{**}}\left(r_0^{\text{**}},0,0\right)\right)\right)=1$ , we obtain

$\text{codim}\left(\text{Range}\left(L_{1\left(u,w\right)}^{\text{**}}\left(r_0^{\text{**}},0,0\right)\right)\right)=\text{dim}\left(\text{ker}\left(L_{1\left(u,w\right)}^{\text{**}}\left(r_0^{\text{**}},0,0\right)\right)\right)=1.$

Next, we verify the transversality condition for the bifurcation. The partial derivative of the operator $L_{1\left(u,w\right)}\left(r_0^{\text{**}},0,0\right)$ with respect to the parameter $r_0$ can be expressed as

$L_{1r_0\left(u,w\right)}\left(r_0^{\text{**}},u,w\right){\left(\varphi ,\psi \right)}^{\top }=\left(\begin{array}{c}\frac{\varphi }{1+k_0\alpha \left(w+v_s^{d_2}\right)}-\frac{k_0\alpha u\psi }{{\left(1+k_0\alpha \left(w+v_s^{d_2}\right)\right)}^2}\\ 0\end{array}\right)$

Letting $u=0$ and $w=0$ , we substitute to obtain

$L_{1r_0\left(u,w\right)}\left(r_0^{\text{**}},0,0\right){\left({\varphi }_{\text{**}},{\psi }_{\text{**}}\right)}^{\top }=\left(\begin{array}{c}\frac{{\varphi }_{\text{**}}}{1+k_0\alpha v_s^{d_2}}\\ 0\end{array}\right)$

Calculating the inner product of the above partial derivative operator at $\left({\varphi }_{\text{**}},{\psi }_{\text{**}}\right)$ with itself gives

$\langle L_{1r_0\left(u,w\right)}\left(r_0^{\text{**}},0,0\right)\left(\begin{array}{c}{\varphi }_{\text{**}}\\ {\psi }_{\text{**}}\end{array}\right),\left(\begin{array}{c}{\varphi }_{\text{**}}\\ {\psi }_{\text{**}}\end{array}\right)\rangle ={\displaystyle {\int }_{\Omega }\frac{{\varphi }_{\text{**}}^2}{1+k_0\alpha v_s^{d_2}}\text{d}x}>0.$

Since this inner product is strictly greater than 0, we conclude that

$L_{1r_0\left(u,w\right)}\left(r_0^{\text{**}},0,0\right)\left(\begin{array}{c}{\varphi }_{\text{**}}\\ {\psi }_{\text{**}}\end{array}\right)\notin \text{Range}\left(L_{1\left(u,w\right)}\left(r_0^{\text{**}},0,0\right)\right).$

Thus, the transversality condition is satisfied. Combining the results from the bifurcation theory in [19], we conclude that the point $\left(r_0^{\text{**}},0,v_s^{d_2}\right)$ is the unique bifurcation point from the semi-trivial solution branch ${\Theta }_v$ of model (2) producing a positive solution branch, and all positive solutions of model (2) in the neighborhood of this point lie on a smooth curve ${\theta }_2$ , with the parametric expression given by

${\theta }_2=\left\{\left(r_0\left(\theta \right),u\left(\theta \right),v\left(\theta \right)\right)=\left(r_0^{**}+\theta {{\widehat{r}}^{\prime }}_0\left(0\right)+o\left(\left|\theta \right|\right),\theta {\varphi }_{**}+o\left(\left|\theta \right|\right),v_s^{d_2}+\theta {\psi }_{**}+o\left(\left|\theta \right|\right)\right):0<\theta <{\theta }_2\right\},$

where ${\theta }_2>0$ is a constant.

Next, we derive the expression for ${\widehat{r}}_{0\text{'}}\left(0\right)$ . Define the operator

${\widehat{L}}_{1\left(u,w\right)}\left(r_0^{\text{**}},u,w\right)=L_{1\left(u,w\right)}\left(r_0^{\text{**}},u,w\right)\left(\begin{array}{l}\varphi \hfill \\ \psi \hfill \end{array}\right),$

then

${\widehat{L}}_{1\left(u,w\right)}\left(r_0^{\text{**}},u,w\right)\left(\begin{array}{l}\varphi \hfill \\ \psi \hfill \end{array}\right)=\left(\begin{array}{l}C\hfill \\ D\hfill \end{array}\right),$

where the components $C,D$ are given by

$\begin{array}{c}C=-2a{\varphi }^2+\frac{2pq\left(w+v_s^{d_2}\right)}{{\left(1+qu\right)}^3}{\varphi }^2+2\alpha \nabla \varphi \nabla \psi +2\alpha \varphi \Delta \psi \\ -\frac{2r_0{k}_0\alpha }{{\left(1+k_0\alpha \left(w+v_s^{d_2}\right)\right)}^2}\varphi \psi -\frac{2p}{{\left(1+qu\right)}^2}\varphi \psi +\frac{2r_0{k}_0^2{\alpha }^{2}u}{{\left(1+k_0\alpha \left(w+v_s^{d_2}\right)\right)}^2}{\psi }^2,\\ D=-\frac{m\left(x\right){\left(w+v_s^{d_2}\right)}^2}{{\left(r+u\right)}^3}{\varphi }^2-\frac{cpq\left(w+v_s^{d_2}\right)}{{\left(1+qu\right)}^3}{\varphi }^2+\frac{4m\left(x\right)\left(w+v_s^{d_2}\right)}{{\left(r+u\right)}^2}\varphi \psi \\ +\frac{2cp}{{\left(1+qu\right)}^2}\varphi \psi -\frac{2m\left(x\right)}{r+u}{\psi }^2.\end{array}$

According to relevant conclusions in [21], the expression for ${{\widehat{r}}^{\prime }}_0\left(0\right)$ can be derived from the following inner product relation

$\begin{array}{c}{{\widehat{r}}^{\prime }}_0\left(0\right)=-\frac{\langle {\widehat{L}}_{1\left(u,w\right)}\left(r_0^{**},0,0\right)\left(\begin{array}{c}{\varphi }_{**}\\ {\psi }_{**}\end{array}\right),l_1\rangle }{\langle 2L_{1r_0\left(u,w\right)}\left(r_0^{**},0,0\right)\left(\begin{array}{c}{\varphi }_{**}\\ {\psi }_{**}\end{array}\right),l_1\rangle }\\ ={\displaystyle {\int }_{\Omega }[\left(a-pq{v}_s^{d_2}\right){\varphi }_{**}^3-\left(\alpha \nabla {\varphi }_{**}\nabla {\psi }_{**}+\alpha {\varphi }_{**}\Delta {\psi }_{**}\right){\varphi }_{**}\underset{}{\overset{\stackrel{}{}}{}}}\\ +\frac{r_0{k}_0\alpha }{{\left(1+k_0\alpha v_s^{d_2}\right)}^2}{\varphi }_{**}^2{\psi }_{**}+p{\varphi }_{**}^2{\psi }_{**}]\text{d}x/{\displaystyle {\int }_{\Omega }\frac{{\varphi }_{**}^2}{1+k_0\alpha v_s^{d_2}}\text{d}x}.\end{array}$

Here, $l_1$ is a linear functional on the space $Y$ satisfying $\langle \left(\begin{array}{l}f\left(x\right)\hfill \\ g\left(x\right)\hfill \end{array}\right),l_1\rangle ={\displaystyle {\int }_{\Omega }f\left(x\right){\varphi }_{\text{**}}\text{d}x}$ . Thus, if $a>\tilde{a}$ and ${{\widehat{r}}^{\prime }}_0\left(0\right)>0$ , the solution branch ${\theta }_2$ undergoes a supercritical bifurcation at the equilibrium point $\left(r_0^{\text{**}},0,v_s^{d_2}\right)$ ; if $a<\tilde{a}$ and ${{\widehat{r}}^{\prime }}_0\left(0\right)<0$ , then this bifurcation is subcritical.

Theorem 3.2 The following conclusions hold

• If $-\frac{cp}{q}<s<0$ , there exists a small positive number ${\theta }_1$ such that the positive solutions $\left(r_0\left(\theta \right),u\left(\theta \right),v\left(\theta \right)\right)$ bifurcating from $\left(r_0^{\text{*}},\frac{r_0^{\text{*}}-d}{a},0\right)$ in model (2) are non-degenerate for $\theta \in \left(0,{\theta }_1\right)$ . Furthermore, if $p>\tilde{p}$ , then $\left(u\left(\theta \right),v\left(\theta \right)\right)$ is locally asymptotically stable; if $p<\tilde{p}$ , then $\left(u\left(\theta \right),v\left(\theta \right)\right)$ is unstable.

• If $0<s<{\lambda }_1^D\left(d_2,{\Omega }_0\right)$ , there exists a small positive number ${\theta }_2$ such that the positive solutions $\left(r_0\left(\theta \right),u\left(\theta \right),v\left(\theta \right)\right)$ bifurcating from $\left(r_0^{\text{**}},0,v_s^{d_2}\right)$ in model (2) are non-degenerate for $\theta \in \left(0,{\theta }_2\right)$ . Moreover, if $\underset{\theta \to 0^{+}}{\text{lim}}\frac{\nabla v\left(\theta \right)}{\theta }\ne \infty$ , then when $a<\tilde{a}$ , $\left(u\left(\theta \right),v\left(\theta \right)\right)$ is locally asymptotically stable; otherwise, it is unstable.

Proof. We provide a detailed proof for (ii) since (i) can be proven similarly.

Firstly, from the proof of Theorem 3.1, there exists a constant ${\theta }_2>0$ such that model (2) has a positive solution branch $\left(r_0\left(\theta \right),u\left(\theta \right),v\left(\theta \right)\right)$ bifurcating from $\left(r_0^{\text{**}},0,v_s^{d_2}\right)$ , where $\theta \in \left(0,{\theta }_2\right)$ . Linearizing the model (2) at the positive solution $\left(u,v\right)=\left(u\left(\theta \right),v\left(\theta \right)\right)$ and substituting $r_0=r_0\left(\theta \right)$ gives

$\chi \left(\theta \right)\left(\begin{array}{c}\varphi \left(\theta \right)\\ \psi \left(\theta \right)\end{array}\right)=\beta \left(\theta \right)\left(\begin{array}{c}\varphi \left(\theta \right)\\ \psi \left(\theta \right)\end{array}\right)\text{where}\chi \left(\theta \right)=\left(\begin{array}{ll}{\chi }_1\left(\theta \right)\hfill & {\chi }_2\left(\theta \right)\hfill \\ {\chi }_3\left(\theta \right)\hfill & {\chi }_4\left(\theta \right)\hfill \end{array}\right)$ (12)

$\begin{array}{l}{\chi }_1\left(\theta \right)=d_1\Delta +\alpha \nabla v\left(\theta \right)\nabla +\alpha \Delta v\left(\theta \right)+\frac{r_0\left(\theta \right)}{1+k_0\alpha v\left(\theta \right)}-d-2au\left(\theta \right)-\frac{pv\left(\theta \right)}{{\left(1+qu\left(\theta \right)\right)}^2},\\ {\chi }_2\left(\theta \right)=\alpha \nabla u\left(\theta \right)\nabla +\alpha u\left(\theta \right)\Delta -\frac{r_0\left(\theta \right)k_0\alpha u\left(\theta \right)}{{\left(1+k_0\alpha v\left(\theta \right)\right)}^2}-\frac{pu\left(\theta \right)}{1+qu\left(\theta \right)},\\ {\chi }_3\left(\theta \right)=\frac{m\left(x\right)v^2\left(\theta \right)}{{\left(r+u\left(\theta \right)\right)}^2}+\frac{cpv\left(\theta \right)}{{\left(1+qu\left(\theta \right)\right)}^2},\\ {\chi }_4\left(\theta \right)=d_2\Delta +s-\frac{2m\left(x\right)v\left(\theta \right)}{r+u\left(\theta \right)}+\frac{cpu\left(\theta \right)}{1+qu\left(\theta \right)}.\end{array}$

As $\left(r_0\left(\theta \right),u\left(\theta \right),v\left(\theta \right)\right)\to \left(r_0^{\text{**}},0,v_s^{d_2}\right)$ when $\theta \to 0^{+}$ , the operator $\chi \left(\theta \right)$ satisfies $\chi \left(\theta \right)\to {\chi }_0 \left(\theta \to 0^{+}\right)$ , where ${\chi }_0$ is given by

${\chi }_0=\left(\begin{array}{cc}{d}_1\Delta +\alpha \nabla v_s^{d_2}\nabla +\alpha \Delta v_s^{d_2}+\frac{r_0\left(\theta \right)}{1+k_0\alpha v_s^{d_2}}-d-p{v}_s^{d_2}& 0\\ \frac{m\left(x\right)v_s^{d_2}}{r^2}+cp{v}_s^{d_2}& d_2\Delta +s-\frac{2m\left(x\right)v_s^{d_2}}{r}\end{array}\right)$

It is not difficult to verify that zero is the smallest eigenvalue of the operator ${\chi }_0$ , with the corresponding positive eigenfunction being $\left({\varphi }_{\text{**}},{\psi }_{\text{**}}\right)$ , which has been defined in the proof of Theorem 3.1 (ii). Moreover, all other eigenvalues of the operator ${\chi }_0$ have positive real parts that are uniformly away from zero. Therefore, by combining linear operator perturbation theory [22], we can conclude that when $\theta >0$ is sufficiently small, the operator $\chi \left(\theta \right)$ has a unique principal eigenvalue $\beta \left(\theta \right)$ and corresponding eigenfunction $\left(\varphi \left(\theta \right),\psi \left(\theta \right)\right)$ , satisfying

$\beta \left(\theta \right)\to 0,\left(\varphi \left(\theta \right),\psi \left(\theta \right)\right)\to \left({\varphi }_{\text{**}},{\psi }_{\text{**}}\right),\theta \to 0^{+}.$

Furthermore, all other eigenvalues of $\chi \left(\theta \right)$ also have positive real parts and are uniformly away from zero. Thus, there exists a sufficiently small positive number ${\theta }_2$ such that the positive solution branch $\left(r_0\left(\theta \right),u\left(\theta \right),v\left(\theta \right)\right)$ bifurcating from the trivial steady state $\left(r_0^{\text{**}},0,v_s^{d_2}\right)$ is non-degenerate in the interval $\theta \in \left(0,{\theta }_2\right)$ .

Now, we study the stability of $\left(u\left(\theta \right),v\left(\theta \right)\right)$ . From (12), we find that $\varphi \left(\theta \right)$ satisfies

$\begin{array}{c}\beta \left(\theta \right)\varphi \left(\theta \right)=d_1\Delta \varphi \left(\theta \right)+\alpha \nabla v\left(\theta \right)\nabla \varphi \left(\theta \right)+\alpha \varphi \left(\theta \right)\Delta v\left(\theta \right)+\frac{r_0\left(\theta \right)\varphi \left(\theta \right)}{1+k_0\alpha v\left(\theta \right)}\\ -d\varphi \left(\theta \right)-2au\left(\theta \right)\varphi \left(\theta \right)-\frac{pv\left(\theta \right)}{{\left(1+qu\left(\theta \right)\right)}^2}\varphi \left(\theta \right)+\alpha \nabla u\left(\theta \right)\nabla \psi \left(\theta \right)\\ +\alpha u\left(\theta \right)\Delta \psi \left(\theta \right)-\frac{r_0\left(\theta \right)k_0\alpha u\left(\theta \right)}{{\left(1+k_0\alpha v\left(\theta \right)\right)}^2}\psi \left(\theta \right)-\frac{pu\left(\theta \right)}{1+qu\left(\theta \right)}\psi \left(\theta \right),\end{array}$ (13)

Multiplying both sides of (13) by $u\left(\theta \right)$ and integrating over Ω, we obtain the equation corresponding to $u\left(\theta \right)$

$\begin{array}{c}\beta \left(\theta \right){\displaystyle {\int }_{\Omega }\varphi \left(\theta \right)u\left(\theta \right) \text{d}x}=\alpha {\displaystyle {\int }_{\Omega }(\varphi \left(\theta \right)\nabla u\left(\theta \right)\nabla v\left(\theta \right)-u\left(\theta \right)\nabla \varphi \left(\theta \right)\nabla v\left(\theta \right)}\\ -u\left(\theta \right)\varphi \left(\theta \right)\Delta v\left(\theta \right)-u\left(\theta \right)\nabla u\left(\theta \right)\nabla \psi \left(\theta \right)\\ -u^2\left(\theta \right)\Delta \psi \left(\theta \right))\text{d}x\\ +{\displaystyle {\int }_{\Omega }\left(a{u}^2\left(\theta \right)\varphi \left(\theta \right)-\frac{pq{u}^2\left(\theta \right)v\left(\theta \right)\varphi \left(\theta \right)}{{\left(1+qu\left(\theta \right)\right)}^2}\right)\text{d}x}\\ +{\displaystyle {\int }_{\Omega }\left(\frac{r_0\left(\theta \right)k_0\alpha u\left(\theta \right)}{{\left(1+k_0\alpha v\left(\theta \right)\right)}^2}+\frac{p{u}^2\left(\theta \right)\psi \left(\theta \right)}{1+qu\left(\theta \right)}\right)\text{d}x},\end{array}$ (14)

Dividing both sides of (14) by ${\theta }^2$ and letting $\theta \to 0$ , we obtain

$\underset{\theta \to 0^{+}}{\text{lim}}\frac{\beta \left(\theta \right)}{\theta }=\frac{N}{{\displaystyle {\int }_{\Omega }{\varphi }_{\text{**}}^2\text{d}x}}$

where

$N={\displaystyle {\int }_{\Omega }\left(a{\varphi }_{**}^3-pq{v}_s^{d_2}{\varphi }_{**}^3-\alpha {\varphi }_{**}^2\Delta {\psi }_{**}-\alpha {\varphi }_{**}\nabla {\varphi }_{**}\nabla {\psi }_{**}+\left(\frac{r_0^{**}{k}_0\alpha }{{\left(1+k_0\alpha v_s^{d_2}\right)}^2}+p\right){\varphi }_{**}^2{\psi }_{**}\right)\text{d}x}.$

Therefore, when $a>\tilde{a}$ , we have $N>0$ , which implies that ${\beta }^{\prime }\left(0\right)>0$ , indicating that the bifurcated solution is unstable. Conversely, when $a<\tilde{a}$ , we find $N<0$ , leading to ${\beta }^{\prime }\left(0\right)<0$ , which signifies that the bifurcated solution is locally stable.

4. Discussion

This paper focuses on a class of Holling-Tanner predator-prey models that incorporate predator attraction behavior and fear effects. Using the prey population growth rate $r_0$ as a bifurcation parameter, we theoretically derive the necessary and sufficient conditions for the local bifurcation of the model from semi-trivial solutions and analyze the linear stability of the local bifurcated solutions. The results indicate that the stability of the bifurcated solutions depends on the relationship between internal system parameters (such as the intraspecific competition coefficient of prey $a$ , the intensity of predator interference $p$ , etc.) and their corresponding critical threshold values: when the parameter values lead to a supercritical bifurcation, small amplitude positive steady-state solutions are linearly stable, suggesting that stable spatially heterogeneous patterns can form in the vicinity of the critical parameters. Conversely, when the bifurcation is subcritical, the corresponding positive steady-state solutions are linearly unstable, indicating that the system may exhibit delayed phenomena or branching behavior of large amplitude solutions. These theoretical results not only enhance the qualitative theoretical framework of predator-prey models incorporating fear effects and chemotaxis in spatially heterogeneous environments but also provide a solid mathematical foundation for understanding the spatial distribution patterns of populations driven by fear effects in actual ecological systems.

Future research will focus on the following directions: 1) in-depth analysis of the global bifurcation structure; 2) bifurcation analysis under the combined effects of multiple parameters; 3) asymptotic behavior in degenerate critical cases; 4) numerical simulations to validate theoretical results.