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 S , u = S k 1 + S + β 1 u , f 2 S , v = S k 2 + S + β 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?

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.

In this paper, we consider following the unstirred chemostat system with the B.–D. functional response: S t = d S x x − m 1 u f 1 S , u − m 2 v f 2 S , v , x ∈ 0,1 , t > 0 , u t = d u x x + m 1 u f 1 S , u , x ∈ 0,1 , t > 0 , v t = d v x x + m 2 v f 2 S , v , x ∈ 0,1 , t > 0 , S x 0 , t = − S 0 , S x 1 , t + γ S 1 , t = 0 , t > 0 , u x 0 , t = u x 1 , t + γ u 1 , t = 0 , t > 0 , v x 0 , t = v x 1 , t + γ v 1 , t = 0 , t > 0 , S x , 0 = S 0 x ≥ 0 , x ∈ 0,1 , u x , 0 = u 0 x ≥ , ≢ 0 , v x , 0 = v 0 x ≥ , ≢ 0 , x ∈ 0,1 , (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]: f ̂ 1 S , u = f 1 S , u , S ≥ 0 , u ≥ 0 , 0 , o t h e r s , f ̂ 2 S , v = f 2 S , v , S ≥ 0 , v ≥ 0 , 0 , o t h e r s .

For convenience, we still denote f ̂ i ( S , u ) ( i = 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 ϕ ( x ) = S 0 ( 1 + γ γ − x ) (see [14], Lemma 2.1); that is, lim t → ∞ S x , t + u x , t + v x , t = ϕ x u n i f o r m l y f o r x ∈ 0,1 .

Hence, we apply the classical internal chain transitive theory [[28], Lemma 2.1] to reduce system (2.1) into the following limiting system: u t = d u x x + m 1 f 1 ϕ x − u − v , u u , x ∈ 0,1 , t > 0 , v t = d v x x + m 2 f 2 ϕ x − u − v , v v , x ∈ 0,1 , t > 0 , u x 0 , t = u x 1 , t + γ u 1 , t = 0 , v x 0 , t = v x 1 , t + γ v 1 , t = 0 , t > 0 , u x , 0 = u 0 x ≥ , ≢ 0 , v x , 0 = v 0 x ≥ , ≢ 0 , u 0 x + v 0 x ≤ ϕ x x ∈ 0,1 . (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: d u x x + m 1 f 1 ϕ x − u − v , u u = 0 , x ∈ 0,1 , d v x x + m 2 f 2 ϕ x − u − v , v v = 0 , x ∈ 0,1 , u x 0 = u x 1 + γ u 1 = 0 , v x 0 = v x 1 + γ v 1 = 0 . (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.

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

We first consider a linear eigenvalue problem d φ x x + q x φ = μ φ , x ∈ 0,1 , φ x 0 = ϕ x 1 + γ φ 1 = 0 , (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 max x ∈ [ 0,1 ] φ 1 = 1 . Furthermore, by the variation characterization of the principal eigenvalue [29], we have μ 1 q x = − inf φ ∈ H 1 0,1 , φ ≠ 0 d ∫ 0 1 φ x 2 d x + d γ φ 2 1 − ∫ 0 1 q x φ 2 d x ∫ 0 1 φ 2 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, + ∞).

We consider the following single-species model: ω t = d ω x x + m f ϕ − ω , ω ω , x ∈ 0,1 , t > 0 , ω x 0 , t = ω x 1 , t + γ ω 1 , t = 0 , t > 0 , ω x , 0 = ω 0 x ≥ , ≢ 0 , x ∈ 0,1 , (3.3) where d, m > 0 are constants and f ( ϕ − ω , ω ) = ϕ − ω k + ϕ + ( β − 1 ) ω . For fixed d, k > 0, let μ ₁(mf(ϕ, 0)) be the principal eigenvalue of d φ x x + m f ϕ , 0 φ = μ φ , x ∈ 0,1 , t > 0 , φ x 0 = φ x 1 + γ φ 1 = 0 , where f ( ϕ , 0 ) = ϕ k + ϕ . Then, we can conclude from Lemma 3.1(iii) that μ ₁(mf(ϕ, 0)) is strictly increasing with respect to m in (0, + ∞). Moreover, lim m → 0 + μ 1 ( m f ( ϕ , 0 ) ) = μ 1 ( 0 ) < 0 by Lemma 3.1(iii), and lim m → + ∞ μ 1 ( m f ( ϕ , 0 ) ) = + ∞ by Eq. 3.2. For fixed d, k > 0, there exists a unique critical value m ^⋆ such that μ 1 m f ϕ , 0 < 0 i f 0 < m < m ∗ , μ 1 m f ϕ , 0 = 0 i f m = m ∗ , μ 1 m f ϕ , 0 > 0 i f m > m ∗ . (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 lim t → ∞ ω ( x , t ) = ω ∗ ( ⋅ ; m , β ) 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,

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 d P m ″ + m f ϕ − ω ∗ , ω ∗ − k + β ϕ ω ∗ k + ϕ + β − 1 ω ∗ 2 P m = − f ϕ − ω ∗ , ω ∗ ω ∗ , x ∈ 0,1 , P m ′ 0 = P m ′ 1 + γ P m 1 = 0 . (3.7)

We define X = { ψ ∈ C 2 [ 0,1 ] : ψ ′ ( 0 ) = ψ ′ ( 1 ) + γ ψ ( 1 ) = 0 } and denote L 1 ( ψ ) = d ψ ′ ′ + m [ f ( ϕ − ω ∗ , ω ∗ ) − ( k + β ϕ ) ω ∗ [ k + ϕ + ( β − 1 ) ω ∗ ] 2 ] ψ . It is easy to see that μ ₁(L ₁) < μ ₁(mf(ϕ − ω _∗, ω _∗)) = 0 by Lemma 3.1(iii). Noting that P m ( x ) ∈ 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 ω ∗ ′ ′ = − m d f ( ϕ − ω ∗ , ω ∗ ) ω ∗ < 0 , ω ∗ ′ is decreasing on [0,1]. Note that the boundary conditions ω ∗ ′ ( 0 ) = 0 , ω ∗ ′ ( 1 ) = − γ ω ∗ ( 1 ) . Then, ω ∗ ′ ( 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 lim m → m ∗ + ω ∗ ⋅ ; m , β = ω 1 , lim m → + ∞ ω ∗ ⋅ ; m , β = ω 2 i n C 0,1 .

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, lim m → ( m ∗ ) + ω ∗ ( ⋅ ; m , β ) = ω 1 weakly in W ², p (0, 1) and the convergence also holds in C ¹[0, 1] by the Sobolev embedding theorem. Then, ω ₁ satisfies d ω 1 ′ ′ + m ∗ f ϕ − ω 1 , ω 1 ω 1 = 0 , x ∈ 0,1 , ω 1 ′ 0 = ω 1 ′ 1 + γ ω 1 1 = 0 . (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 d ω ∗ ′ ′ + m f ϕ − ω ∗ , ω ∗ ω ∗ = 0 , x ∈ 0,1 , ω ∗ ′ 0 = ω ∗ ′ 1 + γ ω ∗ 1 = 0 . (3.9)

Dividing the first equation of Eq. 3.9 by mω _∗ and integrating over (0,1), d m ∫ 0 1 | ω ∗ ′ | 2 ω ∗ 2 d x − d γ m + ∫ 0 1 f ϕ − ω ∗ , ω ∗ d x = 0 , which implies 0 ≤ ∫ 0 1 f ϕ − ω ∗ , ω ∗ d x ≤ d γ m .

Note lim m → + ∞ ω ∗ ( ⋅ ; m , β ) = ω 2 in C[0, 1]. Taking m → +∞, we have ∫ 0 1 f ( ϕ − ω 2 , ω 2 ) 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 lim β → + ∞ ω ∗ ( ⋅ ; m , β ) = 0 u n i f o r m l y o n [ 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.

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 ( ϕ , 0 ) = ϕ k i + ϕ ( i = 1,2 ) , by the similar arguments as in (3.4), we can define m i ∗ such that μ 1 m i f i ϕ , 0 < 0 i f 0 < m i < m i ∗ , μ 1 m i f i ϕ , 0 = 0 i f m i = m i ∗ , μ 1 m i f i ϕ , 0 > 0 i f m i > m i ∗ . (4.1)

Clearly, m i ∗ ( i = 1,2 ) are dependent on k _i but not on β _i .

Proposition 4.1

For fixed d > 0, the following statements hold: ( i ) k 2 k 1 < m 2 ∗ m 1 ∗ < 1 if k 1 > k 2 > 0 ; ( ii ) 1 < m 2 ∗ m 1 ∗ < k 2 k 1 if k 2 > k 1 > 0 ; ( iii ) m 1 ∗ = m 2 ∗ if k 1 = k 2 > 0 .

Proof

(i) Since μ 1 m 1 ∗ ϕ k 1 + ϕ = 0 a n d μ 1 m 2 ∗ ϕ k 2 + ϕ = 0 , (4.2) it is easy to check that m 1 ∗ > m 2 ∗ by k ₁ > k ₂ > 0 and Lemma 3.1(iii). Note that Eq. 4.2) is equivalent to μ 1 ( m 1 ∗ ϕ / k 1 1 + ϕ / k 1 ) = 0 and μ 1 ( m 2 ∗ ϕ / k 2 1 + ϕ / k 2 ) = 0 . Since ϕ 1 + ϕ / k 1 > ϕ 1 + ϕ / k 2 based on k ₁ > k ₂ > 0, we have m 1 ∗ k 1 < m 2 ∗ k 2 ; that is, k 2 k 1 < m 2 ∗ m 1 ∗ . Therefore, k 2 k 1 < m 2 ∗ m 1 ∗ < 1 . (ii) can be obtained similarly. For (iii), it is easy to obtain m 1 ∗ = m 2 ∗ by the fact of k ₁ = k ₂ > 0 and μ 1 ( m i ∗ f i ( ϕ , 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 ∗ ; semi-trivial solution (0, ω _∗(⋅; m ₂, β ₂)) exists if and only if m 2 > m 2 ∗ . For convenience, we denote u ̂ = ω ∗ ( ⋅ ; m 1 , β 1 ) , v ̂ = ω ∗ ( ⋅ ; m 2 , β 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 ) ≤ u ̂ and 0 < v ( x ) ≤ v ̂ 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 ( u ̂ , 0 ) and ( 0 , v ̂ ) . First, the linearized operator of system (2.3) at ( u ̂ , 0 ) is given by d φ x x + m 1 f 1 ϕ − u ̂ , u ̂ + u ̂ f 1 u ϕ − u ̂ , u ̂ φ + m 1 u ̂ f 1 v ϕ − u ̂ , u ̂ ψ = μ φ , x ∈ 0,1 , d ψ x x + m 2 f 2 ϕ − u ̂ , 0 ψ = μ ψ , x ∈ 0,1 , φ x 0 = φ x 1 + γ φ 1 = 0 , ψ x 0 = ψ x 1 + γ ψ 1 = 0 , (4.3) where f 1 u ( ϕ − u ̂ , u ̂ ) = − k 1 + β 1 ϕ [ k 1 + ϕ + ( β 1 − 1 ) u ̂ ] 2 < 0 , f 1 v ( ϕ − u ̂ , u ̂ ) = − k 1 + β 1 u ̂ [ k 1 + ϕ + ( β 1 − 1 ) u ̂ ] 2 < 0 . By the Riesz–Schauder theory, the eigenvalues of Eq. 4.3 consist of the eigenvalues of the following two operators: B 1 m 1 = d d 2 d x + m 1 f 1 ϕ − u ̂ , u ̂ + u ̂ f 1 u ϕ − u ̂ , u ̂ , B 2 m 2 = d d 2 d x + m 2 f 2 ϕ − u ̂ , 0 , (4.4) subject to the corresponding boundary conditions. It follows from Lemma 3.1 (iii) that μ 1 ( m 1 ( f 1 ( ϕ − u ̂ , u ̂ ) + u ̂ f 1 u ( ϕ − u ̂ , u ̂ ) ) ) < μ 1 ( m 1 f 1 ( ϕ − u ̂ , u ̂ ) ) . Moreover, μ 1 ( m 1 f 1 ( ϕ − u ̂ , u ̂ ) ) = 0 with eigenfunction φ 1 ( m 1 f 1 ( ϕ − u ̂ , u ̂ ) ) = u ̂ , which implies μ 1 ( m 1 f 1 ( ϕ − u ̂ , u ̂ ) + u ̂ f 1 u ( ϕ − u ̂ , u ̂ ) ) < 0 . Hence, the stability of ( u ̂ , 0 ) is determined by the sign of the principal eigenvalue μ 1 ( m 2 f 2 ( ϕ − u ̂ , 0 ) ) of B 2 ( m 2 ) . More precisely, ( u ̂ , 0 ) is asymptotically stable if μ 1 ( m 2 f 2 ( ϕ − u ̂ , 0 ) ) < 0 , while ( u ̂ , 0 ) is unstable if μ 1 ( m 2 f 2 ( ϕ − u ̂ , 0 ) ) > 0 .

The linearized operator of system (2.3) at ( 0 , v ̂ ) is given by d φ x x + m 1 f 1 ϕ − v ̂ , 0 φ = μ φ , x ∈ 0,1 , d ψ x x + m 2 f 2 ϕ − v ̂ , v ̂ + v ̂ f 2 v ϕ − v ̂ , v ̂ ψ + m 2 v ̂ f 2 u ϕ − v ̂ , v ̂ φ = μ ψ , x ∈ 0,1 , φ x 0 = φ x 1 + γ φ 1 = 0 , ψ x 0 = ψ x 1 + γ ψ 1 = 0 , (4.5) where f 2 u ( ϕ − v ̂ , v ̂ ) = − k 2 + β 2 v ̂ [ k 2 + ϕ + ( β 2 − 1 ) v ̂ ] 2 < 0 , f 2 v ( ϕ − v ̂ , v ̂ ) = − k 2 + β 2 ϕ [ k 2 + ϕ + ( β 2 − 1 ) v ̂ ] 2 < 0 . We denote B 3 m 2 = d d 2 d x + m 1 f 1 ϕ − v ̂ , 0 , B 4 m 2 = d d 2 d x + m 2 f 2 ϕ − v ̂ , v ̂ + v ̂ f 2 v ϕ − v ̂ , v ̂ . (4.6)

Similarly, ( 0 , v ̂ ) is asymptotically stable if μ 1 ( m 1 f 1 ( ϕ − v ̂ , 0 ) ) < 0 , while ( 0 , v ̂ ) is unstable if μ 1 ( m 1 f 1 ( ϕ − 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 ≤ m 1 ⋆ and m 2 ≤ m 2 ⋆ , then

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 ∗ , m 2 ≥ max { k 2 k 1 , 1 } m 1 , and β ₁ > 0 fixed, ( u ̂ , 0 ) is unstable.

Note that u ̂ satisfies d u ̂ x x + m 1 f 1 ϕ − u ̂ , u ̂ u ̂ = 0 , x ∈ 0,1 , u ̂ x 0 = u ̂ x 1 + γ u ̂ 1 = 0 , (4.12) which implies μ 1 ( m 1 f 1 ( ϕ − u ̂ , u ̂ ) ) = 0 . We recall that ( u ̂ , 0 ) is asymptotically stable if μ 1 ( m 2 f 2 ( ϕ − u ̂ , 0 ) ) < 0 and it is unstable if μ 1 ( m 2 f 2 ( ϕ − u ̂ , 0 ) ) > 0 . Then, we conclude from m 2 ≥ max { k 2 k 1 , 1 } m 1 and 0 < u ̂ < ϕ on [0,1] that m 1 f 1 ϕ − u ̂ , u ̂ − m 2 f 2 ϕ − u ̂ , 0 < m 1 f 1 ϕ − u ̂ , 0 − m 2 f 2 ϕ − u ̂ , 0 = m 1 k 2 − m 2 k 1 + m 1 − m 2 ϕ − u ̂ k 1 + ϕ − u ̂ k 2 + ϕ − u ̂ ϕ − u ̂ ≤ 0 , , which means that μ 1 ( m 2 f 2 ( ϕ − u ̂ , 0 ) ) > μ 1 ( m 1 f 1 ( ϕ − u ̂ , u ̂ ) ) = 0 by Lemma 3.2(iii). That is, ( u ̂ , 0 ) is unstable.

Claim 2. (1) For m 1 > m 1 ∗ , m 2 > max { m 2 ∗ , max { k 2 k 1 , 1 } m 1 } , and β ₁ > 0 fixed, ( 0 , v ̂ ) is asymptotically stable when 0 < β 2 ≤ β 2 0 , where β 2 0 is defined by Eq. 4.10.

(2) There exists a unique β 2 ∗ > β 2 0 such that ( 0 , v ̂ ) is asymptotically stable when 0 < β 2 < β 2 ∗ , and ( 0 , v ̂ ) is unstable when β 2 > β 2 ∗ .

For (1), we recall that ( 0 , v ̂ ) is asymptotically stable if μ 1 ( m 1 f 1 ( ϕ − v ̂ , 0 ) ) < 0 and unstable if μ 1 ( m 1 f 1 ( ϕ − v ̂ , 0 ) ) > 0 . Similarly, we can conclude from m 2 > max { m 2 ∗ , max { k 2 k 1 , 1 } m 1 } and 0 < v ̂ < ϕ ( x ) < S 0 ( 1 + γ ) γ on [0,1] that m 1 f 1 ϕ − v ̂ , 0 − m 2 f 2 ϕ − v ̂ , v ̂ = m 1 k 2 − m 2 k 1 + m 1 − m 2 ϕ − v ̂ + m 1 β 2 v ̂ k 2 + ϕ + β 2 − 1 v ̂ k 1 + ϕ − v ̂ ϕ − v ̂ < 0 , provided 0 < β 2 ≤ β 2 0 ≔ ( m 2 k 1 − m 1 k 2 ) γ m 1 S 0 ( 1 + γ ) . It follows from Lemma 3.2(iii) that μ 1 ( m 1 f 1 ( ϕ − v ̂ , 0 ) ) < μ 1 ( m 2 f 2 ( ϕ − v ̂ , v ̂ ) ) = 0 when 0 < β 2 ≤ β 2 0 . That is, ( 0 , v ̂ ) is asymptotically stable when 0 < β 2 ≤ β 2 0 .