Dynamics of an unstirred chemostat model with Beddington–DeAngelis functional response

This paper deals with an unstirred competitive chemostat model with the Beddington–DeAngelis functional response. With the help of the linear eigenvalue theory and the monotone dynamical system theory, we establish a relatively clear dynamic classification of this system in terms of the growth rates of two species. The results indicate that there exist several critical curves, which may classify the dynamics of this system into three scenarios: 1) extinction; 2) competitive exclusion; and 3) coexistence. Comparing with the classical chemostat model [26], our theoretical results reveal that under the weak–strong competition cases, the role of intraspecific competition can lead to species coexistence. Moreover, the simulations suggest that under different competitive cases, coexistence can occur for suitably small diffusion rates and some intermediate diffusion rates. These new phenomena indicate that the intraspecific competition and diffusion have a great influence on the dynamics of the unstirred chemostat model of two species competing with the Beddington–DeAngelis functional response.

1 Introduction

It is well known that the chemostat is a laboratory apparatus used for the continuous culture of microorganisms, while the chemostat models are extensively applied in ecology to simulate the growth of single-celled algal plankton in oceans and lakes [1–4]. Most of the earlier chemostat models assume the well-stirring of culture, which leads to chemostat models generally described by ordinary differential equation models (see, e.g., [2, 4, 5]). However, this idealized mixing is quite different from the real environment in which microbial populations live. Since the ability of microorganisms to move in a random fashion plays an important role in determining the survival and extinction of populations, many unstirred chemostat models have sprung up in which populations and resources are distributed in spatially variable habitats; please refer to [6–10] for small sampling of such works.

There are various types of response functions; among them, Holling types I–IV [11] are usually introduced to model the growth of microorganisms. Particularly, the various chemostat models with Holling type II functional response have been extensively studied (see, e.g., [4, 10, 12, 13]). As far as we know, for the unstirred competitive chemostat models with Holling type II functional response, So and Waltman [14] first obtained the local coexistence by standard bifurcation theorems. Later, Hsu and Waltman [6] obtained the asymptotic behavior of solutions by the theory of uniform persistence in an infinite-dimensional dynamical system and the theory of strongly order-preserving semi-dynamical system. To explore the effect of diffusion, Shi et al. [8] further studied this model and confirmed that stable coexistence solutions only occur at the intermediate diffusion rates. In addition, a diffusive predator–prey chemostat model with Holling type II functional response was studied by Nie et al. [7], and their analytical and numerical results show that a relatively small diffusion is conducive to the coexistence of species.

However, in nature, it is known that there is not only competition between two species but also mutual interference in species. Therefore, it is necessary to consider mutual interference in species. To this end, Beddington [15] and DeAngelis et al. [16] (simplified as B.–D.) proposed the following B.–D. functional response:

$f_1\left(S,u\right)=\frac{S}{k_1+S+{\beta }_{1}u},f_2\left(S,v\right)=\frac{S}{k_2+S+{\beta }_{2}v},$(1.1)

where k_i > 0 (i = 1, 2) are the Michaelis–Menten constants, S represents the density of the resources, u and v represent the density of two species, respectively, and β_i > 0 (i = 1, 2) model the mutual interference between two species.

As illustrated by Harrision [17], the B.–D. functional response with intraspecific interference competition was superior to well-known Holling type II functional response in modeling the resource uptake of species. Therefore, there appear successively many works to describe the population dynamics by using the B.–D. functional response. For instance, Jiang et al. [18] discussed a competition model with the B.–D. functional response, and they applied the fixed-point index theory to obtain the sufficient conditions for the existence of positive solutions. In addition, a predator–prey model with a heterogeneous environment and the B.–D. functional response was constructed by Zhang and Wang [19], and the existence of positive stationary solutions was obtained by using the fixed-point index theory. We also refer the recent works [20–22] about population models with the B.–D. functional response.

Particularly, the unstirred chemostat models with the B.–D. functional response have also received considerable attention in the past decades. Wang et al. [23] obtained the sufficient conditions for the existence of positive steady-state solutions and studied the effect of parameter β₁ on coexistence states by the fixed-point index theory, the perturbation technique, and the bifurcation theory. Meanwhile, Nie and Wu [24] studied the unstirred chemostat model with the B.–D. functional response and inhibitor, and the uniqueness, multiplicity, and stability of the coexistence solutions were obtained by the degree theory in cones, bifurcation theory, and perturbation technique. More works on chemostat models with the B.–D. functional response can be found in [25–27] and the references therein.

Mathematically speaking, these sufficient conditions for the existence of coexistence solutions are usually established in terms of the principal eigenvalues of the corresponding linearized eigenvalue problems at trivial or semi-trivial steady states (see, e.g., [18, 19, 23, 24]). It is worth noting that these principal eigenvalues depend heavily on the model parameters, which motivates us to explore how these model parameters affect the existence of coexistence solutions and establish the dynamics classification of this system in terms of these model parameters. Moreover, it should be noted that studying the asymptotic analysis of steady states of chemostat models is non-trivial, and some new techniques need to be introduced. Overall, for the unstirred chemostat system with the B.–D. functional response, we are concerned with the following questions:

(1) How do parameters such as diffusion rates, growth rates, and intraspecific competition parameters affect the dynamics of the unstirred chemostat system with the B.–D. functional response?

(2) Can we establish a clear dynamic classification of the unstirred chemostat system with the B.–D. functional response in terms of these parameters?

(3) Will there arise a new phenomenon if one introduces the B.–D. functional response into the unstirred chemostat model?

The purpose of this paper is to address these problems. We hope that the approaches in this paper might provide some new insights on the dynamical behavior of the unstirred chemostat models.

This paper is organized as follows. In Section 2, we introduce an unstirred chemostat model with the B.–D. functional response and its corresponding limiting system. In Section 3, some preliminary results are given. In Section 4, we aim to investigate the dynamics of this limiting system and obtain a relatively clear dynamic classification of this limiting system in the m₁ − m₂ plane. In Section 5, the coexistence solution for this limiting system is established by a bifurcation argument. In Section 6, we study the effect of diffusion on system dynamics by numerical approaches. In Section 7, a discussion is presented from the opinion of analytic and numerical results. Finally, the proofs of some theoretical results are deferred to the Supplementary Appendix in Supplementary Section S8.

2 The model

In this paper, we consider following the unstirred chemostat system with the B.–D. functional response:

$\left\{\begin{array}{cc}{S}_t=d{S}_{xx}-m_{1}u{f}_1\left(S,u\right)-m_{2}v{f}_2\left(S,v\right),\hfill & x\in \left(\mathrm{0,1}\right),t>0,\hfill \\ u_t=d{u}_{xx}+m_{1}u{f}_1\left(S,u\right),\hfill & x\in \left(\mathrm{0,1}\right),t>0,\hfill \\ v_t=d{v}_{xx}+m_{2}v{f}_2\left(S,v\right),\hfill & x\in \left(\mathrm{0,1}\right),t>0,\hfill \\ S_x\left(0,t\right)=-S^0,S_x\left(1,t\right)+\gamma S\left(1,t\right)=0,\hfill & t>0,\hfill \\ u_x\left(0,t\right)=u_x\left(1,t\right)+\gamma u\left(1,t\right)=0,\hfill & t>0,\hfill \\ v_x\left(0,t\right)=v_x\left(1,t\right)+\gamma v\left(1,t\right)=0,\hfill & t>0,\hfill \\ S\left(x,0\right)=S_0\left(x\right)\ge 0,\hfill & x\in \left[\mathrm{0,1}\right],\hfill \\ u\left(x,0\right)=u_0\left(x\right)\ge ,\not\equiv 0,v\left(x,0\right)=v_0\left(x\right)\ge ,\not\equiv 0,\hfill & x\in \left[\mathrm{0,1}\right],\hfill \end{array}\right.$(2.1)

where S(x, t) is the concentration of the nutrient and u(x, t) and v(x, t) represent the population density for the two competing microorganisms with location x and time t, respectively. The positive constants m₁ and m₂ are corresponding to the growth rates of species u and v with nutrient concentration S. d > 0 is the diffusion rate of the nutrient and microorganisms. The initial data S₀(x), u₀(x), and v₀(x) are non-negative non-trivial continuous functions. In the reactor, the nutrients are pumped with the rate of S⁰ > 0 at position x = 0, and the mixed cultures containing nutrients and microorganisms are pumped out with the rate of γ > 0 at the position x = 1, which results in the Robin boundary conditions at x = 1 [6]. Here, f₁(S, u), f₂(S, v) satisfying Eq. 1.1 are the nutrient uptake of species u and v at nutrient concentration S. Moreover, we redefine f₁(S, u), f₂(S, v) as follows [10]:

${\widehat{f}}_1\left(S,u\right)=\left\{\begin{array}{cc}{f}_1\left(S,u\right),\hfill & S\ge 0,u\ge 0,\hfill \\ 0,\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{s},\hfill \end{array}\right.{\widehat{f}}_2\left(S,v\right)=\left\{\begin{array}{cc}{f}_2\left(S,v\right),\hfill & S\ge 0,v\ge 0,\hfill \\ 0,\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{s}.\hfill \end{array}\right.$

For convenience, we still denote ${\widehat{f}}_i(S,u)(i=\mathrm{1,2})$ as f_i(S, u), throughout this paper.

It is worth pointing out that system (2.1) satisfies the conservation law [4]. In other words, the total biomass concentration S + u + v in the chemostat approaches asymptotically a steady state $\varphi (x)=S^0(\frac{1+\gamma }{\gamma }-x)$ (see [14], Lemma 2.1); that is,

$\underset{t\to \infty }{\lim }\left(S\left(x,t\right)+u\left(x,t\right)+v\left(x,t\right)\right)=\varphi \left(x\right)\mathrm{u}\mathrm{n}\mathrm{i}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{l}\mathrm{y}\mathrm{f}\mathrm{o}\mathrm{r}x\in \left[\mathrm{0,1}\right].$

Hence, we apply the classical internal chain transitive theory [[28], Lemma 2.1] to reduce system (2.1) into the following limiting system:

$\left\{\begin{array}{cc}{u}_t=d{u}_{xx}+m_1{f}_1\left(\varphi \left(x\right)-u-v,u\right)u,\hfill & x\in \left(\mathrm{0,1}\right),t>0,\hfill \\ v_t=d{v}_{xx}+m_2{f}_2\left(\varphi \left(x\right)-u-v,v\right)v,\hfill & x\in \left(\mathrm{0,1}\right),t>0,\hfill \\ u_x\left(0,t\right)=u_x\left(1,t\right)+\gamma u\left(1,t\right)=0,v_x\left(0,t\right)=v_x\left(1,t\right)+\gamma v\left(1,t\right)=0,\hfill & t>0,\hfill \\ u\left(x,0\right)=u_0\left(x\right)\ge ,\not\equiv 0,v\left(x,0\right)=v_0\left(x\right)\ge ,\not\equiv 0,u_0\left(x\right)+v_0\left(x\right)\le \varphi \left(x\right)\hfill & x\in \left[\mathrm{0,1}\right].\hfill \end{array}\right.$(2.2)

In this paper, we are mainly concerned with the dynamics classification of system (2.2). Based on the competition relationship of two species, system (2.2) generates a strictly monotone dynamical system in the partial competitive order induced by the cone K = {(u, v) ∈ C[0, 1] × C[0, 1]: u ≥ 0, v ≤ 0} (see [27], Proposition 1.3 in Chapter 8). Since the dynamics of Eq. 2.2 are related to the stability of non-negative steady states [29], we also focus on the following steady-state system:

$\left\{\begin{array}{cc}d{u}_{xx}+m_1{f}_1\left(\varphi \left(x\right)-u-v,u\right)u=0,\hfill & x\in \left(\mathrm{0,1}\right),\hfill \\ d{v}_{xx}+m_2{f}_2\left(\varphi \left(x\right)-u-v,v\right)v=0,\hfill & x\in \left(\mathrm{0,1}\right),\hfill \\ u_x\left(0\right)=u_x\left(1\right)+\gamma u\left(1\right)=0,v_x\left(0\right)=v_x\left(1\right)+\gamma v\left(1\right)=0.\hfill \end{array}\right.$(2.3)

The contribution of this paper is to explore the effect of these model parameters on the dynamics of system (2.2). Precisely, we first apply the linear eigenvalue theory and the monotone dynamical system theory to establish the threshold dynamics of system (2.2) in terms of growth rates and intraspecific competition parameters (see Theorems 4.1, 4.2). Moreover, we give a relatively clear dynamic classification of system (2.2) in the m₁ − m₂ plane (Figure 2). Finally, by numerical simulations, we further investigate the effect of diffusion on the dynamics of system (2.2) (Figures 3–5). Particularly, the numerical results show that under the different competitive cases, coexistence occurs for suitably small diffusion rates and some intermediate diffusion rates, which reveals that the dynamics of system (2.2) are relatively complicated.

3 Preliminaries

In this section, some preliminary results are presented, which are helpful in the following analysis.

We first consider a linear eigenvalue problem

$\left\{\begin{array}{cc}d{\phi }_{xx}+q\left(x\right)\phi =\mu \phi ,\hfill & x\in \left(\mathrm{0,1}\right),\hfill \\ {\phi }_x\left(0\right)={\varphi }_x\left(1\right)+\gamma \phi \left(1\right)=0,\hfill \end{array}\right.$(3.1)

where γ is a positive constant and q(x) ∈ C[0, 1]. For fixed d > 0, it is well known that problem (3.1) admits a principal eigenvalue μ₁(q(x)) [29], which corresponds to a positive eigenfunction φ₁(⋅, q(x)) normalized by $\underset{x\in [\mathrm{0,1}]}{\max }{\phi }_1=1$. Furthermore, by the variation characterization of the principal eigenvalue [29], we have

${\mu }_1\left(q\left(x\right)\right)=-\underset{\phi \in H^1\left(\mathrm{0,1}\right),\phi \ne 0}{\inf }\frac{d{\int }_0^1{\left({\phi }_x\right)}^2\mathrm{d}x+d\gamma {\phi }^2\left(1\right)-{\int }_0^{1}q\left(x\right){\phi }^2\mathrm{d}x}{{\int }_0^1{\phi }^2\mathrm{d}x}.$(3.2)

Moreover, the principal eigenvalue μ₁(q(x)) has the following properties.

Lemma 3.1 (See [21], Lemma 2.1). The following statements on the principal eigenvalue μ₁(q(x)) are true:

(i) μ₁(q(x)) depends continuously and differentially on parameter d in (0, + ∞), and it is strictly decreasing with respect to d in (0, + ∞).

(ii) q_n(x) → q(x) in C[0, 1] implies μ₁(q_n(x)) → μ₁(q(x)).

(iii) q₁(x) ≥ q₂(x) implies that μ₁(q₁(x)) ≥ μ₁(q₂(x)), and the equality holds only if q₁(x) ≡ q₂(x). Particularly, μ₁(0) < 0.

We consider the following single-species model:

$\left\{\begin{array}{cc}{\omega }_t=d{\omega }_{xx}+mf\left(\varphi -\omega ,\omega \right)\omega ,\hfill & x\in \left(\mathrm{0,1}\right),t>0,\hfill \\ {\omega }_x\left(0,t\right)={\omega }_x\left(1,t\right)+\gamma \omega \left(1,t\right)=0,\hfill & t>0,\hfill \\ \omega \left(x,0\right)={\omega }_0\left(x\right)\ge ,\not\equiv 0,\hfill & x\in \left[\mathrm{0,1}\right],\hfill \end{array}\right.$(3.3)

where d, m > 0 are constants and $f(\varphi -\omega ,\omega )=\frac{\varphi -\omega }{k+\varphi +(\beta -1)\omega }$. For fixed d, k > 0, let μ₁(mf(ϕ, 0)) be the principal eigenvalue of

$\left\{\begin{array}{cc}d{\phi }_{xx}+mf\left(\varphi ,0\right)\phi =\mu \phi ,\hfill & x\in \left(\mathrm{0,1}\right),t>0,\hfill \\ {\phi }_x\left(0\right)={\phi }_x\left(1\right)+\gamma \phi \left(1\right)=0,\hfill \end{array}\right.$

where $f(\varphi ,0)=\frac{\varphi }{k+\varphi }$. Then, we can conclude from Lemma 3.1(iii) that μ₁(mf(ϕ, 0)) is strictly increasing with respect to m in (0, + ∞). Moreover, $\underset{m\to 0^{+}}{\lim }{\mu }_1(mf(\varphi ,0))={\mu }_1(0)<0$ by Lemma 3.1(iii), and $\underset{m\to +\infty }{\lim }{\mu }_1(mf(\varphi ,0))=+\infty$ by Eq. 3.2. For fixed d, k > 0, there exists a unique critical value m^⋆ such that

$\left\{\begin{array}{c}{\mu }_1\left(mf\left(\varphi ,0\right)\right)<0\mathrm{i}\mathrm{f}0<m<m^{\ast },\hfill \\ {\mu }_1\left(mf\left(\varphi ,0\right)\right)=0\mathrm{i}\mathrm{f}m=m^{\ast },\hfill \\ {\mu }_1\left(mf\left(\varphi ,0\right)\right)>0\mathrm{i}\mathrm{f}m>m^{\ast }.\hfill \end{array}\right.$(3.4)

To stress the dependence of the unique positive steady state of system (3.3) on m and β, let us denote it by ω_∗(⋅; m, β).

Lemma 3.2 Suppose d, m, β, k > 0. Let ω(x, t) be the solution of system (3.3). Then,

(i) if m > m^∗, system (3.3) admits a unique positive steady state 0 < ω_∗(⋅; m, β) < ϕ(x) for x ∈ [0, 1], and $\underset{t\to \infty }{\lim }\omega (x,t)={\omega }_{\ast }(\cdot ;m,\beta )$ uniformly on [0,1];

(ii) if m ≤ m^∗, system (3.3) has no positive steady state and $\underset{t\to \infty }{\lim }\omega (x,t)=0$ uniformly on [0,1].

The proof of Lemma 3.2 is similar to the arguments in [6], Theorem 3.2. So, we omit it here.

We next give some asymptotic properties of the unique positive steady state ω_∗(⋅; m, β) of system (3.3) by taking m and β as the variable parameters.

Lemma 3.3 suppose that m > m^∗ holds. The following statements about the positive solution ω_∗(; m, β) will hold.

(i) For fixed d, k, β > 0, there exists positive solution ω_∗(; m, β), which is continuously differentiable with respect to m in (m^∗, + ∞), and it is point-wise strictly increasing in m ∈ (m^∗, + ∞). Moreover,

$\underset{m\to {\left(m^{\ast }\right)}^{+}}{\lim }{\omega }_{\ast }\left(\cdot ;m,\beta \right)=0,\underset{m\to +\infty }{\lim }{\omega }_{\ast }\left(\cdot ;m,\beta \right)=\varphi \left(x\right)\text{uniformly}\text{on}\left[\mathrm{0,1}\right].$(3.5)

(ii) For fixed d, k > 0 and m > m^∗, there exists positive solution ω_∗(⋅; m, β), which is continuously differentiable with respect to β in (0, + ∞), and it is point-wise strictly decreasing in β ∈ (0, + ∞). Moreover,

$\underset{\beta \to +\infty }{\lim }{\omega }_{\ast }\left(\cdot ;m,\beta \right)=0\text{uniformly}\text{on}\left[\mathrm{0,1}\right].$(3.6)

Proof For (i), it follows from Lemma 3.2 that ω_∗(⋅; m, β) exists if and only if m > m^∗. Moreover, ω_∗(⋅; m, β) is continuously differentiable with respect to m in (m^∗, + ∞) refering to the arguments in [30], Lemma 5.4(ii). Differentiating the equation of ω_∗(⋅; m, β) with respect to m and denoting P_m(x) = ∂ω_∗(⋅; m, β)/∂m, P_m(x) satisfies

$\left\{\begin{array}{c}d{P}_m^{″}+m\left[f\left(\varphi -{\omega }_{\ast },{\omega }_{\ast }\right)-\frac{\left(k+\beta \varphi \right){\omega }_{\ast }}{{\left[k+\varphi +\left(\beta -1\right){\omega }_{\ast }\right]}^2}\right]\hfill \\ P_m=-f\left(\varphi -{\omega }_{\ast },{\omega }_{\ast }\right){\omega }_{\ast },x\in \left(\mathrm{0,1}\right),\hfill \\ P_m^{\prime }\left(0\right)=P_m^{\prime }\left(1\right)+\gamma P_m\left(1\right)=0.\hfill \end{array}\right.$(3.7)

We define $\mathbb{X}=\{\psi \in C^2[\mathrm{0,1}]:{\psi }^{\prime }(0)={\psi }^{\prime }(1)+\gamma \psi (1)=0\}$ and denote $L_1(\psi )=d{\psi }^{\prime \prime }+m[f(\varphi -{\omega }_{\ast },{\omega }_{\ast })-\frac{(k+\beta \varphi ){\omega }_{\ast }}{{[k+\varphi +(\beta -1){\omega }_{\ast }]}^2}]\psi$. It is easy to see that μ₁(L₁) < μ₁(mf(ϕ − ω_∗, ω_∗)) = 0 by Lemma 3.1(iii). Noting that $P_m(x)\in \mathbb{X}$ and L₁(P_m) = −f(ϕ − ω_∗, ω_∗)ω_∗ < 0, we have P_m(x) > 0 on [0,1] by the generalized maximum principle [23, Theorem 5], which implies that ω_∗(⋅; m, β) is point-wise strictly increasing in m ∈ (m^∗, + ∞).

Since 0 < ω_∗(⋅; m, β) < ϕ and ${\omega }_{\ast }^{\prime \prime }=-\frac{m}{d}f(\varphi -{\omega }_{\ast },{\omega }_{\ast }){\omega }_{\ast }<0$, ${\omega }_{\ast }^{\prime }$ is decreasing on [0,1]. Note that the boundary conditions ${\omega }_{\ast }^{\prime }(0)=0,{\omega }_{\ast }^{\prime }(1)=-\gamma {\omega }_{\ast }(1)$. Then, ${\omega }_{\ast }^{\prime }(x)$ is uniformly bounded for x ∈ [0, 1]. It follows from the Arzela–Ascoli theorem that there exist ω₁, ω₂ ∈ C[0, 1] with 0 ≤ ω₁ ≤ ϕ and 0 ≤ ω₂ ≤ ϕ such that

$\underset{m\to {\left(m^{\ast }\right)}^{+}}{\lim }{\omega }_{\ast }\left(\cdot ;m,\beta \right)={\omega }_1,\underset{m\to +\infty }{\lim }{\omega }_{\ast }\left(\cdot ;m,\beta \right)={\omega }_2\mathrm{i}\mathrm{n}C\left[\mathrm{0,1}\right].$

To prove ω₁ = 0 on [0,1], we assume by contradiction that ω₁≢0 on [0,1]. Since 0 < f(ϕ − ω_∗, ω_∗) < 1, the standard L^p-estimate implies that ω_∗(⋅; m, β) is uniformly bounded in W^2,p(0, 1) with p ∈ (1, ∞) for m ∈ (m^∗, M], where M is a fixed constant larger than m^∗. Therefore, $\underset{m\to {(m^{\ast })}^{+}}{\lim }{\omega }_{\ast }(\cdot ;m,\beta )={\omega }_1$ weakly in W², p (0, 1) and the convergence also holds in C¹[0, 1] by the Sobolev embedding theorem. Then, ω₁ satisfies

$\left\{\begin{array}{cc}d{\omega }_1^{\prime \prime }+m^{\ast }f\left(\varphi -{\omega }_1,{\omega }_1\right){\omega }_1=0,\hfill & x\in \left(\mathrm{0,1}\right),\hfill \\ {\omega }_1^{\prime }\left(0\right)={\omega }_1^{\prime }\left(1\right)+\gamma {\omega }_1\left(1\right)=0.\hfill \end{array}\right.$(3.8)

Since ω₁≢0 on [0,1], we have ω₁ > 0 on [0,1] by the strong maximum principle. It is easy to see that μ₁(m^∗f(ϕ, 0)) > μ₁(m^∗f(ϕ − ω₁, ω₁)) = 0 by Lemma 3.1(iii), a contradiction to the definition of m^∗ (Eq. 3.4). Thus, ω₁ = 0.

We next prove ω₂ = ϕ(x) on [0,1]. We recall that ω_∗(⋅; m, β) satisfies

$\left\{\begin{array}{cc}d{\omega }_{\ast }^{\prime \prime }+mf\left(\varphi -{\omega }_{\ast },{\omega }_{\ast }\right){\omega }_{\ast }=0,\hfill & x\in \left(\mathrm{0,1}\right),\hfill \\ {\omega }_{\ast }^{\prime }\left(0\right)={\omega }_{\ast }^{\prime }\left(1\right)+\gamma {\omega }_{\ast }\left(1\right)=0.\hfill \end{array}\right.$(3.9)

Dividing the first equation of Eq. 3.9 by mω_∗ and integrating over (0,1),

$\frac{d}{m}{\int }_0^1\frac{|{\omega }_{\ast }^{\prime }{|}^2}{{\omega }_{\ast }^2}\mathrm{d}x-\frac{d\gamma }{m}+{\int }_0^{1}f\left(\varphi -{\omega }_{\ast },{\omega }_{\ast }\right)\mathrm{d}x=0,$

which implies

$0\le {\int }_0^{1}f\left(\varphi -{\omega }_{\ast },{\omega }_{\ast }\right)\mathrm{d}x\le \frac{d\gamma }{m}.$

Note $\underset{m\to +\infty }{\lim }{\omega }_{\ast }(\cdot ;m,\beta )={\omega }_2$ in C[0, 1]. Taking m → +∞, we have ${\int }_0^{1}f(\varphi -{\omega }_2,{\omega }_2)\mathrm{d}x=0$, which means ω₂ = ϕ(x) on [0,1] by 0 ≤ ω₂ ≤ ϕ.

(ii) The monotonicity of ω_∗(⋅; m, β) with respect to β in (0, + ∞) can be proved by the similar arguments as in the proof of (i) and $\underset{\beta \to +\infty }{\lim }{\omega }_{\ast }(\cdot ;m,\beta )=0\mathrm{u}\mathrm{n}\mathrm{i}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{l}\mathrm{y}\mathrm{o}\mathrm{n}[\mathrm{0,1}]$ holds (see [23], Remark 1.2).

It is clear that system (2.2) generates a monotone dynamical system in the partial competitive order induced by the cone K = {(u, v) ∈ C[0, 1] × C[0, 1]: u ≥ 0, v ≤ 0} (see [27], Proposition 1.3 in Chapter 8). Hence, we can recall the well-known results on the monotone dynamical system as follows.

Lemma 3.4 [9]. For the monotone dynamical system,

(i) if two semi-trivial steady states are asymptotically stable, then it has at least one unstable coexistence steady state.

(ii) if two semi-trivial steady states are unstable, then it has at least one stable coexistence steady state. Furthermore, if its coexistence steady states are all linearly stable, then there is a unique coexistence steady state that is globally asymptotically stable.

(iii) if there is no coexistence steady state and if one semi-trivial solution is linearly unstable, the other semi-trivial solution is globally asymptotically stable.

4 The dynamics analysis of system (2.2)

As we already know, the local dynamics of system (2.2) are related to the stability of semi-trivial solutions [29]. Hence, we next establish the stability of semi-trivial solutions, including local and some global stability results. Recalling $f_i(\varphi ,0)=\frac{\varphi }{k_i+\varphi }(i=\mathrm{1,2})$, by the similar arguments as in (3.4), we can define $m_i^{\ast }$ such that

$\left\{\begin{array}{c}{\mu }_1\left(m_i{f}_i\left(\varphi ,0\right)\right)<0\mathrm{i}\mathrm{f}0<m_i<m_i^{\ast },\hfill \\ {\mu }_1\left(m_i{f}_i\left(\varphi ,0\right)\right)=0\mathrm{i}\mathrm{f}{m}_i=m_i^{\ast },\hfill \\ {\mu }_1\left(m_i{f}_i\left(\varphi ,0\right)\right)>0\mathrm{i}\mathrm{f}{m}_i>m_i^{\ast }.\hfill \end{array}\right.$(4.1)

Clearly, $m_i^{\ast }(i=\mathrm{1,2})$ are dependent on k_i but not on β_i.

Proposition 4.1 For fixed d > 0, the following statements hold:

$\left(\text{i}\right)\frac{k_2}{k_1}<\frac{m_2^{\ast }}{m_1^{\ast }}<1\mathit{\text{if}}{k}_1>k_2>0;\left(\text{ii}\right)1<\frac{m_2^{\ast }}{m_1^{\ast }}<\frac{k_2}{k_1}\mathit{\text{if}}{k}_2>k_1>0;\left(\text{iii}\right)m_1^{\ast }=m_2^{\ast }\mathit{\text{if}}{k}_1=k_2>0.$

Proof (i) Since

${\mu }_1\left(\frac{m_1^{\ast }\varphi }{k_1+\varphi }\right)=0\mathrm{a}\mathrm{n}\mathrm{d}{\mu }_1\left(\frac{m_2^{\ast }\varphi }{k_2+\varphi }\right)=0,$(4.2)

it is easy to check that $m_1^{\ast }>m_2^{\ast }$ by k₁ > k₂ > 0 and Lemma 3.1(iii). Note that Eq. 4.2) is equivalent to ${\mu }_1(\frac{m_1^{\ast }\varphi /k_1}{1+\varphi /k_1})=0$ and ${\mu }_1(\frac{m_2^{\ast }\varphi /k_2}{1+\varphi /k_2})=0$. Since $\frac{\varphi }{1+\varphi /k_1}>\frac{\varphi }{1+\varphi /k_2}$ based on k₁ > k₂ > 0, we have $\frac{m_1^{\ast }}{k_1}<\frac{m_2^{\ast }}{k_2}$; that is, $\frac{k_2}{k_1}<\frac{m_2^{\ast }}{m_1^{\ast }}$. Therefore, $\frac{k_2}{k_1}<\frac{m_2^{\ast }}{m_1^{\ast }}<1$. (ii) can be obtained similarly. For (iii), it is easy to obtain $m_1^{\ast }=m_2^{\ast }$ by the fact of k₁ = k₂ > 0 and ${\mu }_1(m_i^{\ast }{f}_i(\varphi ,0))=0$.

As the consequence of Lemma 3.2, system (2.2) admits the following trivial and semi-trivial solutions: trivial solution (0,0); semi-trivial solution (ω_∗(⋅; m₁, β₁), 0) exists if and only if $m_1>m_1^{\ast }$; semi-trivial solution (0, ω_∗(⋅; m₂, β₂)) exists if and only if $m_2>m_2^{\ast }$. For convenience, we denote $\widehat{u}={\omega }_{\ast }(\cdot ;m_1,{\beta }_1)$, $\widehat{v}={\omega }_{\ast }(\cdot ;m_2,{\beta }_2)$, and next, we give an a priori estimate for the positive steady-state solution of system (2.2).

Lemma 4.1 Suppose that (u(x), v(x)) is a non-negative solution of system (2.2) with u≢0 and v≢0 on [0,1]. Then,

(i) $0<u(x)\le \widehat{u}$ and $0<v(x)\le \widehat{v}$ for x ∈ [0, 1]

(ii) 0 < u(x) + v(x) < ϕ(x) for x ∈ [0, 1]

The proof of Lemma 4.1 is exactly similar to that in [10], Lemma 4.2; hence, it is omitted here.

We next establish the linear stability of $(\widehat{u},0)$ and $(0,\widehat{v})$. First, the linearized operator of system (2.3) at $(\widehat{u},0)$ is given by

$\left\{\begin{array}{cc}d{\phi }_{xx}+m_1\left[f_1\left(\varphi -\widehat{u},\widehat{u}\right)+\widehat{u}{f}_1^u\left(\varphi -\widehat{u},\widehat{u}\right)\right]\phi +m_1\widehat{u}{f}_1^v\left(\varphi -\widehat{u},\widehat{u}\right)\psi =\mu \phi ,\hfill & x\in \left(\mathrm{0,1}\right),\hfill \\ d{\psi }_{xx}+m_2{f}_2\left(\varphi -\widehat{u},0\right)\psi =\mu \psi ,\hfill & x\in \left(\mathrm{0,1}\right),\hfill \\ {\phi }_x\left(0\right)={\phi }_x\left(1\right)+\gamma \phi \left(1\right)=0,\hfill \\ {\psi }_x\left(0\right)={\psi }_x\left(1\right)+\gamma \psi \left(1\right)=0,\hfill \end{array}\right.$(4.3)

where $f_1^u(\varphi -\widehat{u},\widehat{u})=-\frac{k_1+{\beta }_1\varphi }{{[k_1+\varphi +({\beta }_1-1)\widehat{u}]}^2}<0,f_1^v(\varphi -\widehat{u},\widehat{u})=-\frac{k_1+{\beta }_1\widehat{u}}{{[k_1+\varphi +({\beta }_1-1)\widehat{u}]}^2}<0$. By the Riesz–Schauder theory, the eigenvalues of Eq. 4.3 consist of the eigenvalues of the following two operators:

${\mathcal{B}}_1\left(m_1\right)=d\frac{{\mathrm{d}}^2}{\mathrm{d}x}+m_1\left(f_1\left(\varphi -\widehat{u},\widehat{u}\right)+\widehat{u}{f}_1^u\left(\varphi -\widehat{u},\widehat{u}\right)\right),{\mathcal{B}}_2\left(m_2\right)=d\frac{{\mathrm{d}}^2}{\mathrm{d}x}+m_2{f}_2\left(\varphi -\widehat{u},0\right),$(4.4)

subject to the corresponding boundary conditions. It follows from Lemma 3.1 (iii) that ${\mu }_1(m_1(f_1(\varphi -\widehat{u},\widehat{u})+\widehat{u}{f}_1^u(\varphi -\widehat{u},\widehat{u})))<{\mu }_1(m_1{f}_1(\varphi -\widehat{u},\widehat{u}))$. Moreover, ${\mu }_1(m_1{f}_1(\varphi -\widehat{u},\widehat{u}))=0$ with eigenfunction ${\phi }_1(m_1{f}_1(\varphi -\widehat{u},\widehat{u}))=\widehat{u}$, which implies ${\mu }_1(m_1\left(f_1(\varphi -\widehat{u},\widehat{u})+\widehat{u}{f}_1^u(\varphi -\widehat{u},\widehat{u})\right))<0$. Hence, the stability of $(\widehat{u},0)$ is determined by the sign of the principal eigenvalue ${\mu }_1(m_2{f}_2(\varphi -\widehat{u},0))$ of ${\mathcal{B}}_2(m_2)$. More precisely, $(\widehat{u},0)$ is asymptotically stable if ${\mu }_1(m_2{f}_2(\varphi -\widehat{u},0))<0$, while $(\widehat{u},0)$ is unstable if ${\mu }_1(m_2{f}_2(\varphi -\widehat{u},0))>0$.

The linearized operator of system (2.3) at $(0,\widehat{v})$ is given by

$\left\{\begin{array}{cc}d{\phi }_{xx}+m_1{f}_1\left(\varphi -\widehat{v},0\right)\phi =\mu \phi ,\hfill & x\in \left(\mathrm{0,1}\right),\hfill \\ d{\psi }_{xx}+m_2\left[f_2\left(\varphi -\widehat{v},\widehat{v}\right)+\widehat{v}{f}_2^v\left(\varphi -\widehat{v},\widehat{v}\right)\right]\psi +m_2\widehat{v}{f}_2^u\left(\varphi -\widehat{v},\widehat{v}\right)\phi =\mu \psi ,\hfill & x\in \left(\mathrm{0,1}\right),\hfill \\ {\phi }_x\left(0\right)={\phi }_x\left(1\right)+\gamma \phi \left(1\right)=0,\hfill \\ {\psi }_x\left(0\right)={\psi }_x\left(1\right)+\gamma \psi \left(1\right)=0,\hfill \end{array}\right.$(4.5)

where $f_2^u(\varphi -\widehat{v},\widehat{v})=-\frac{k_2+{\beta }_2\widehat{v}}{{[k_2+\varphi +({\beta }_2-1)\widehat{v}]}^2}<0,f_2^v(\varphi -\widehat{v},\widehat{v})=-\frac{k_2+{\beta }_2\varphi }{{[k_2+\varphi +({\beta }_2-1)\widehat{v}]}^2}<0$. We denote

${\mathcal{B}}_3\left(m_2\right)=d\frac{{\mathrm{d}}^2}{\mathrm{d}x}+m_1{f}_1\left(\varphi -\widehat{v},0\right),{\mathcal{B}}_4\left(m_2\right)=d\frac{{\mathrm{d}}^2}{\mathrm{d}x}+m_2\left[f_2\left(\varphi -\widehat{v},\widehat{v}\right)+\widehat{v}{f}_2^v\left(\varphi -\widehat{v},\widehat{v}\right)\right].$(4.6)

Similarly, $(0,\widehat{v})$ is asymptotically stable if ${\mu }_1(m_1{f}_1(\varphi -\widehat{v},0))<0$, while $(0,\widehat{v})$ is unstable if ${\mu }_1(m_1{f}_1(\varphi -\widehat{v},0))>0$.

Theorem 4.1 We consider d, k₁, k₂ > 0 fixed. Let (u(x, t), v(x, t)) be the solution of (2.2) with any non-negative non-trivial initial condition. The following statements hold:

(i) We consider β₁, β₂ > 0 fixed.

(i.1) If $m_1\le m_1^{\star }$ and $m_2\le m_2^{\star }$, then

$\underset{t\to +\infty }{\lim }u\left(x,t\right)=0,\underset{t\to +\infty }{\lim }v\left(x,t\right)=0\text{uniformly on}\left[\mathrm{0,1}\right].$(4.7)

(i.2) If $m_1\le m_1^{\star }$ and $m_2>m_2^{\star }$, then

$\underset{t\to +\infty }{\lim }u\left(x,t\right)=0,\underset{t\to +\infty }{\lim }v\left(x,t\right)=\widehat{v}\text{uniformly on}\left[\mathrm{0,1}\right].$(4.8)

(i.3) If $m_1>m_1^{\star }$ and $m_2\le m_2^{\star }$, then

$\underset{t\to +\infty }{\lim }u\left(x,t\right)=\widehat{u},\underset{t\to +\infty }{\lim }v\left(x,t\right)=0\text{uniformly on}\left[\mathrm{0,1}\right].$(4.9)

(ii) We consider $m_1>m_1^{\star }$, $m_2>\max \{m_2^{\star },\max \{\frac{k_2}{k_1},1\}{m}_1\}$, and β₁ > 0 fixed. Then, $(\widehat{u},0)$ is unstable. Moreover,

(ii.1) (4.8) holds provided

$0<{\beta }_2\le {\beta }_2^0≔\frac{\left(m_2{k}_1-m_1{k}_2\right)\gamma }{m_1{S}^0\left(1+\gamma \right)},$(4.10)

(ii.2) there exists a unique ${\beta }_2^{\star }\in ({\beta }_2^0,+\infty )$ such that $(0,\widehat{v})$ is asymptotically stable when $0<{\beta }_2<{\beta }_2^{\star }$; $(0,\widehat{v})$ is unstable when ${\beta }_2>{\beta }_2^{\star }$, and system (2.2) admits at least one stable coexistence steady state when ${\beta }_2>{\beta }_2^{\star }$.

(iii) We consider $m_2>m_2^{\star }$, $m_1>\max \{m_1^{\star },\max \{\frac{k_1}{k_2},1\}{m}_2\}$, and β₂ > 0 fixed. Then, $(0,\widehat{v})$ is unstable. Moreover,

(iii.1) Eq. 4.9 holds provided

$0<{\beta }_1\le {\beta }_1^0≔\frac{\left(m_1{k}_2-m_2{k}_1\right)\gamma }{m_2{S}^0\left(1+\gamma \right)},$(4.11)

(iii.2) there exists a unique ${\beta }_1^{\star }\in ({\beta }_1^0,+\infty )$ such that $(\widehat{u},0)$ is asymptotically stable when $0<{\beta }_1<{\beta }_1^{\star }$; $(\widehat{u},0)$ is unstable when ${\beta }_1>{\beta }_1^{\star }$, and system (2.2) admits at least one stable coexistence steady state when ${\beta }_1>{\beta }_1^{\star }$.

Proof (i) can be proved by the similar arguments as in [6], Theorems 3.5, 3.6, and we omit it here. Next, we only prove (ii), since (iii) can be proved by similar arguments.

Claim 1. For $m_1>m_1^{\ast }$, $m_2\ge \max \{\frac{k_2}{k_1},1\}{m}_1$, and β₁ > 0 fixed, $(\widehat{u},0)$ is unstable.

Note that $\widehat{u}$ satisfies

$\left\{\begin{array}{cc}d{\widehat{u}}_{xx}+m_1{f}_1\left(\varphi -\widehat{u},\widehat{u}\right)\widehat{u}=0,\hfill & x\in \left(\mathrm{0,1}\right),\hfill \\ {\widehat{u}}_x\left(0\right)={\widehat{u}}_x\left(1\right)+\gamma \widehat{u}\left(1\right)=0,\hfill \end{array}\right.$(4.12)

which implies ${\mu }_1(m_1{f}_1(\varphi -\widehat{u},\widehat{u}))=0$. We recall that $(\widehat{u},0)$ is asymptotically stable if ${\mu }_1(m_2{f}_2(\varphi -\widehat{u},0))<0$ and it is unstable if ${\mu }_1(m_2{f}_2(\varphi -\widehat{u},0))>0$. Then, we conclude from $m_2\ge \max \{\frac{k_2}{k_1},1\}{m}_1$ and $0<\widehat{u}<\varphi$ on [0,1] that

$\begin{array}{cc}\hfill m_1{f}_1\left(\varphi -\widehat{u},\widehat{u}\right)-m_2{f}_2\left(\varphi -\widehat{u},0\right)& <m_1{f}_1\left(\varphi -\widehat{u},0\right)-m_2{f}_2\left(\varphi -\widehat{u},0\right)\hfill \\ \hfill & =\frac{m_1{k}_2-m_2{k}_1+\left(m_1-m_2\right)\left(\varphi -\widehat{u}\right)}{\left(k_1+\varphi -\widehat{u}\right)\left(k_2+\varphi -\widehat{u}\right)}\left(\varphi -\widehat{u}\right)\hfill \\ \hfill & \le 0,\hfill \end{array},$

which means that ${\mu }_1(m_2{f}_2(\varphi -\widehat{u},0))>{\mu }_1(m_1{f}_1(\varphi -\widehat{u},\widehat{u}))=0$ by Lemma 3.2(iii). That is, $(\widehat{u},0)$ is unstable.

Claim 2. (1) For $m_1>m_1^{\ast }$, $m_2>\max \{m_2^{\ast },\max \{\frac{k_2}{k_1},1\}{m}_1\}$, and β₁ > 0 fixed, $(0,\widehat{v})$ is asymptotically stable when $0<{\beta }_2\le {\beta }_2^0$, where ${\beta }_2^0$ is defined by Eq. 4.10.

(2) There exists a unique ${\beta }_2^{\ast }>{\beta }_2^0$ such that $(0,\widehat{v})$ is asymptotically stable when $0<{\beta }_2<{\beta }_2^{\ast }$, and $(0,\widehat{v})$ is unstable when ${\beta }_2>{\beta }_2^{\ast }$.

For (1), we recall that $(0,\widehat{v})$ is asymptotically stable if ${\mu }_1(m_1{f}_1(\varphi -\widehat{v},0))<0$ and unstable if ${\mu }_1(m_1{f}_1(\varphi -\widehat{v},0))>0$. Similarly, we can conclude from $m_2>\max \{m_2^{\ast },\max \{\frac{k_2}{k_1},1\}{m}_1\}$ and $0<\widehat{v}<\varphi (x)<\frac{S^0(1+\gamma )}{\gamma }$ on [0,1] that

$\begin{array}{cc}\hfill & m_1{f}_1\left(\varphi -\widehat{v},0\right)-m_2{f}_2\left(\varphi -\widehat{v},\widehat{v}\right)\hfill \\ \hfill =& \frac{m_1{k}_2-m_2{k}_1+\left(m_1-m_2\right)\left(\varphi -\widehat{v}\right)+m_1{\beta }_2\widehat{v}}{\left[k_2+\varphi +\left({\beta }_2-1\right)\widehat{v}\right]\left(k_1+\varphi -\widehat{v}\right)}\left(\varphi -\widehat{v}\right)\hfill \\ \hfill <& 0\hfill \end{array},$

provided $0<{\beta }_2\le {\beta }_2^0≔\frac{(m_2{k}_1-m_1{k}_2)\gamma }{m_1{S}^0(1+\gamma )}$. It follows from Lemma 3.2(iii) that ${\mu }_1(m_1{f}_1(\varphi -\widehat{v},0))<{\mu }_1(m_2{f}_2(\varphi -\widehat{v},\widehat{v}))=0$ when $0<{\beta }_2\le {\beta }_2^0$. That is, $(0,\widehat{v})$ is asymptotically stable when $0<{\beta }_2\le {\beta }_2^0$.

For (2), since $\widehat{v}(\cdot ;m_2,{\beta }_2)$ is point-wise strictly decreasing in β₂ ∈ (0, + ∞) and $\underset{{\beta }_2\to +\infty }{\lim }\widehat{v}(\cdot ;m_2,{\beta }_2)=0$ uniformly on [0,1] (see Lemma 3.3(ii)), we can obtain from Lemma 3.2(iii) that ${\mu }_1(m_1{f}_1(\varphi -\widehat{v}(\cdot ;m_2,{\beta }_2),0))$ is strictly increasing in β₂ ∈ (0, + ∞). Furthermore,

$\underset{{\beta }_2\to +\infty }{\lim }{\mu }_1\left(m_1{f}_1\left(\varphi -\widehat{v}\left(\cdot ;m_2,{\beta }_2\right),0\right)\right)={\mu }_1\left(m_1{f}_1\left(\varphi ,0\right)\right)>0,$

based on $m_1>m_1^{\ast }$ (Eq. 4.1). Moreover, ${\mu }_1(m_1{f}_1(\varphi -\widehat{v}(\cdot ;m_2,{\beta }_2),0))<0$ when $0<{\beta }_2\le {\beta }_2^0$. Therefore, there exists a unique ${\beta }_2^{\ast }>{\beta }_2^0$ such that

${\mu }_1\left(m_1{f}_1\left(\varphi -\widehat{v}\left(\cdot ;m_2,{\beta }_2\right),0\right)\right)\left\{\begin{array}{c}<0\mathrm{i}\mathrm{f}0<{\beta }_2<{\beta }_2^{\ast },\hfill \\ =0\mathrm{i}\mathrm{f}{\beta }_2={\beta }_2^{\ast },\hfill \\ >0\mathrm{i}\mathrm{f}{\beta }_2>{\beta }_2^{\ast },\hfill \end{array}\right.$

which means that $(0,\widehat{v})$ is asymptotically stable when $0<{\beta }_2<{\beta }_2^{\ast }$, while it is unstable when ${\beta }_2>{\beta }_2^{\ast }$.

Claim 3. For $m_1>m_1^{\ast }$, $m_2>\max \{m_2^{\ast },\max \{\frac{k_2}{k_1},1\}{m}_1\}$, and β₁ > 0 fixed, system (2.2) has no positive steady states when $0<{\beta }_2\le {\beta }_2^0$.

We assume by contradiction that system (2.2) admits a positive steady state $(\bar{u},\bar{v})$, which satisfies

$\left\{\begin{array}{cc}d{\bar{u}}_{xx}+m_1{f}_1\left(\varphi -\bar{u}-\bar{v},\bar{u}\right)\bar{u}=0,\hfill & x\in \left(\mathrm{0,1}\right),\hfill \\ d{\bar{v}}_{xx}+m_2{f}_2\left(\varphi -\bar{u}-\bar{v},\bar{v}\right)\bar{v}=0,\hfill & x\in \left(\mathrm{0,1}\right),\hfill \\ {\bar{u}}_x\left(0\right)={\bar{u}}_x\left(1\right)+\gamma \bar{u}\left(1\right)=0,\hfill \\ {\bar{v}}_x\left(0\right)={\bar{v}}_x\left(1\right)+\gamma \bar{v}\left(1\right)=0.\hfill \end{array}\right.$(4.13)

Multiplying the first equation of (4.13) by $\bar{v}$ and the second equation by $\bar{u}$, integrating over (0,1), and then, subtracting the resulting equations, we have

${\int }_0^1\frac{m_1{k}_2-m_2{k}_1+\left(m_1-m_2\right)\left(\varphi -\bar{u}-\bar{v}\right)+m_1{\beta }_2\bar{v}-m_2{\beta }_1\bar{u}}{\left[k_1+\varphi +\left({\beta }_1-1\right)\bar{u}-\bar{v}\right]\left[k_2+\varphi +\left({\beta }_2-1\right)\bar{v}-\bar{u}\right]}\bar{u}\bar{v}\left(\varphi -\bar{u}-\bar{v}\right)\mathrm{d}x=0.$(4.14)

Since $\bar{u}+\bar{v}<\varphi (x)<\frac{S^0(1+\gamma )}{\gamma }$ on [0,1], we can conclude from $m_1>m_1^{\ast }$ and $m_2>\max \{m_2^{\ast },\max \{\frac{k_2}{k_1},1\}{m}_1\}$ that the left side of (4.14) is less than 0, when $0<{\beta }_2\le {\beta }_2^0≔\frac{(m_2{k}_1-m_1{k}_2)\gamma }{m_1{S}^0(1+\gamma )}$. That is a contradiction.

In conclusion, we can deduce that (ii.1) holds from Claim 1, Claim 2(1), Claim 3, and Lemma 3.4(iii). In addition, (ii.2) is the direct result of Claim 1, Claim 2(2), and Lemma 3.4(ii). The proof is completed.

Remark 4.1 Theorem 4.1(i) implies that both species with sufficiently small growth rates are washed out, while competition exclusion occurs and the species with a sufficiently faster growth rate will finally win the competition. In particular, when both species admit sufficient fast growth rates, Theorem 4.1(ii.1) suggests that the species v with stronger growth ability ($\frac{m_2}{k_2}$ is large) and weaker intraspecific competition (β₂ is small) will finally win the competition. This is consistent with the biological intuition that the species with stronger growth ability and weaker intraspecific competition has more competitive advantages. Theorem 4.1(iii.1) illustrates the similar biological phenomenon.

We next investigate the local dynamics of system (2.2). Note that the stability of $(\widehat{u}(m_1,{\beta }_1),0)$ is determined by the sign of ${\mu }_1(m_2{f}_2(\varphi -\widehat{u}(m_1,{\beta }_1),0))$, and the stability of $(0,\widehat{v}(m_2,{\beta }_2))$ is determined by the sign of ${\mu }_1(m_1{f}_1(\varphi -\widehat{v}(m_1,{\beta }_1),0))$. Clearly, ${\mu }_1(m_2{f}_2(\varphi -\widehat{u}(m_1,{\beta }_1),0))$ depends on m₁, m₂, and β₁, and ${\mu }_1(m_1{f}_1(\varphi -\widehat{v}(m_2,{\beta }_2),0))$ depends on m₁, m₂, and β₂. To this end, we define

${\sigma }_1\left(m_1,m_2,{\beta }_1\right)≔{\mu }_1\left(m_2{f}_2\left(\varphi -\widehat{u}\left(m_1,{\beta }_1\right),0\right)\right)\mathrm{f}\mathrm{o}\mathrm{r}{m}_1>m_1^{\ast },m_2>0,{\beta }_1>0,$

${\tau }_1\left(m_1,m_2,{\beta }_2\right)≔{\mu }_1\left(m_1{f}_1\left(\varphi -\widehat{v}\left(m_2,{\beta }_2\right),0\right)\right)\mathrm{f}\mathrm{o}\mathrm{r}{m}_1>0,m_2>m_2^{\ast },{\beta }_2>0.$

Lemma 4.2 The principal eigenvalues σ₁(m₁, m₂, β₁) and τ₁(m₁, m₂, β₂) have the following properties:

(i) For fixed d, k₁, k₂ > 0 and $m_1>m_1^{\ast }$,

(i.1) σ₁(m₁, m₂, β₁) is strictly decreasing with respect to m₁ in $(m_1^{\ast },+\infty )$,

(i.2) σ₁(m₁, m₂, β₁) is strictly increasing with respect to β₁ in (0, + ∞),

(i.3) σ₁(m₁, m₂, β₁) is strictly increasing with respect to m₂ in $(m_2^{\ast },+\infty )$; moreover,

$\underset{m_2\to {\left(m_2^{\ast }\right)}^{+}}{\lim }{\sigma }_1\left(m_1,m_2,{\beta }_1\right)={\mu }_1\left(m_2^{\ast }{f}_2\left(\varphi -\widehat{u},0\right)\right)<0,\underset{m_2\to +\infty }{\lim }{\sigma }_1\left(m_1,m_2,{\beta }_1\right)=+\infty .$

(ii) For fixed d, k₁, k₂ > 0 and $m_2>m_2^{\ast }$,

(ii.1) τ₁(m₁, m₂, β₂) is strictly increasing with respect to m₁ in $(m_1^{\ast },+\infty )$,

(ii.2) τ₁(m₁, m₂, β₂) is strictly increasing with respect to β₂ in (0, + ∞),

(ii.3) τ₁(m₁, m₂, β₂) is strictly decreasing with respect to m₂ in $(m_2^{\ast },+\infty )$; moreover,

$\underset{m_2\to {\left(m_2^{\ast }\right)}^{+}}{\lim }{\tau }_1\left(m_1,m_2,{\beta }_2\right)={\mu }_1\left(m_1{f}_1\left(\varphi ,0\right)\right)>0,\underset{m_2\to +\infty }{\lim }{\tau }_1\left(m_1,m_2,{\beta }_2\right)={\mu }_1\left(0\right)<0.$

Proof For (i), (i.1) can be obtained by Lemma 3.1(iii) and Lemma 3.3(i). Similarly, (i.2) is followed by Lemma 3.1(iii) and Lemma 3.3(ii). To prove (i.3), it is obvious that σ₁(m₁, m₂, β₁) is strictly increasing with respect to m₂ in $(m_2^{\ast },+\infty )$ by Lemma 3.1(iii) and $\underset{m_2\to +\infty }{\lim }{\sigma }_1(m_1,m_2,{\beta }_1)=+\infty$ by (3.2). We recall ${\mu }_1(m_2^{\ast }{f}_2(\varphi ,0))=0$. Then, we can conclude from Lemma 3.1(ii) (iii) that

$\underset{m_2\to {\left(m_2^{\ast }\right)}^{+}}{\lim }{\sigma }_1\left(m_1,m_2,{\beta }_1\right)={\mu }_1\left(m_2^{\ast }{f}_2\left(\varphi -\widehat{u},0\right)\right)<{\mu }_1\left(m_2^{\ast }{f}_2\left(\varphi ,0\right)\right)=0.$

For (ii), (ii.1) can be obtained by Lemma 3.1(iii), and (ii.2) can be proved by Lemma 3.1(iii) and Lemma 3.3(ii). We then prove (ii.3). Since $\widehat{v}(m_2,{\beta }_2)$ is point-wise strictly increasing in $m_2\in (m_2^{\ast },+\infty )$ by Lemma 3.3(i), it follows from Lemma 3.1(iii) that τ₁(m₁, m₂, β₂) is strictly decreasing with respect to m₂ in $(m_2^{\ast },+\infty )$. Noting that $\underset{m_2\to {(m_2^{\ast })}^{+}}{\lim }\widehat{v}(m_2,{\beta }_2)=0$ by (3.5), we can conclude from Lemma 3.1(ii) that

$\underset{m_2\to {\left(m_2^{\ast }\right)}^{+}}{\lim }{\tau }_1\left(m_1,m_2,{\beta }_2\right)={\mu }_1\left(m_1{f}_1\left(\varphi ,0\right)\right)>0,$

based on $m_1>m_1^{\ast }$ (see (4.1)). Moreover, since $\underset{m_2\to +\infty }{\lim }\widehat{v}(m_2,{\beta }_2)=\varphi (x)$ on [0,1] (see (3.5)), we can obtain from Lemma 3.1(ii) (iii) that

$\underset{m_2\to +\infty }{\lim }{\tau }_1\left(m_1,m_2,{\beta }_2\right)={\mu }_1\left(m_1{f}_1\left(\varphi -\widehat{v}\left(m_2,{\beta }_2\right),0\right)\right)={\mu }_1\left(0\right)<0.$

The proof is completed.

Clearly, both σ₁(m₁, m₂, β₁) and τ₁(m₁, m₂, β₂) depend on m₁ and m₂. To investigate the local dynamics of system (2.2) in the m₁ − m₂ plane, we fix β₁, β₂ > 0 and denote them by σ₁(m₁, m₂) and τ₁(m₁, m₂).

Lemma 4.3 Suppose $m_i>m_i^{\ast }(i=\mathrm{1,2})$. For fixed d, β₁, β₂, k₁, k₂ > 0, there exist two continuous critical curves

${\mathrm{\Gamma }}_1:m_2={\bar{m}}_2\left(m_1\right)\mathit{\text{for}}{m}_1\in \left(m_1^{\ast },+\infty \right),$

${\mathrm{\Gamma }}_2:m_2={\tilde{m}}_2\left(m_1\right)\mathit{\text{for}}{m}_1\in \left(m_1^{\ast },+\infty \right),$

where ${\bar{m}}_2(m_1)$ and ${\tilde{m}}_2(m_1)$ are differentially dependent on m₁ and uniquely determined by

${\mu }_1\left({\bar{m}}_2\left(m_1\right)f_2\left(\varphi -\widehat{u},0\right)\right)=0\mathit{\text{and}}{\mu }_1\left(m_1{f}_1\left(\varphi -\widehat{v}\left(\cdot ;{\tilde{m}}_2\left(m_1\right)\right),0\right)\right)=0,$(4.15)

respectively. Then,

(i) the semi-trivial steady state $(\widehat{u},0)$ is locally asymptotically stable if $(m_1,m_2)\in (m_1^{\ast },+\infty )\times (m_2^{\ast },{\bar{m}}_2(m_1))$, neutrally stable if $(m_1,m_2)\in (m_1^{\ast },+\infty )\times \{{\bar{m}}_2(m_1)\}$, and unstable if.

(ii) the semi-trivial steady state $(0,\widehat{v})$ is locally asymptotically stable if $(m_1,m_2)\in (m_1^{\ast },+\infty )\times \{{\tilde{m}}_2(m_1),+\infty \}$, neutrally stable if $(m_1,m_2)\in (m_1^{\ast },+\infty )\times \{{\tilde{m}}_2(m_1)\}$, and unstable if $(m_1,m_2)\in (m_1^{\ast },+\infty )\times (m_2^{\ast },{\tilde{m}}_2(m_1))$.

Proof (i) By Lemma 4.2(i), we conclude that for any $m_1\in (m_1^{\ast },+\infty )$ given, there exists a unique ${\bar{m}}_2(m_1)>m_2^{\ast }$ such that

${\sigma }_1\left(m_1,m_2\right)\left\{\begin{array}{c}<0\mathrm{i}\mathrm{f}{m}_2^{\ast }<m_2<{\bar{m}}_2\left(m_1\right),\hfill \\ =0\mathrm{i}\mathrm{f}{m}_2={\bar{m}}_2\left(m_1\right),\hfill \\ >0\mathrm{i}\mathrm{f}{m}_2>{\bar{m}}_2\left(m_1\right).\hfill \end{array}\right.$

Therefore, $(\widehat{u},0)$ is locally asymptotically stable for $m_2^{\ast }<m_2<{\bar{m}}_2(m_1)$, neutrally stable for $m_2={\bar{m}}_2(m_1)$, and unstable for $m_2>{\bar{m}}_2(m_1)$.

(ii) Similarly, we can conclude from Lemma 4.2(ii) that for any $m_1\in (m_1^{\ast },+\infty )$ given, there exists a unique ${\tilde{m}}_2(m_1)>m_2^{\ast }$ such that

${\tau }_1\left(m_1,m_2\right)\left\{\begin{array}{c}>0\mathrm{i}\mathrm{f}{m}_2^{\ast }<m_2<{\tilde{m}}_2\left(m_1\right),\hfill \\ =0\mathrm{i}\mathrm{f}{m}_2={\tilde{m}}_2\left(m_1\right),\hfill \\ <0\mathrm{i}\mathrm{f}{m}_2>{\tilde{m}}_2\left(m_1\right).\hfill \end{array}\right.$

Therefore, $(0,\widehat{v})$ is unstable for $m_2^{\ast }<m_2<{\tilde{m}}_2(m_1)$, neutrally stable for $m_2={\tilde{m}}_2(m_1)$, and locally asymptotically stable for $m_2>{\tilde{m}}_2(m_1)$. The proof is completed.

Combining with Lemma 3.4 and Lemma 4.3, we obtain the following results.

Theorem 4.2 Suppose $m_i>m_i^{\ast }(i=\mathrm{1,2})$. For fixed d, k₁, k₂, β₁, β₂ > 0, ${\bar{m}}_2$ and ${\tilde{m}}_2$ are defined by Eq. 4.15, respectively, and the following statements hold.

(i) Suppose ${\bar{m}}_2<{\tilde{m}}_2$. If $m_2<{\bar{m}}_2$, then $(\widehat{u},0)$ is locally asymptotically stable, and $(0,\widehat{v})$ is unstable if it exists. If $m_2>{\tilde{m}}_2$, then $(\widehat{u},0)$ is unstable, and $(0,\widehat{v})$ is locally asymptotically stable. If $m_2\in ({\bar{m}}_2,{\tilde{m}}_2)$, then $(\widehat{u},0)$ and $(0,\widehat{v})$ are both unstable, and system (2.2) admits at least one stable coexistence steady state.

(ii) Suppose ${\bar{m}}_2>{\tilde{m}}_2$. If $m_2<{\tilde{m}}_2$, then $(\widehat{u},0)$ is locally asymptotically stable, and $(0,\widehat{v})$ is unstable. If $m_2>{\bar{m}}_2$, then $(\widehat{u},0)$ is unstable, and $(0,\widehat{v})$ is locally asymptotically stable. If $m_2\in ({\tilde{m}}_2,{\bar{m}}_2)$, then $(\widehat{u},0)$ and $(0,\widehat{v})$ are both stable, and system (2.2) admits at least one unstable coexistence steady state.

(iii) Suppose ${\bar{m}}_2={\tilde{m}}_2$. If $m_2<{\tilde{m}}_2$, then $(\widehat{u},0)$ is locally asymptotically stable, and $(0,\widehat{v})$ is unstable. If $m_2>{\tilde{m}}_2$, then $(\widehat{u},0)$ is unstable, and $(0,\widehat{v})$ is locally asymptotically stable.

For fixed d, k₁, k₂, β₁, β₂ > 0, Lemma 4.3 and Theorem 4.2 imply that there exist two critical curves Γ₁ and Γ₂ in the m₁ − m₂ plane, which divide the local dynamics of Eq. 2.2 into competitive exclusion, bi-stability, and coexistence. To further characterize classification on the dynamics of system (2.2) in the m₁ − m₂ plane, we next give some properties of critical curves Γ₁ and Γ₂. We recall

${\mu }_1\left({\bar{m}}_2{f}_2\left(\varphi -\widehat{u}\left(\cdot ;m_1,{\beta }_1\right),0\right)\right)=0\mathrm{a}\mathrm{n}\mathrm{d}{\mu }_1\left(m_1{f}_1\left(\varphi -\widehat{v}\left(\cdot ;{\tilde{m}}_2,{\beta }_2\right),0\right)\right)=0.$

Clearly, ${\bar{m}}_2$ depends on m₁ and β₁ and ${\tilde{m}}_2$ depends on m₁ and β₂. To emphasize these parameters, we denote ${\bar{m}}_2$ and ${\tilde{m}}_2$ by ${\bar{m}}_2(m_1,{\beta }_1)$ and ${\tilde{m}}_2(m_1,{\beta }_2)$, respectively.

Proposition 4.2 We consider d, k₁, k₂ > 0 fixed. The critical curve ${\bar{m}}_2(m_1,{\beta }_1)$ has the following properties.

(i) For any β₁ > 0 given, ${\bar{m}}_2(m_1,{\beta }_1)$ is strictly increasing with respect to m₁ in $(m_1^{\ast },+\infty )$. Moreover,

$\underset{m_1\to {\left(m_1^{\ast }\right)}^{+}}{\lim }{\bar{m}}_2\left(m_1,{\beta }_1\right)=m_2^{\ast }.$

(ii) For any $m_1>m_1^{\ast }$ given, ${\bar{m}}_2(m_1,{\beta }_1)$ is strictly decreasing with respect to β₁ in (0, + ∞) and

$m_2^{\ast }<{\bar{m}}_2\left(m_1,{\beta }_1\right)<\max \left\{\frac{k_2}{k_1},1\right\}{m}_1\text{and}\underset{{\beta }_1\to +\infty }{\lim }{\bar{m}}_2\left(m_1,{\beta }_1\right)=m_2^{\ast }.$(4.16)

Particularly,

$\max \left\{m_2^{\ast },\min \left\{\frac{k_2}{k_1},1\right\}{m}_1\right\}\le {\bar{m}}_2\left(m_1,{\beta }_1\right)<\max \left\{\frac{k_2}{k_1},1\right\}{m}_1\mathrm{\text{provided}}0<{\beta }_1\le {\beta }_1^0.$(4.17)

(iii) For any $m_1>m_1^{\ast }$ given,

$\underset{{\beta }_1\to +\infty }{\lim }{\dot{\stackrel{̄}{m}}}_2\left(m_1,{\beta }_1\right)=0,$(4.18)

where ${\dot{\stackrel{̄}{m}}}_2(m_1,{\beta }_1)$ is the derivative of ${\bar{m}}_2(m_1,{\beta }_1)$ with respect to m₁ in $(m_1^{\ast },+\infty )$.

Proof (i) For any β₁ > 0 given, we can conclude from Lemma 4.2(i.1) (i.3) that ${\mu }_1(m_2{f}_2(\varphi -\widehat{u}),0)$ is strictly decreasing with respect to m₁ in $(m_1^{\ast },+\infty )$ and strictly increasing with respect to m₂ in $(m_2^{\ast },+\infty )$. Then, it follows from ${\mu }_1({\bar{m}}_2(m_1,{\beta }_1)f_2(\varphi -\widehat{u},0))=0$ and the implicit function theorem that ${\bar{m}}_2(m_1,{\beta }_1)$ is strictly increasing with respect to m₁ in $(m_1^{\ast },+\infty )$. Since ${\mu }_1({\bar{m}}_2(m_1,{\beta }_1)f_2(\varphi -\widehat{u}),0)=0$, we can obtain $\underset{m_1\to {(m_1^{\ast })}^{+}}{\lim }{\bar{m}}_2(m_1,{\beta }_1)=m_2^{\ast }$ based on ${\mu }_1(m_2^{\ast }{f}_2(\varphi ,0))=0$ and $\underset{m_1\to {(m_1^{\ast })}^{+}}{\lim }\widehat{u}=0$ uniformly on [0,1] (see Eq. 3.5).

(ii) For any $m_1>m_1^{\ast }$ given, it follows from Lemma 4.2(i.2) (i.3) that ${\mu }_1(m_2{f}_2(\varphi -\widehat{u}),0)$ is strictly increasing with respect to β₁ in (0, + ∞) and strictly increasing with respect to m₂ in $(m_2^{\ast },+\infty )$. Then, we can obtain from ${\mu }_1({\bar{m}}_2(m_1,{\beta }_1)f_2(\varphi -\widehat{u},0))=0$ and the implicit function theorem that ${\bar{m}}_2(m_1,{\beta }_1)$ is strictly decreasing with respect to β₁ in (0, + ∞). Then,

$\underset{{\beta }_1\to +\infty }{\lim }{\bar{m}}_2\left(m_1,{\beta }_1\right)<{\bar{m}}_2\left(m_1,{\beta }_1\right)<\underset{{\beta }_1\to 0^{+}}{\lim }{\bar{m}}_2\left(m_1,{\beta }_1\right)\mathrm{f}\mathrm{o}\mathrm{r}{m}_1>m_1^{\ast }.$

Substituting β₁ → 0⁺ in ${\mu }_1({\bar{m}}_2(m_1,{\beta }_1)f_2(\varphi -\widehat{u},0))=0$, we have ${\mu }_1({\bar{m}}_2(m_1,0)f_2(\varphi -u_0,0))=0$ by Lemma 3.1(ii). Here, u₀ satisfies

$\left\{\begin{array}{cc}d{\left(u_0\right)}_{xx}+\frac{m_1\left(\varphi -u_0\right)u_0}{k_1+\varphi -u_0}=0,\hfill & x\in \left(\mathrm{0,1}\right),\hfill \\ {\left(u_0\right)}_x\left(0\right)={\left(u_0\right)}_x\left(1\right)+\gamma u_0\left(1\right)=0.\hfill \end{array}\right.$

Then, we can deduce from [26], Theorem 2.1 that $\underset{{\beta }_1\to 0^{+}}{\lim }{\bar{m}}_2(m_1,{\beta }_1)\le \max \{\frac{k_2}{k_1},1\}{m}_1$. Substituting β₁ → +∞ in ${\mu }_1({\bar{m}}_2(m_1,{\beta }_1)f_2(\varphi -\widehat{u},0))=0$, we have $\underset{{\beta }_1\to +\infty }{\lim }{\bar{m}}_2(m_1,{\beta }_1)=m_2^{\ast }$ based on ${\mu }_1(m_2^{\ast }{f}_2(\varphi ,0))=0$ and $\underset{{\beta }_1\to +\infty }{\lim }\widehat{u}=0$ uniformly on [0,1] (see (3.6)). The estimate of ${\bar{m}}_2(m_1,{\beta }_1)$ in Eq. 4.16 is obtained.

Finally, when $m_1>m_1^{\ast },m_2^{\ast }<m_2<\min \{\frac{k_2}{k_1},1\}{m}_1$, Theorem 4.1(iii.1) shows that $(\widehat{u},0)$ is globally asymptotically stable provided $0<{\beta }_1\le {\beta }_1^0$. Therefore, ${\bar{m}}_2(m_1,{\beta }_1)\ge \max \{m_2^{\ast },\min \{\frac{k_2}{k_1},1\}{m}_1\}$ provided $0<{\beta }_1\le {\beta }_1^0$. Combining with Eq. 4.16, 4.17 is obtained.

(iii) Let ψ₁ > 0 with ‖ψ₁‖_∞ = 1 be the corresponding principal eigenfunction of ${\mu }_1({\bar{m}}_2{f}_2(\varphi -\widehat{u},0))=0$. Then, ψ₁ satisfies

$\left\{\begin{array}{cc}d{\left({\psi }_1\right)}_{xx}+{\bar{m}}_2{f}_2\left(\varphi -\widehat{u},0\right){\psi }_1=0,\hfill & x\in \left(\mathrm{0,1}\right),\hfill \\ {\left({\psi }_1\right)}_x\left(0\right)={\left({\psi }_1\right)}_x\left(1\right)+\gamma {\psi }_1\left(1\right)=0.\hfill \end{array}\right.$(4.19)

By differentiating Eq. 4.19 with respect to m₁, denoting $\frac{\partial }{\partial m_1}=\stackrel{̇}{}$, we have

$\left\{\begin{array}{c}d{\left(\dot{{\psi }_1}\right)}_{xx}+{\dot{\stackrel{̄}{m}}}_2{f}_2\left(\varphi -\widehat{u},0\right){\psi }_1+{\bar{m}}_2\hfill \\ \left[f_2\left(\varphi -\widehat{u},0\right)\dot{{\psi }_1}-\frac{k_2}{{\left(k_2+\varphi -\widehat{u}\right)}^2}\cdot \frac{\partial \widehat{u}\left(m_1,{\beta }_1\right)}{\partial m_1}{\psi }_1\right]=0,& x\in \left(\mathrm{0,1}\right),\hfill \\ {\left(\dot{{\psi }_1}\right)}_x\left(0\right)={\left(\dot{{\psi }_1}\right)}_x\left(1\right)+\gamma \dot{{\psi }_1}\left(1\right)=0,\hfill \end{array}\right.$(4.20)

where $\frac{\partial \widehat{u}(m_1,{\beta }_1)}{\partial m_1}$ is the derivative of $\widehat{u}(m_1,{\beta }_1)$ with respect to m₁ in $(m_1^{\ast },+\infty )$. Multiplying Eq. 4.19 by $\dot{{\psi }_1}$ and Eq. 4.20 by ψ₁, integrating over (0,1), and then, subtracting the two resulting equations,

${\dot{\stackrel{̄}{m}}}_2\left(m_1,{\beta }_1\right)=\frac{{\bar{m}}_2{\int }_0^1\frac{k_2}{{\left(k_2+\varphi -\widehat{u}\right)}^2}\cdot \frac{\partial \widehat{u}\left(m_1,{\beta }_1\right)}{\partial m_1}{\psi }_1^2\mathrm{d}x}{{\int }_0^1{f}_2\left(\varphi -\widehat{u},0\right){\psi }_1^2\mathrm{d}x}.$(4.21)

We next show that $\underset{{\beta }_1\to +\infty }{\lim }\frac{\partial \widehat{u}(m_1,{\beta }_1)}{\partial m_1}=\frac{\partial }{\partial m_1}\underset{{\beta }_1\to +\infty }{\lim }\widehat{u}(m_1,{\beta }_1)=0$. We choose an increasing sequence {β_1,n} with $\underset{n\to +\infty }{\lim }{\beta }_{1,n}=+\infty$. Then, ${\widehat{u}}_n≔\widehat{u}(m_1,{\beta }_{1,n})$ is the unique positive solution of

$\left\{\begin{array}{cc}d{\left({\widehat{u}}_n\right)}_{xx}+m_1{f}_1^{\left(n\right)}\left(\varphi -{\widehat{u}}_n,{\widehat{u}}_n\right){\widehat{u}}_n=0,\hfill & x\in \left(\mathrm{0,1}\right),\hfill \\ {\left({\widehat{u}}_n\right)}_x\left(0\right)={\left({\widehat{u}}_n\right)}_x\left(1\right)+\gamma \left({\widehat{u}}_n\right)\left(1\right)=0,\hfill \end{array}\right.$(4.22)

where $f_1^{(n)}(\varphi -{\widehat{u}}_n,{\widehat{u}}_n)=\frac{\varphi -{\widehat{u}}_n}{k_1+\varphi +({\beta }_{1,n}-1){\widehat{u}}_n}$. Note that $0<{\widehat{u}}_n<\varphi$ and $f_1^{(n)}(\varphi -{\widehat{u}}_n,{\widehat{u}}_n)$ is uniformly bounded on [0,1]. By L^p estimates for p ∈ (1, + ∞) and the Sobolev embedding theorem, we may assume that $\underset{n\to +\infty }{\lim }{\widehat{u}}_n=U_1(x)$ in C¹[0, 1] by passing to a subsequence if necessary. Since ${\widehat{u}}_n$ is continuously differentiable with respect to m₁ in $(m_1^{\ast },+\infty )$ (see Lemma 3.3(i)), we differentiate (4.22) with respect to m₁, denote $\frac{\partial }{\partial m_1}=\stackrel{̇}{}$, and obtain

$\left\{\begin{array}{c}d{\left({\dot{\widehat{u}}}_n\right)}_{xx}+m_1\hfill \\ \left[f_1^{\left(n\right)}\left(\varphi -{\widehat{u}}_n,{\widehat{u}}_n\right)-\frac{\left(k_1+{\beta }_{1,n}\varphi \right){\widehat{u}}_n}{{\left[k_1+\varphi +\left({\beta }_{1,n}-1\right){\widehat{u}}_n\right]}^2}\right]{\dot{\widehat{u}}}_n=-f_1^{\left(n\right)}\left(\varphi -{\widehat{u}}_n,{\widehat{u}}_n\right){\widehat{u}}_n,& x\in \left(\mathrm{0,1}\right),\hfill \\ {\left({\dot{\widehat{u}}}_n\right)}_x\left(0\right)={\left({\dot{\widehat{u}}}_n\right)}_x\left(1\right)+\gamma {\dot{\widehat{u}}}_n\left(1\right)=0.\hfill \end{array}\right.$(4.23)

Since $0<{\widehat{u}}_n<\varphi$, $f_1^{(n)}(\varphi -{\widehat{u}}_n,{\widehat{u}}_n){\widehat{u}}_n$, and $f_1^{(n)}(\varphi -{\widehat{u}}_n,{\widehat{u}}_n)-\frac{(k_1+{\beta }_{1,n}\varphi ){\widehat{u}}_n}{{[k_1+\varphi +({\beta }_{1,n}-1){\widehat{u}}_n]}^2}$ are uniformly bounded on [0,1], we can use L^p estimates for p ∈ (1, + ∞) and the Sobolev embedding theorem again to assume that $\underset{n\to +\infty }{\lim }{\dot{\widehat{u}}}_n=U_2(x)$ in C¹[0, 1] by passing to a subsequence if necessary. Therefore, $\underset{{\beta }_1\to +\infty }{\lim }\frac{\partial \widehat{u}(m_1,{\beta }_1)}{\partial m_1}=\frac{\partial }{\partial m_1}\underset{{\beta }_1\to +\infty }{\lim }\widehat{u}(m_1,{\beta }_1)=0$ based on $\underset{{\beta }_1\to +\infty }{\lim }\widehat{u}(m_1,{\beta }_1)=0$ uniformly on [0,1] (see Lemma 3.3(ii)). Finally, taking β₁ → +∞ in (4.21) and noting that $\underset{{\beta }_1\to +\infty }{\lim }{\bar{m}}_2(m_1,{\beta }_1)=m_2^{\ast }$ (see (4.16)), we have $\underset{{\beta }_1\to +\infty }{\lim }{\dot{\stackrel{̄}{m}}}_2(m_1,{\beta }_1)=0$. The proof is completed.

Proposition 4.3 We consider d, k₁, k₂ > 0 fixed. The critical curve ${\tilde{m}}_2(m_1,{\beta }_2)$ has the following properties:

(i) For any β₂ > 0 given, ${\tilde{m}}_2(m_1,{\beta }_2)$ is strictly increasing with respect to m₁ in $(m_1^{\ast },+\infty )$. Moreover,

$\underset{m_1\to {\left(m_1^{\ast }\right)}^{+}}{\lim }{\tilde{m}}_2\left(m_1,{\beta }_2\right)=m_2^{\ast }.$

(ii) For any $m_1>m_1^{\ast }$ given, ${\tilde{m}}_2(m_1,{\beta }_2)$ is strictly increasing with respect to β₂ in (0, + ∞) and

$\max \left\{m_2^{\ast },\min \left\{\frac{k_2}{k_1},1\right\}{m}_1\right\}<{\tilde{m}}_2\left(m_1,{\beta }_2\right)<+\infty .$(4.24)

Particularly,

$\max \left\{m_2^{\ast },\min \left\{\frac{k_2}{k_1},1\right\}{m}_1\right\}<{\tilde{m}}_2\left(m_1,{\beta }_2\right)\le \max \left\{\frac{k_2}{k_1},1\right\}{m}_1\text{provided}0<{\beta }_2\le {\beta }_2^0.$(4.25)

(iii) For any $m_1>m_1^{\ast }$ given,

$\underset{{\beta }_2\to +\infty }{\lim }{\dot{\stackrel{̃}{m}}}_2\left(m_1,{\beta }_2\right)=+\infty ,$(4.26)

where ${\dot{\stackrel{̃}{m}}}_2(m_1,{\beta }_2)$ is the derivative of ${\tilde{m}}_2(m_1,{\beta }_2)$ with respect to m₁ in $(m_1^{\ast },+\infty )$.

Proof (i) For any β₂ > 0 given, we can conclude from Lemma 4.2(ii.1) (ii.3) that ${\mu }_1(m_1{f}_1(\varphi -\widehat{v}(\cdot ;m_2,{\beta }_2),0))$ is strictly increasing with respect to m₁ in $(m_1^{\ast },+\infty )$ and strictly decreasing with respect to m₂ in $(m_2^{\ast },+\infty )$. Then, ${\tilde{m}}_2(m_1,{\beta }_2)$ is strictly increasing with respect to m₁ in $(m_1^{\ast },+\infty )$ by ${\mu }_1(m_1{f}_1(\varphi -\widehat{v}(\cdot ;{\tilde{m}}_2,{\beta }_2),0))=0$ and the implicit function theorem. Substituting $m_1\to {(m_1^{\ast })}^{+}$ in ${\mu }_1(m_1{f}_1(\varphi -\widehat{v}(\cdot ;{\tilde{m}}_2,{\beta }_2),0))=0$, we have $\underset{m_1\to {(m_1^{\ast })}^{+}}{\lim }\widehat{v}(\cdot ;{\tilde{m}}_2,{\beta }_2)=0$ based on ${\mu }_1(m_1^{\ast }{f}_1(\varphi ,0))=0$. Therefore, $\underset{m_1\to {(m_1^{\ast })}^{+}}{\lim }{\tilde{m}}_2(m_1,{\beta }_2)=m_2^{\ast }$ by $\underset{{\tilde{m}}_2\to {(m_2^{\ast })}^{+}}{\lim }\widehat{v}(\cdot ;{\tilde{m}}_2,{\beta }_2)=0$ uniformly on [0,1] (Eq. 3.5).

(ii) For any $m_1>m_1^{\ast }$ given, it follows from Lemma 4.2(ii.2) (ii.3) that ${\mu }_1(m_1{f}_1(\varphi -\widehat{v}(\cdot ;m_2,{\beta }_2),0))$ is strictly increasing with respect to β₂ in (0, + ∞) and strictly decreasing with respect to m₂ in $(m_2^{\ast },+\infty )$. Then, we can conclude from ${\mu }_1(m_1{f}_1(\varphi -\widehat{v}(\cdot ;{\tilde{m}}_2,{\beta }_2),0))=0$ and the implicit function theorem that ${\tilde{m}}_2(m_1,{\beta }_2)$ is strictly increasing with respect to β₂ in (0, + ∞). Therefore,

$\underset{{\beta }_2\to 0^{+}}{\lim }{\tilde{m}}_2\left(m_1,{\beta }_2\right)<{\tilde{m}}_2\left(m_1,{\beta }_2\right)<\underset{{\beta }_2\to +\infty }{\lim }{\tilde{m}}_2\left(m_1,{\beta }_2\right)\mathrm{f}\mathrm{o}\mathrm{r}{m}_1>m_1^{\ast }.$

Similarly, substituting β₂ → 0⁺ in ${\mu }_1(m_1{f}_1(\varphi -\widehat{v}(\cdot ;{\tilde{m}}_2,{\beta }_2),0))=0$, we can deduce from [8], Theorem 2.1 again that $\max \{m_2^{\ast },\min \{\frac{k_2}{k_1},1\}{m}_1\}\le \underset{{\beta }_2\to 0^{+}}{\lim }{\tilde{m}}_2(m_1,{\beta }_2)$. Moreover, $\underset{{\beta }_2\to +\infty }{\lim }{\tilde{m}}_2(m_1,{\beta }_2)=+\infty$. Hence, the inequality about ${\tilde{m}}_2(m_1,{\beta }_2)$ in Eq. 4.24 is obtained.

When $m_1>m_1^{\ast }$, $m_2>\max \{m_2^{\ast },\max \{\frac{k_2}{k_1},1\}{m}_1\}$, Theorem 4.1(ii.1) implies that $(0,\widehat{v})$ is globally asymptotically stable provided $0<{\beta }_2\le {\beta }_2^0$. Therefore, ${\tilde{m}}_2(m_1,{\beta }_2)\le \max \{\frac{k_2}{k_1},1\}{m}_1$ provided $0<{\beta }_2\le {\beta }_2^0$. Combining with 4.24, 4.25 holds.

(iii) Let φ₁ with ‖φ₁‖_∞ = 1 be the corresponding principal eigenfunction of ${\mu }_1(m_1{f}_1(\varphi -\widehat{v}({\tilde{m}}_2,{\beta }_2),0))=0$. Then, φ₁ satisfies

$\left\{\begin{array}{cc}d{\left({\phi }_1\right)}_{xx}+m_1{f}_1\left(\varphi -\widehat{v}\left({\tilde{m}}_2,{\beta }_2\right),0\right){\phi }_1=0,\hfill & x\in \left(\mathrm{0,1}\right),\hfill \\ {\left({\phi }_1\right)}_x\left(0\right)={\left({\phi }_1\right)}_x\left(1\right)+\gamma {\phi }_1\left(1\right)=0.\hfill \end{array}\right.$(4.27)

By differentiating Eq. 4.27 with respect to m₁ and denoting $\frac{\partial }{\partial m_1}=\stackrel{̇}{}$, we have

$\left\{\begin{array}{c}d{\left(\dot{{\phi }_1}\right)}_{xx}+f_1\left(\varphi -\widehat{v},0\right){\phi }_1+m_1\hfill \\ \left[f_1\left(\varphi -\widehat{v},0\right)\dot{{\phi }_1}-\frac{k_1}{{\left(k_1+\varphi -\widehat{v}\right)}^2}\cdot \frac{\partial \widehat{v}\left({\tilde{m}}_2,{\beta }_2\right)}{\partial m_2}\times {\dot{\stackrel{̃}{m}}}_2\left(m_1\right){\phi }_1\right]=0,& x\in \left(\mathrm{0,1}\right),\hfill \\ {\left(\dot{{\phi }_1}\right)}_x\left(0\right)={\left(\dot{{\phi }_1}\right)}_x\left(1\right)+\gamma \dot{{\phi }_1}\left(1\right)=0,\hfill \end{array}\right.$(4.28)

where $\frac{\partial \widehat{v}({\tilde{m}}_2,{\beta }_2)}{\partial m_2}$ is the derivative of $\widehat{v}$ with respect to m₂ at $m_2={\tilde{m}}_2$. Multiplying (4.27) by $\dot{{\phi }_1}$ and (4.28) by φ₁, integrating over (0,1), and then, subtracting the two resulting equations,

${\dot{\stackrel{̃}{m}}}_2\left(m_1,{\beta }_2\right)=\frac{{\int }_0^1{f}_1\left(\varphi -\widehat{v},0\right){\phi }_1^2\mathrm{d}x}{m_1{\int }_0^1\frac{k_1}{{\left(k_1+\varphi -\widehat{v}\right)}^2}\cdot \frac{\partial \widehat{v}\left({\tilde{m}}_2,{\beta }_2\right)}{\partial m_2}{\phi }_1^2\mathrm{d}x}.$(4.29)

Similar to Proposition 4.2(iii), we can show $\underset{{\beta }_2\to +\infty }{\lim }\frac{\partial \widehat{v}({\tilde{m}}_2,{\beta }_2)}{\partial m_2}=\frac{\partial }{\partial m_2}{\lim }_{{\beta }_2\to +\infty }\widehat{v}({\tilde{m}}_2,{\beta }_2)=0$ based on $\underset{{\beta }_2\to +\infty }{\lim }\widehat{v}({\tilde{m}}_2,{\beta }_2)=0$ uniformly on [0,1] (see Lemma 3.3(ii)). Substituting β₂ → +∞ in (4.29), it is easy to obtain $\underset{{\beta }_2\to +\infty }{\lim }{\dot{\stackrel{̃}{m}}}_2(m_1,{\beta }_2)=+\infty$.

We assume k₁ > k₂ > 0 without loss of generality. Then, there exist six critical curves: $m_1=m_1^{\ast },m_2=m_2^{\ast }$,

$L_1:m_2=\frac{k_2}{k_1}{m}_1,L_2:m_1=m_2,{\mathrm{\Gamma }}_1:m_2={\bar{m}}_2\left(m_1\right),{\mathrm{\Gamma }}_2:m_2={\tilde{m}}_2\left(m_1\right),$

in the m₁ − m₂ plane (Figure 2), which classify the dynamics of system (2.2) into extinction of both species, competitive exclusion and coexistence. Clearly, it follows from Proposition 4.1(i) that line $m_2=\frac{m_2^{\ast }}{m_1^{\ast }}{m}_1$ is located above line L₁ and below line L₂ under the assumption k₁ > k₂ > 0. Propositions 4.2(i) and 4.3(i) suggest that Γ₁ and Γ₂ are increasing with respect to $m_1\in (m_1^{\ast },+\infty )$, respectively. Moreover, $\underset{m_1\to {(m_1^{\ast })}^{+}}{\lim }{\bar{m}}_2(m_1)=\underset{m_1\to {(m_1^{\ast })}^{+}}{\lim }{\tilde{m}}_2(m_1)=m_2^{\ast }$, which implies that Γ₁ and Γ₂ intersect at point $(m_1^{\ast },m_2^{\ast })$. In addition, due to the effect of β₁ and β₂, Propositions 4.2(ii) (iii) and 4.3(ii) (iii) indicate that the locations of Γ₁ and Γ₂ in the m₁ − m₂ plane have the following four occasions shown in Figures 2A–D. Here, we assume that Γ₂ is always located above Γ₁ in this region (Figure 2). We set

${\mathrm{\Pi }}_0≔\left(0,m_1^{\ast }\right]\times \left(0,m_2^{\ast }\right];{\mathrm{\Pi }}_1≔\left(m_1^{\ast },+\infty \right)\times \left(0,m_2^{\ast }\right];{\mathrm{\Pi }}_2≔\left(0,m_1^{\ast }\right]\times \left(m_2^{\ast },+\infty \right).$

Now, we are ready to illustrate the dynamical classification of system (2.2) in the m₁ − m₂ plane under the assumption k₁ > k₂ > 0, by dividing the following four cases.

Case I: $0<{\beta }_1\le {\beta }_1^0$, $0<{\beta }_2\le {\beta }_2^0$ (Figure 2A). Then, (4.17) holds provided $0<{\beta }_1\le {\beta }_1^0$; that is, $\frac{k_2}{k_1}{m}_1<{\bar{m}}_2(m_1,{\beta }_1)<m_1$ under the assumption k₁ > k₂ > 0. Similarly, (4.25) holds provided $0<{\beta }_2\le {\beta }_2^0$, which means $\frac{k_2}{k_1}{m}_1<{\tilde{m}}_2(m_1,{\beta }_1)<m_1$. These suggest that both Γ₁ and Γ₂ lie between lines L₁ and L₂ (Figure 2A). It follows from Theorem 4.1(i.1) that (0,0) is g.a.s in region (m₁, m₂) ∈ Π₀. In particular, the phase portrait graph of system (2.2) with (m₁, m₂) = (0.2, 0.1) ∈ Π₀ is illustrated in Figure 1A, which shows that (0,0) is g.a.s in this case. Then, $(\widehat{u},0)$ is g.a.s when (m₁, m₂) ∈ Π₁ ∪ Π₃ by Theorem 4.1(i.3), (iii.1). Particularly, the phase portrait graph of system (2.2) with (m₁, m₂) = (1, 0.1) ∈ Π₁ ∪ Π₃ is shown in Figure 1B. $(0,\widehat{v})$ is g.a.s in region (m₁, m₂) ∈ Π₂ ∪ Π₄ by Theorem 4.1(i.2) and (ii.1). Moreover, the specific phase portrait graph with (m₁, m₂) = (0.2, 1) ∈ Π₂ ∪ Π₄ is displayed in Figure 1C. Furthermore, by Lemma 4.3, $(\widehat{u},0)$ is locally asymptotically stable in region Π₅ and unstable in region Π₆ ∪ Π₇; $(0,\widehat{v})$ is locally asymptotically stable in region Π₆ and unstable in region Π₅ ∪ Π₇. Then, we can conclude from Theorem 4.2(i) that there exist stable coexistence steady states in Π₇, and the specific phase portrait graph with (m₁, m₂) = (1, 0.545) ∈ Π₇ is presented in Figure 1D.

Case II: $0<{\beta }_1\le {\beta }_1^0$, ${\beta }_2>{\beta }_2^{\ast }$ (Figure 2B). Then, (4.17) holds and Γ₁ still lies between lines L₁ and L₂ when $0<{\beta }_1\le {\beta }_1^0$. Moreover, (4.24) holds when ${\beta }_2>{\beta }_2^{\ast }$; that is, ${\tilde{m}}_2(m_1,{\beta }_2)>\max \{m_2^{\ast },\frac{k_2}{k_1}{m}_1\}$ under the assumption k₁ > k₂ > 0. This implies that Γ₂ lies above lines L₁ and $m_2=m_2^{\ast }$. Note that $\underset{m_1\to {(m_1^{\ast })}^{+}}{\lim }{\tilde{m}}_2(m_1,{\beta }_2)=m_2^{\ast }$ by Proposition 4.3(i) and $m_1^{\ast }>m_2^{\ast }$ if k₁ > k₂ > 0 (see Proposition 4.1). Then, ${\tilde{m}}_2(m_1,{\beta }_2)<m_1$ when m₁ is near $m_1^{\ast }$, which means that there exists sufficiently small ϵ > 0 such that Γ₂ lies below line L₂ for $m_1\in (m_1^{\ast },m_1^{\ast }+ϵ)$.

Next, we illustrate two-fold that Γ₂ will first intersect L₂ and then lie above L₂ as m₁ increases (Figure 2B). On one hand, Theorem 4.1(ii) indicates that when $m_1>m_1^{\ast }$ and $m_2>\max \{m_2^{\ast },m_1\}$, there exists a unique ${\beta }_2^{\ast }$ such that $(0,\widehat{v})$ is unstable when ${\beta }_2>{\beta }_2^{\ast }$. This means that there exist regions above L₂ such that $(0,\widehat{v})$ is unstable when ${\beta }_2>{\beta }_2^{\ast }$. Note that Γ₂ is strictly increasing with respect to β₂ in (0, + ∞) by Proposition 4.3(ii). Hence, when ${\beta }_2>{\beta }_2^{\ast }$, we can conclude from Lemma 4.3(ii) that the critical curve Γ₂ for the stability change of $(0,\widehat{v})$ must intersect L₂ and then lie above L₂ as m₁ is suitably large. This implies that there exist regions above line L₂ and below Γ₂ such that $(0,\widehat{v})$ is unstable when ${\beta }_2>{\beta }_2^{\ast }$. On the other hand, Proposition 4.3(iii) implies $\underset{{\beta }_2\to +\infty }{\lim }{\dot{\stackrel{̃}{m}}}_2(m_1,{\beta }_2)=+\infty$, which means that the slope of Γ₂ will be bigger than 1 for suitably large β₂. This, in turn, suggests that for large β₂, there exists ${\bar{m}}_1^{\ast }>m_1^{\ast }$ such that Γ₂ must intersect L₂ and then locates above L₂ for all $m_1>{\bar{m}}_1^{\ast }$. Therefore, the region ${\mathrm{\Pi }}_7^1$ occurs for large β₂. This region is denoted as ${\mathrm{\Pi }}_7^1$ in Figure 2B. Then, we deduce from Theorem 4.2(i) that there exist stable coexistence steady states in ${\mathrm{\Pi }}_7={\mathrm{\Pi }}_7^0\cup {\mathrm{\Pi }}_7^1$. The other regions Π₀ − Π₆ can be similarly defined as in Case I.

Case III: ${\beta }_1>{\beta }_1^{\ast }$, $0<{\beta }_2\le {\beta }_2^0$ (see Figure 2C). Similar to Case II, (4.25) holds provided $0<{\beta }_2\le {\beta }_2^0$; that is, line Γ₂ lies between lines L₁ and L₂ when $0<{\beta }_2\le {\beta }_2^0$. Moreover, (4.16) holds when ${\beta }_1>{\beta }_1^{\ast }$; that is, $m_2^{\ast }<{\bar{m}}_2(m_1,{\beta }_1)<m_1$ under the assumption k₁ > k₂ > 0. This implies that Γ₁ lies below line L₂ and above line $m_2=m_2^{\ast }$. Furthermore, note that $\underset{m_1\to {(m_1^{\ast })}^{+}}{\lim }{\bar{m}}_2(m_1,{\beta }_1)=m_2^{\ast }$ by Proposition 4.2(i) and $m_2^{\ast }>\frac{k_2}{k_1}{m}_1^{\ast }$ if k₁ > k₂ > 0 by Proposition 4.1. Then, ${\bar{m}}_2(m_1,{\beta }_1)>\frac{k_2}{k_1}{m}_1$, when m₁ is near $m_1^{\ast }$, which means that there exists sufficiently small ɛ > 0 such that Γ₁ lies above L₁ for $m_1\in (m_1^{\ast },m_1^{\ast }+\epsilon )$.

Similarly, we next illustrate that Γ₁ will first intersect L₁ and then lie below L₁ as m₁ increases (see Figure 2C). On one hand, Theorem 4.1(iii) suggests that when $m_1>m_1^{\ast }$ and $m_2^{\ast }<m_2<\frac{k_2}{k_1}{m}_1$, there exists a unique ${\beta }_1^{\ast }$ such that $(\widehat{u},0)$ is unstable when ${\beta }_1>{\beta }_1^{\ast }$. This means that there exist regions below L₁ such that $(\widehat{u},0)$ is unstable when ${\beta }_1>{\beta }_1^{\ast }$. Note that Γ₁ is strictly decreasing with respect to β₁ in (0, + ∞) by Proposition 4.2(ii). Hence, when ${\beta }_1>{\beta }_1^{\ast }$, we can deduce from Lemma 4.3(i) that the critical curve Γ₁ for the stability change of $(\widehat{u},0)$ must intersect L₁ and then lie below L₁ as m₁ is suitably large. This implies that there exist regions below L₁ and above Γ₁ such that $(\widehat{u},0)$ is unstable when ${\beta }_1>{\beta }_1^{\ast }$.

On the other hand, Proposition 4.2(iii) gives $\underset{{\beta }_1\to +\infty }{\lim }{\dot{\stackrel{̄}{m}}}_2(m_1,{\beta }_1)=0$, which implies that the slope of Γ₁ will be smaller than $\frac{k_2}{k_1}$ for suitably large β₁. This, in turn, suggests that, for large β₁, there exists ${\bar{m}}_1^{\ast \ast }>m_1^{\ast }$ such that Γ₁ must intersect L₁ and then lie below L₁ for all $m_1>{\bar{m}}_1^{\ast \ast }$. Therefore, the region ${\mathrm{\Pi }}_7^2$ occurs for large β₁. This region is denoted as ${\mathrm{\Pi }}_7^2$ in Figure 2C. Then, we obtain from Theorem 4.2(i) that there exist stable coexistence steady states in ${\mathrm{\Pi }}_7={\mathrm{\Pi }}_7^0\cup {\mathrm{\Pi }}_7^2$. The other regions Π₀ − Π₆ can be similarly illustrated as in Case I.

Case IV: ${\beta }_1>{\beta }_1^{\ast }$, ${\beta }_2>{\beta }_2^{\ast }$ (Figure 2D). Combining the analysis of ${\beta }_1>{\beta }_1^{\ast }$ in Case II and ${\beta }_2>{\beta }_2^{\ast }$ in Case III, we know that both the regions ${\mathrm{\Pi }}_7^1$ and ${\mathrm{\Pi }}_7^2$ exist. It follows from Theorem 4.2(i) that there exist stable coexistence steady states in ${\mathrm{\Pi }}_7={\mathrm{\Pi }}_7^0\cup {\mathrm{\Pi }}_7^1\cup {\mathrm{\Pi }}_7^2$. The other regions Π₀ − Π₆ can be similarly illustrated as in Case I.

We next make a comparison with the results in [8]. When the intraspecific competition is relatively weak (i.e., for the case of $0<{\beta }_1\le {\beta }_1^0$ and $0<{\beta }_2\le {\beta }_2^0$), the competitive dynamics of system (2.2) (Figure 2A) are similar to the unstirred chemostat model with Holling type II functional response (see [8], Theorems 2.1, 2.4] and Figure 1 in [8]), which suggests that the weak intraspecific competition has little effect on the competition outcomes of species.

However, when the intraspecific competition becomes strong, some new phenomena may occur. For instance, for the standard unstirred chemostat models, competition exclusion always happens for the weak–strong competition of two species (see [26], Theorem 2.1), while coexistence may occur in the unstirred chemostat model with B.–D. functional response with the increase of intraspecific competition, under the weak–strong competition cases (see Theorem 4.1 and Figures 2B–D). More precisely, under the weak–strong competition case $m_2>\max \{\frac{k_2}{k_1},1\}{m}_1$ (i.e., species v has a stronger growth ability compared to species u), the competitive ability of species v becomes weak with the increase of β₂, which may result in the coexistence of the two species (see Theorem 4.1(ii.2) and Figure 2B). Similarly, under the weak–strong competition case $m_2<\min \{\frac{k_2}{k_1},1\}{m}_1$ (i.e., species u has a stronger growth ability compared to species v), the competitive ability of species u becomes weak with the increase of β₁, which may cause the coexistence of the two species (see Theorem 4.1(iii.2) and Figure 2C). In particular, Theorem 4.1(ii.2), (iii.2) and Figure 2D suggest that coexistence is more likely to happen, when the intraspecific competition of two species is strong.

In summary, these theoretical results indicate that for the weak–strong competition cases, if the intraspecific competition parameter of the species with stronger growth ability is suitably large, we can observe different results from [8] that coexistence may occur. This new phenomenon suggests that the intraspecific competition parameters β₁ and β₂ have a great influence on the competitive outcomes of two species.

FIGURE 1

FIGURE 1. Phase portrait graphs of system (2.2) for different growth rates m₁ and m₂. Here, we take β₁ = 0.01, β₂ = 0.01, L = 1, S⁰ = 1, d = 0.5, γ = 0.5, k₁ = 1, and k₂ = 0.4. As shown, (0,0) is globally asymptotically stable (simplified as g.a.s) in (A) with m₁ = 0.2, m₂ = 0.1; $(\widehat{u},0)$ is g.a.s. in (B) with m₁ = 1, m₂ = 0.1; $(0,\widehat{v})$ is g.a.s. in (C) with m₁ = 0.2, m₂ = 1; and there exist stable coexistence steady states in (D) with m₁ = 1, m₂ = 0.545.

FIGURE 2

FIGURE 2. Illustration of the dynamics of system (2.2) in the m₁ − m₂ plane for the case of k₁ > k₂ > 0. More precisely, (A) $0<{\beta }_1\le {\beta }_1^0$ and $0<{\beta }_2\le {\beta }_2^0$; (B) $0<{\beta }_1\le {\beta }_1^0$ and ${\beta }_2>{\beta }_2^{\star }$; (C) ${\beta }_1>{\beta }_1^{\star }$ and $0<{\beta }_2\le {\beta }_2^0$; and (D) ${\beta }_1>{\beta }_1^{\star }$ and ${\beta }_2>{\beta }_2^{\star }$. Then, (0,0) is g.a.s in region Π₀; $(\widehat{u},0)$ is g.a.s in region Π₁ ∪ Π₃; $(0,\widehat{v})$ is g.a.s in region Π₂ ∪ Π₄; $(\widehat{u},0)$ is locally asymptotically stable in region Π₅ and unstable in region Π₆ ∪ Π₇; $(0,\widehat{v})$ is locally asymptotically stable in region Π₆ and unstable in region Π₅ ∪ Π₇; and there exist stable coexistence steady states in Π₇. Here, we note that (B) ${\mathrm{\Pi }}_7={\mathrm{\Pi }}_7^0\cup {\mathrm{\Pi }}_7^1$, (C) ${\mathrm{\Pi }}_7={\mathrm{\Pi }}_7^0\cup {\mathrm{\Pi }}_7^2$, and (D) ${\mathrm{\Pi }}_7={\mathrm{\Pi }}_7^0\cup {\mathrm{\Pi }}_7^1\cup {\mathrm{\Pi }}_7^2$.

5 Positive solution branches of system (2.2)

We define X = W^2,p(0, 1) × W^2,p(0, 1) and Y = L^p(0, 1) × L^p(0, 1), where p > 1. For fixed d, k₁, k₂, β₁, β₂ > 0 and $m_1>m_1^{\ast }$, system (2.2) admits three branches of trivial or semi-trivial solutions in the space ${\mathbb{R}}_{+}\times X$: ${\mathrm{\Gamma }}_0=\{(m_2,\mathrm{0,0}):m_2>0\},{\mathrm{\Gamma }}_u=\{(m_2,\widehat{u},0):m_2>0\}\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{\Gamma }}_v=\{(m_2,0,\widehat{v}):m_2>m_2^{\ast }\}$, where $\widehat{u}={\omega }^{\ast }(\cdot ;m_1,{\beta }_1)$ and $\widehat{v}={\omega }^{\ast }(\cdot ;m_2,{\beta }_2)$. In this section, we will regard m₂ as the bifurcation parameter and study separately positive solutions bifurcating from the semi-trivial branches Γ_u and Γ_v by the Crandall–Rabinowitz bifurcation theorem in [31].

We first show that there exists a positive solution branch that bifurcates from the semi-trivial solution $(\widehat{u},0)$. Moreover, the bifurcation of positive solutions from $(\widehat{u},0)$ can only occur at $m_2={\bar{m}}_2$ by Lemma 4.3.

Theorem 5.1 For fixed d, k₁, k₂, β₁, β₂ > 0 and $m_1>m_1^{\ast }$, there is a smooth non-constant solutions curve Γ₁ = {(m₂(s), u(s), v(s)): s ∈ (−ϵ, ϵ)} such that (m₂(s), u(s), v(s)) is a positive solution of system (2.2) for s ∈ (0, ϵ) and satisfies $m_2(0)={\bar{m}}_2$, $u(s)=\widehat{u}+s{\phi }_0+o(s)$, and v(s) = sψ₀ + o(s). Here, ψ₀ > 0 is the principal eigenfunction corresponding to the eigenvalue ${\mu }_1({\bar{m}}_2{f}_2(\varphi -\widehat{u},0))=0$, which satisfies

$\left\{\begin{array}{cc}{\mathcal{B}}_2\left({\bar{m}}_2\right){\psi }_0=0,\hfill & x\in \left(\mathrm{0,1}\right),\hfill \\ {\psi }_0^{\prime }\left(0\right)={\psi }_0^{\prime }\left(1\right)+\gamma {\psi }_0\left(1\right)=0,\hfill \end{array}\right.$(5.1)

and φ₀ < 0 satisfies

$\left\{\begin{array}{cc}{\mathcal{B}}_1\left(m_1\right){\phi }_0=-m_1\widehat{u}{f}_1^v\left(\varphi -\widehat{u},\widehat{u}\right){\psi }_0,\hfill & x\in \left(\mathrm{0,1}\right),\hfill \\ {\phi }_0^{\prime }\left(0\right)={\phi }_0^{\prime }\left(1\right)+\gamma {\phi }_0\left(1\right)=0.\hfill \end{array}\right.$(5.2)

Moreover,

$m_2^{\prime }\left(0\right)=-\frac{{\int }_0^1{\bar{m}}_2\left[f_2^u\left(\varphi -\widehat{u},0\right){\phi }_0+f_2^v\left(\varphi -\widehat{u},0\right){\psi }_0\right]{\psi }_0^2\text{d}x}{{\int }_0^1{f}_2\left(\varphi -\widehat{u},0\right){\psi }_0^2\text{d}x},$(5.3)

where $f_2^u(\varphi -\widehat{u},0)=-\frac{k_2+{\beta }_2\widehat{u}}{{(k_2+\varphi -\widehat{u})}^2}$ and $f_2^v(\varphi -\widehat{u},0)=-\frac{k_2+{\beta }_2\varphi }{{(k_2+\varphi -\widehat{v})}^2}$.

Theorem 5.1. can be proved by the similar arguments as in [35], Theorem 6.2. For completeness, we defer the proof of Theorem 5.1 to the Supplementary Appendix.

Theorem 5.2 For fixed d, k₁, k₂, β₁, β₂ > 0 and $m_2>m_2^{\ast }$, there is a smooth non-constant solutions curve Γ₂ = {(m₂(s), u(s), v(s)): s ∈ (−ϵ, ϵ)} such that (m₂(s), u(s), v(s)) is a positive solution of (2.2) for s ∈ (0, ϵ) and satisfies $m_2(0)={\tilde{m}}_2$, $u(s)=s{\tilde{\phi }}_0+o(s)$, and $v(s)=\widehat{v}+s{\tilde{\psi }}_0+o(s)(\widehat{v}=\widehat{v}(\cdot ;{\tilde{m}}_2,{\beta }_2))$. Here, ${\tilde{\phi }}_0>0$ is the principal eigenfunction corresponding to ${\mu }_1(m_1{f}_1(\varphi -\widehat{v}(\cdot ;{\tilde{m}}_2,{\beta }_2),0))=0$, which satisfies

$\left\{\begin{array}{cc}{\mathcal{B}}_3\left({\tilde{m}}_2\right){\tilde{\phi }}_0=0,\hfill & x\in \left(\mathrm{0,1}\right),\hfill \\ {\tilde{\phi }}_0^{\prime }\left(0\right)={\tilde{\phi }}_0^{\prime }\left(1\right)+\gamma {\tilde{\phi }}_0\left(1\right)=0,\hfill \end{array}\right.$(5.4)

and ${\tilde{\psi }}_0<0$ satisfies

$\left\{\begin{array}{cc}{\mathcal{B}}_4\left({\tilde{m}}_2\right){\tilde{\psi }}_0=-{\tilde{m}}_2\widehat{v}{f}_2^u\left(\varphi -\widehat{v},\widehat{v}\right){\tilde{\phi }}_0,\hfill & x\in \left(\mathrm{0,1}\right),\hfill \\ {\tilde{\psi }}^{\prime }\left(0\right)={\tilde{\psi }}^{\prime }\left(1\right)+\gamma \tilde{\psi }\left(1\right)=0.\hfill \end{array}\right.$(5.5)

Moreover,

$m_2^{\prime }\left(0\right)=-\frac{{\int }_0^1\left[f_1^u\left(\varphi -\widehat{v},0\right){\tilde{\phi }}_0+f_1^v\left(\varphi -\widehat{v},0\right){\tilde{\psi }}_0\right]{\tilde{\phi }}_0^2\mathit{\text{d}}x}{{\int }_0^1{\widehat{v}}^{\prime }{f}_1^v\left(\varphi -\widehat{v},0\right){\tilde{\psi }}_0^2\mathit{\text{d}}x},$(5.6)

where $f_1^u(\varphi -\widehat{v},0)=-\frac{k_1+{\beta }_1\varphi }{{(k_1+\varphi -\widehat{v})}^2}$, $f_1^v(\varphi -\widehat{v},0)=-\frac{k_1+{\beta }_1\widehat{v}}{{(k_1+\varphi -\widehat{v})}^2}$, and ${\widehat{v}}^{\prime }$ is the derivative of $\widehat{v}$ with respect to m₂ at $m_2={\tilde{m}}_2$.

The proof of Theorem 5.2 is similar to the arguments in Theorem 5.1, and we omit it here.

We define $\mathrm{\Omega }=\{(m_2,u,v)\in {\mathbb{R}}_{+}\times X:u>0,v>0,(m_2,u,v)\mathrm{s}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{s}\mathrm{fi}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{y}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{m}(2.3)\}$. The following results show that the two bifurcation continua Γ₁ and Γ₂ in Theorems 5.1 and 5.2 are connected.

Theorem 5.3 We consider d, k₁, k₂, β₁, β₂ > 0 and $m_1>m_1^{\ast }$ fixed. Then, there exists a connected component Γ of Ω, which bifurcates from the semi-trivial solution branch Γ_u at $({\bar{m}}_2,\widehat{u},0)$ and meets the other semi-trivial solution branch Γ_v at $({\tilde{m}}_2,0,\widehat{v})$. In particular, system (2.3) admits a positive solution (u, v) if m₂ lies between ${\bar{m}}_2$ and ${\tilde{m}}_2$.

The proof is motivated by the methods in [13], Theorem 6.4 (see also [8], Theorem 4.10]). For readability, the proof is given in Supplementary Appendix.

6 Numerical descriptions

In this section, we will study the effect of diffusion rates d on the dynamics of system (2.2). It follows from Eqs 4.1, 4.15 that the threshold values $m_1^{\ast }$, $m_2^{\ast }$, ${\bar{m}}_2$, and ${\tilde{m}}_2$ depend on the diffusion rates d. Since these threshold values are dependent on d and they are non-monotone in d, we cannot theoretically establish the threshold dynamics of system (2.2) in terms of d. In order to explore the effect of d on the dynamics of system (2.2), in this subsection, we resort to numerical approaches. We fix L = 1, S⁰ = 1, γ = 0.5, k₁ = 1, k₂ = 0.4, and m₁ = 1 as mentioned before. Then, by presenting the bifurcation diagrams of positive equilibrium solution of system (2.2) with the bifurcation parameter d increasing, the results are divided into the following three cases.

Case I: $m_2>\max \{\frac{k_2}{k_1},1\}{m}_1$. We take m₂ = 2 such that $\frac{m_2}{k_2}=5>\frac{m_1}{k_1}=1$, which means that species v has a stronger growth ability compared to species u. We can call this case the weak–strong competition [13]. To identify the effect of intraspecific competition, we fix β₁ = 0.01 and let β₂ change from β₂ = 0.01 to β₂ = 1.

First, if β₁ = 0.01, β₂ = 0.01, there is no coexistence and the competitive exclusion principle holds (species v with a stronger growth ability will win the competition) when d is sufficiently small (Figure 3A). As d increases, both species go extinct, which is consistent with our biological intuition that the sufficiently large diffusion rates will put species at a disadvantage. These numerical observations in Figure 3A coincide with [26], Theorem 5.4.

FIGURE 3

FIGURE 3. Bifurcation diagrams of the positive steady-state solutions to system (2.2) at t = 2000 with the bifurcation parameter d ranging from 0.01 to 10. Here, we take L = 1, S⁰ = 1, γ = 0.5, k₁ = 1, k₂ = 0.4, m₁ = 1, and m₂ = 2 and (A) β₁ = 0.01, β₂ = 0.01 and (B) β₁ = 0.01, β₂ = 1.

Second, if β₁ = 0.01, β₂ = 1. Clearly, under this weak–strong competition case, though species v has stronger growth ability compared to species u, the increase of β₂ makes the competitive ability of v weaker. This is consistent with our biological intuition that the stronger intraspecific competition will put species at a disadvantage. Moreover, coexistence may occur when d is sufficiently small (Figure 3B), which is different from [8], Theorem 5.4. As d increases, species v wins the competition. As d further increases, the sufficiently large diffusion rates drive both species to extinction.

Case II: $m_1>\max \{\frac{k_1}{k_2},1\}{m}_2$. We take m₂ = 0.2 such that $m_1=1>\max \{\frac{k_1}{k_2},1\}{m}_2=0.5$, which implies that species u has a stronger growth ability compared to species v. Similarly, we call this case the weak–strong competition case [13]. We fix β₂ = 0.01 and let β₁ change from β₁ = 0.01 to β₁ = 1. Similar to Case I, when β₁ = 0.01, β₂ = 0.01, results similar to those in Figure 3A are shown in Figure 4A. When β₁ = 1, β₂ = 0.01, the increase of β₁ makes the competitive ability of u weaker. Then, we can observe that coexistence may occur when d is sufficiently small (see Figure 4B), which is quite different from [26], Theorem 5.4.

FIGURE 4

FIGURE 4. Bifurcation diagrams of the positive steady-state solutions to system (2.2) at t = 2000 with the bifurcation parameter d ranging from 0.01 to 10. Here, we take L = 1, S⁰ = 1, γ = 0.5, k₁ = 1, k₂ = 0.4, m₁ = 1, and m₂ = 0.2 and (A) β₁ = 0.01, β₂ = 0.01 and (B) β₁ = 1, β₂ = 0.01.

Case III: $\frac{k_2}{k_1}{m}_1<m_2<m_1$. We take m₂ = 0.6 such that $\frac{k_2}{k_1}{m}_1=0.4<m_2=0.6<m_1=1$. We call this case the evenly matched competition [13]. Moreover, the fact of $\frac{m_1}{k_1}<\frac{m_2}{k_2}$ suggests that though the growth ability of the two species is evenly matched, the competitive ability of species v is still slightly better than that of species u. Then, for this evenly matched competition case, we fix β₁ = 0.01 and let β₂ change from β₂ = 0.01 to β₂ = 1.

For β₁ = 0.01, β₂ = 0.01, as shown in Figure 5A, the diffusion rates have a significant effect on the dynamics of system (2.2). More precisely, the dynamics of system (2.2) shift between four scenarios with the bifurcation parameter d increasing; that is, 1) competitive exclusion occurs and species v wins the competition, when d is sufficiently small; 2) coexistence occurs as d increases; 3) competitive exclusion occurs again and species u wins the competition, as d further increases; and 4) both species are washed out as d continues to increase. These suggest that system (2.2) may show a trade-off among extinction, exclusion, and coexistence as d increases. Particularly, coexistence occurs at the intermediate diffusion rates, which is in line with the theoretical results in [8].

FIGURE 5

FIGURE 5. Bifurcation diagrams of the positive steady state solutions to system (2.2) at t = 2000 with the bifurcation parameter d ranging from 0.01 to 5. Here we take L = 1, S⁰ = 1, γ = 0.5, k₁ = 1, k₂ = 0.4, m₁ = 1, and m₂ = 0.6 and (A) β₁ = 0.01, β₂ = 0.01 and (B) β₁ = 0.01, β₂ = 1.

For β₁ = 0.01, β₂ = 1, as stated before, the increase of β₂ will make the competitive ability of v weaker. Then, we can observe from Figure 5B that both species can coexist when d is sufficiently small. As d increases, competitive exclusion happens and species u wins the competition. As d further increases, the large diffusion rates drive two species to extinction.

In shorts, for different competition Cases (I)–(III), we investigate the effect of diffusion on the dynamics of system (2.2) by taking different intraspecific competition parameters β₁, β₂. As shown in Figures 3–5, the impacts of diffusion and intraspecific competition on the competitive outcomes of species are complex, which further suggests that diffusion and intraspecific competition play a key role in determining the dynamics of system (2.2).

7 Discussion

In this paper, we investigate an unstirred chemostat model with the Beddington–DeAngelis functional response (see system (2.2)). The analytical and numerical results show that the intraspecific competition and diffusion have an important biological effect on the dynamics of system (2.2).

Theoretically, we first adopt a basic strategy regarding the growth rates as variable/bifurcation parameters to study the effect of growth rates on system (2.2). The results show that there exist six critical curves

$m_1=m_1^{\ast },m_2=m_2^{\ast },L_1:m_2=\frac{k_2}{k_1}{m}_1,L_2:m_1=m_2,{\mathrm{\Gamma }}_1:m_2={\bar{m}}_2\left(m_1\right),{\mathrm{\Gamma }}_2:m_2={\tilde{m}}_2\left(m_1\right)$

in the m₁ − m₂ plane, which may classify the dynamics of system (2.2) into extinction of both species, competitive exclusion and coexistence (see Theorems 4.1, 4.2). To further understand the effect of β_i, (i = 1, 2) on the dynamics of (2.2), we explore the properties of critical curves Γ₁ and Γ₂ (see Propositions 4.2 and 4.3) and get a relatively clear dynamics classification of system (2.2) in the m₁ − m₂ plane (Figure 2).

Numerically, since diffusion plays a key role in determining the competition outcomes of two species, we study the effect of diffusion on the dynamics of system (2.2). More precisely, for two weak–strong competition cases, due to the effect of intraspecific competition parameters β₁ and β₂, the coexistence may occur at sufficiently small diffusion rates (Figures 3B, 4B), while for the evenly matched competition case, the dynamics of system (2.2) shift between different scenarios (competitive exclusion, coexistence, and extinction) when β₁ and β₂ are small and the coexistence only occurs at the intermediate diffusion rates (Figure 5A). When β₂ is larger than β₁, we observe from Figure 5B that coexistence may occur at sufficiently small diffusion rates.

In conclusion, in this paper, the dynamics classification of system (2.2) in the m₁ − m₂ plane is established by the linear eigenvalue theory and the monotone dynamical system theory (Figure 2). Due to the effect of intraspecific competition parameters β₁ and β₂, the dynamics of system (2.2) are more complex than that of the unstirred chemostat model with Holling type II functional response (see Figure 1 of [8]). Numerically, we study the effect of diffusion on system (2.2) and obtain rich numerical results (Figures 3–5). These numerical observations reveal that, under different competition cases, the effects of diffusion and intraspecific competition on the dynamics of system (2.2) are complex. This, in turn, suggests that the B.–D. functional response is more biologically realistic and superior to the well-known Holling type II functional response in modeling the resource uptake of species.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

Author contributions

The theory part and simulation were obtained by the first author WZ. The second and third authors HN and ZW guided the work. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the National Natural Science Foundation of China (12071270 and 12171296) and the Natural Science Basic Research Program of Shaanxi (No. 2023-JC-JQ-03).

Acknowledgments

The authors are very grateful to the referees and the handling co editor-in-chief for their kind and valuable suggestions leading to a substantial improvement of the manuscript.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2023.1205571/full#supplementary-material

References

1. Bayen T, Cazenave-Lacroutz H, Coville J. Stability of the chemostat system including a linear coupling between species. Discrete Contin Dyn Syst Ser B (2023) 28(2023):2104–29. doi:10.3934/dcdsb.2022160

CrossRef Full Text | Google Scholar

2. Hsu SB, Jin Y. The dynamics of a two host-two virus system in a chemostat environment. Discrete Contin Dyn Syst Ser B (2021) 26:415–41. doi:10.3934/dcdsb.2020298

CrossRef Full Text | Google Scholar

3. Nie H, Hsu SB, Grover JP. Algal competition in a water column with excessive dioxide in the atmosphere. J Math Biol (2016) 72:1845–92. doi:10.1007/s00285-015-0926-8

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Smith HL, Waltman P. The theory of the chemostat. Cambridge: Cambridge University, Press (1995).

Google Scholar

5. Dellal M, Bar B. Global analysis of a model of competition in the chemostat with internal inhibitor. Discrete Contin Dyn Syst Ser B (2021) 26:1129–48. doi:10.3934/dcdsb.2020156

CrossRef Full Text | Google Scholar

6. Hsu SB, Waltman P. On a system of reaction-diffusion equations arising from competition in an unstirred chemostat. SIAM J Appl Math (1993) 53:1026–44. doi:10.1137/0153051

CrossRef Full Text | Google Scholar

7. Nie H, Shi Y, Wu JH. The effect of diffusion on the dynamics of a predator-prey chemostat model. SIAM J Appl Math (2022) 3:821–48. doi:10.1137/21m1432090

CrossRef Full Text | Google Scholar

8. Shi JP, Wu YX, Zou XF. Coexistence of competing species for intermediate dispersal rates in a reaction-diffusion chemostat model. J Dynam Differential Equations (2020) 32:1085–112. doi:10.1007/s10884-019-09763-0

CrossRef Full Text | Google Scholar

9. Smith HL. Monotone dynamical systems: An introduction to the theory of competitive and cooperative systems. In: Mathematical surveys and monographs. Providence, RI: American Mathematical Society (1995).

Google Scholar

10. Wu JH. Global bifurcation of coexistence state for the competition model in the chemostat. Nonlinear Anal (2000) 39:817–35. doi:10.1016/s0362-546x(98)00250-8

CrossRef Full Text | Google Scholar

11. Holling CS. The functional response of predator to prey density and its role in mimicry and population regulation. Mem Entonmol Soc Can (1965) 45:1–60.

Google Scholar

12. Zhang W, Nie H, Wu JH. A reaction-diffusion-advection chemostat model in a flowing habitat: Mathematical analysis and numerical simulations. Internat J Bifur Chaos Appl Sci Engrg (2023) 33:1245–74. doi:10.1142/s0218127423500736

CrossRef Full Text | Google Scholar

13. Zhang W, Nie H, Wu JH. Dynamics of a reaction-diffusion-advection model with two species competing in a flow reactor. Discrete Contin Dyn Syst Ser B (2023) 28:3453–86. doi:10.3934/dcdsb.2022226

CrossRef Full Text | Google Scholar

14. So WH, Waltman P. A nonlinear boundary value problem arising from competition in the chemostat. Appl Math Comput (1989) 32:169–83. doi:10.1016/0096-3003(89)90092-1

CrossRef Full Text | Google Scholar

15. Beddington JR. Mutual interference between parasites or predators and its effect on searching efficiency. J Anim Ecol (1975) 44:331–40. doi:10.2307/3866

CrossRef Full Text | Google Scholar

16. DeAngelis DL, Goldstein RA, O’Neill RV. A model for trophic interaction. Ecology (1975) 56:661–892.

CrossRef Full Text | Google Scholar

17. Harrision GW. Comparing predator-prey models to Luckinbill’s experiment with didinium and paramecium. Ecology (1995) 76:357–69.

Google Scholar

18. Jiang HL, Wu JH, Wang LJ, Guo GH. Qualitative analysis for a competition model with B-D functional response and numerical simulation. Numer Methods Partial Differential Equations (2014) 30:1575–94. doi:10.1002/num.21848

CrossRef Full Text | Google Scholar

19. Zhang GH, Wang XL. Extinction and coexistence of species for a diffusive intraguild predation model with B-D functional response. Discrete Contin Dyn Syst Ser B (2018) 23:3755–86. doi:10.3934/dcdsb.2018076

CrossRef Full Text | Google Scholar

20. Meng Q, Yang L. Steady state in a cross-diffusion predator-prey model with the Beddington-DeAngelis functional response. Nonlinear Anal Real World Appl (2019) 45:401–13. doi:10.1016/j.nonrwa.2018.07.012

CrossRef Full Text | Google Scholar

21. He X, Zheng SN. Protection zone in a diffusive predator-prey model with Beddington-DeAngelis functional response. J Math Biol (2017) 75:239–57. doi:10.1007/s00285-016-1082-5

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Feng XZ, Sun C, Yang WB, Li CT. Dynamics of a predator-prey model with nonlinear growth rate and B-D functional response. Nonlinear Anal Real World Appl (2023) 70:103766. in press. doi:10.1016/j.nonrwa.2022.103766

CrossRef Full Text | Google Scholar

23. Wang YE, Wu JH, Guo GH. Coexistence and stability of an unstirred chemostat model with Beddington-DeAngelis function. Comput Math Appl (2010) 60:2497–507. doi:10.1016/j.camwa.2010.08.057

CrossRef Full Text | Google Scholar

24. Nie H, Wu JH. Coexistence of an unstirred chemostat model with Beddington-DeAngelis functional response and inhibitor. Nonlinear Anal Real World Appl (2010) 11:3639–52. doi:10.1016/j.nonrwa.2010.01.010

CrossRef Full Text | Google Scholar

25. Feng XZ, Sun SP, Zhang TQ, An XM. The effect of parameters on positive solutions and asymptotic behavior of an unstirred chemostat model with B-D functional response. Adv Difference Equ (2018) 23:181–204. doi:10.1186/s13662-018-1587-x

CrossRef Full Text | Google Scholar

26. Zhang S, Tan D, Chen LS. Chaotic behavior of a chemostat model with Beddington-DeAngelis functional response and periodically impulsive invasion. Chaos Solitons Fractals (2006) 29:474–82. doi:10.1016/j.chaos.2005.08.026

CrossRef Full Text | Google Scholar

27. Zhou X, Song X, Shi X. Analysis of competitive chemostat models with the Beddington-DeAngelis functional response and impulsive effect. Appl Math Model (2007) 31:2299–312. doi:10.1016/j.apm.2006.08.010

CrossRef Full Text | Google Scholar

28. Hirsch MW, Smith HL, Zhao XQ. Chain transitivity, attractivity, and strong repellors for semidynamical systems. J Dynam Differential Equations (2001) 13:107–31. doi:10.1023/a:1009044515567

CrossRef Full Text | Google Scholar

29. Cantrell R, Cosner C. Spatial ecology via reaction-diffusion equations. Chester, UK: John Wiley and Sons Ltd (2003).

Google Scholar

30. Lou Y, Nie H, Wang YE. Coexistence and bistability of a competition model in open advective environments. Math Biosci (2018) 306:10–9. doi:10.1016/j.mbs.2018.09.013

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Crandall MG, Rabinowitz PH. Bifurcation from simple eigenvalues. J Funct Anal (1971) 8:321–40. doi:10.1016/0022-1236(71)90015-2

CrossRef Full Text | Google Scholar

32. Protter MH, Weinberger HF. Maximum principles in differential equations. New York: Springer-Verlag (1984).

Google Scholar

33. Shi JP. Persistence and bifurcation of degenerate solutions. J Funct Anal (1999) 169:494–531. doi:10.1006/jfan.1999.3483

CrossRef Full Text | Google Scholar

34. Shi JP, Wang XF. On global bifurcation for quasilinear elliptic systems on bounded domains. J Differential Equations (2009) 246:2788–812. doi:10.1016/j.jde.2008.09.009

CrossRef Full Text | Google Scholar

35. L´opez-G´omez J. Global bifurcation for Fredholm operators. Rend Istit Mat Univ Trieste (2016) 48:539–64. doi:10.13137/2464-8728/13172

CrossRef Full Text | Google Scholar