Effects of anisotropic diffusion in a two-dimensional unstirred chemostat

1 Introduction

The chemostat is a widely employed laboratory apparatus for continuous microbial cultures, with significant applications in fields such as fermentation engineering and wastewater treatment. Moreover, it is often used as a simplified model to study more intricate microbial ecosystems, including ponds and lakes. Consequently, various mathematical models have been developed to characterize the interactions between nutrients and microorganisms in chemostats, integrating factors such as species growth, diffusion, and interspecies interactions; see, e.g., [1].

Initial formulations of chemostat dynamics were based on the premise of complete mixing within the culture vessel, thereby maintaining uniform nutrient dispersion. Under this assumption, the system is commonly modeled by ordinary differential equations that describe spatially homogeneous concentrations of nutrients and microbial populations; see, e.g., [2]. However, substantial discrepancies remain between these idealized chemostat models and the complexities of real natural environments. To bridge this gap and enhance ecological realism, subsequent research has led to substantial modifications of classical chemostat models, yielding a wide array of generalized formulations. Prominent examples encompass models incorporating periodic temporal variations [3], [4], [5], delayed responses [6], and the presence of inhibitors [7], [8].

Diffusion is a widespread natural process, which has led to considerable interest in the analysis of unstirred chemostat models. Hsu and Waltman [9] introduced random diffusion into the chemostat model and derived results concerning the coexistence of two competing populations. Their findings revealed that diffusion plays a crucial role in population coexistence. In the context of unstirred chemostat models, Wu [10] studied the existence of coexistence states using global bifurcation theory. Subsequently, Wu and collaborators extended this line of research to unstirred chemostat models incorporating inhibitors [11], as well as models involving two resources [12], [13]. Furthermore, it is also noted in ref. [9] that assuming nutrients and microorganisms have identical diffusion rates is biologically unrealistic. Previous studies have presented various results for different random diffusion rates; see, e.g., [14], [15], [16], [17].

However, population dispersal is often non-random due to heterogeneity in resource distribution, variability in topography and climate, and anthropogenic modifications, leading to directional movement biases. Anisotropic diffusion is one such form that captures these complexities. Bouin et al. [18] proposed a scenario where the population moves horizontally with probability θ 2 and vertically with probability 1 − θ 2 , where θ ∈ [0, 1]. Moreover, they showed that in a model of two competing populations in two-dimensional heterogeneous environments, the optimal dispersal strategy for both populations is to migrate preferentially in the direction exhibiting lower variability in resource distribution. For more advanced theoretical studies concerning anisotropy or other diffusion mechanisms, the reader is referred to refs. [19], [20], [21].

Biologically, anisotropic diffusion may offer more adaptive possibilities for populations, enabling them to better align with available resources. This can enhance the likelihood of coexistence. Exploring the evolution of anisotropic diffusion in more biologically plausible environments would be both intriguing and significant.

Motivated by the works of refs. [9], [10], [18], we investigate in this study an unstirred chemostat model incorporating anisotropic diffusion.

subject to the boundary conditions

and the initial conditions

Here, Δ ≔ ∂ 2 ∂ x 2 + ∂ 2 ∂ y 2 is the Laplace operator and ∇ ⋅ ν is the normal derivative on the boundary, where ν = (ν _x , ν _y ) denotes the outward unit normal vector. S(x, y, t) represents the nutrient concentration, while u(x, y, t) and v(x, y, t) denote the population densities of two competitors within the unstirred chemostat. Ω is a strictly convex, open, bounded subset of R 2 with a smooth boundary ∂Ω. The function D ( θ ) ≔ D ̲ + ( D ̄ − D ̲ ) θ for any θ ∈ [0, 1], where 0 < D ̲ < D ̄ to avoid the degeneracy [18]. In particular, we assume that

The functions

where m _i , a _i > 0. Here, the constants m _i > 0, i = 1, 2, represent the maximal growth rates, and the constants a _i > 0, i = 1, 2, denote the Michaelis-Menten constants. The nutrients are supplied at the rate of H(x, y) at position (x, y), where H(x, y) ≥ 0. The mixture of nutrients and microorganisms is withdrawn at a spatially dependent rate b(x, y) at location (x, y), leading to the Robin boundary conditions. The functions b(x, y), H(x, y) ∈ C^2+α(∂Ω) and b(x, y), H(x, y) ≥, ≢0 on ∂Ω. We assume that ∂Ω is partitioned into two non-empty, pairwise disjoint subsets Γ _i , i = 1, 2, with ∂Ω = Γ₁ ∩ Γ₂, where Γ₁ ≔ {(x, y) ∈ ∂Ω : b(x, y) = 0} and H(x, y) ≥, ≢0 on Γ₁. The boundary conditions (1.2) are intuitive and appropriate for this type of equations; see [9], [16], [22] for the derivation and explanation. Moreover, the initial values ( S 0 , u 0 , v 0 ) ∈ [ C ( Ω ̄ ) ] 3 . For simplicity, the parameters not explicitly specified are treated as constants throughout this paper.

If the initial condition v₀(x, y) ≡ 0, formally equivalent to setting v(x, y, t) ≡ 0 in system (1.1)–(1.3), then this system can be reduced to

Similarly, we obtain another single-species model

Furthermore, it is well known that in the absence of species in this chemostat, the nutrient concentration S(x, y, t) approaches a unique equilibrium, denoted by z; see Lemma 2.5.

Inspired by the seminal contribution of Bouin et al. [18], which established the concavity property of the principal eigenvalue for a certain class of elliptic operators with respect to the anisotropic diffusion parameter (denoted by either p or q), we significantly broaden this result. In particular, we establish the concavity property for a substantially wider class of elliptic eigenvalue problems encompassing more general functional response terms and boundary conditions in reaction-diffusion equations; see Lemma 2.1 for the precise statement.

By introducing anisotropic diffusion, where random diffusion can be seen as a special case when p = q = 1 2 (e.g., [9]), we face several technical challenges. One such challenges is the failure of the conservation law in the chemostat, which prevents us from simplifying the models. As a result, we are compelled to work with higher-dimensional operators. Another challenge arises from the necessity of simultaneously considering both the average diffusion rate, 1 2 ( D ̲ + D ̄ ) (i.e., D₀), and the diffusion strategies, p and q, within these models.

Interestingly, our findings reveal a more intricate structure than that observed in the case of isotropic diffusion. To illustrate this, we examine the steady states of the single-species model (1.4). In the isotropic diffusion setting, the corresponding results (see, e.g., [9] or Lemma 2.7) are represented by the slice p = 1 2 in the schematic bifurcation diagrams shown in Figure 1. From this figure, we observe that regardless of whether positive steady states exist in the isotropic case, variations in the anisotropic diffusion parameter p can induce the emergence or disappearance of such steady states.

Figure 1:

Schematic bifurcation diagrams of (S, u) for system (3.8); see Theorems 3.7–3.9.

We derive several key results concerning the dynamics of systems (1.1)–(1.5) are summarized briefly as follows.

1. There exists a unique connected open interval I _p ⊂ [0, 1] such that the steady state (z, 0) of system (1.4) is globally asymptotically stable for p ∈ I _p (see Theorem 3.3), and system (1.4) is uniformly persistent for p ∈ [ 0,1 ] \ I p ̄ (see Theorem 3.4). Moreover, system (1.4) admits at least one positive steady state for p ∈ [ 0,1 ] \ I p ̄ , and the positive steady state is locally asymptotically stable when the parameter p is near the outer boundary of the interval I _p (see Theorems 3.7–3.9 and Remarks there).

2. Similarly, there exists a unique connected open interval I _q ⊂ [0, 1] such that the steady state (z, 0) of system (1.5) is globally asymptotically stable for p ∈ I _q , and system (1.5) is uniformly persistent for q ∈ [ 0,1 ] \ I q ̄ (see Theorem 3.5). Moreover, system (1.5) admits at least one positive steady state for q ∈ [ 0,1 ] \ I q ̄ , and the positive steady state is locally asymptotically stable when the parameter q is near the outer boundary of the interval I _q (see Theorem 3.11 and Remark there).

3. There exists a unique connected region Σ₀ ⊂ [0, 1]² such that the steady state (z, 0, 0) of system (1.1)–(1.3) is globally asymptotically stable for (p, q) ∈ Σ₀ (see Theorem 4.2).

4. Competitive exclusion between the two species u, v in system (1.1)–(1.3) may occur if the average rate of anisotropic diffusion 1 2 ( D ̲ + D ̄ ) (i.e., D₀) is relatively low (see Theorem 4.3), where (p, q) ∈ [0, 1]²\Σ₀.

5. There may exist a region Σ₁ ⊂ [0, 1]² located on the boundary of [0, 1]² such that system (1.1)–(1.3) is uniformly persistent for (p, q) ∈ Σ₁ (see Theorem 4.4). Moreover, we derive the conditions for the linear stability of several special solutions (see Theorem 5.1). Additionally, we conclude that under certain additional conditions, system (1.1)–(1.3) admits at least one positive steady state when (p, q) lies within a specific region of Σ₁ (see Theorems 5.8 and 5.10), and the positive steady state is locally asymptotically stable when (p, q) is located in a very small region of Σ₁ (see Lemmas 5.5 and 5.6).

6. Moreover, we investigate the impact of anisotropic diffusion on system (1.1)–(1.3) using numerical simulations, with the goal of providing insights and support (see Section 6).

Based on the above analysis of the models, we found that it is more advantageous for the species to disperse in a single direction. This is because our conclusion shows that the diffusion strategies leading to extinction or competitive exclusion are typically positioned in the middle of their feasible area, exhibiting no distinct directional bias. In contrast, anisotropic diffusion can significantly enhance the potential for species coexistence, which is consistent with the conclusions drawn in ref. [18].

The organization of this paper is as follows. In Section 2, we present several findings that will be utilized in the subsequent analysis. Section 3 analyzes the dynamical behavior of the single-species models (1.4) and (1.5), and examines the existence and local stability of their positive steady states. The aim of Section 4 is to investigate the dynamic behavior of the two-species model (1.1)–(1.3). In Section 5, we further study the structure of coexistence steady states of system (1.1)–(1.3) using bifurcation theory. Furthermore, we discuss the local stability of some coexistence steady states using perturbation theory. Section 6 presents numerical simulations examining the effects of dispersal strategies. Finally, a brief discussion is provided in Section 7.

2 Preliminaries

In this section, we present some preliminary results that will be useful for the subsequent analysis.

Consider the following linear elliptic eigenvalue problem

where h ∈ L^∞(Ω). This problem admits a unique principal eigenvalue [23], denoted by λ₁(p, h), with a positive eigenfunction φ₁(p, h). The corresponding eigenfunction φ₁(p, h) is unique up to a scalar multiple. We can normalize the positive eigenfunction φ₁(p, h) by enforcing the condition ∫ Ω φ 1 2 = 1 . The variational characterization of λ₁(p, h) could be given by

Furthermore, if φ₁ ∈ H¹(Ω), then in fact φ₁ ∈ W^2,r(Ω) for any r ∈ (1, ∞); see, e.g., [24], Theorem 2.1 and Remarks there].

It is well known that λ₁(p, h) varies smoothly with p. The implicit function theorem further implies that the corresponding principal eigenfuntion φ₁(p, h) > 0 on Ω ̄ also depends continuously and differentiably on p; see, e.g., [24]. For convenience, we write ∂ λ 1 ( p , h ) ∂ p and ∂ 2 λ 1 ( p , h ) ∂ p 2 as λ 1 ′ and λ 1 ′ ′ , respectively. Similarly, we denote ∂ φ 1 ( p , h ) ∂ p by φ 1 ′ . In particular, when h satisfies the following assumption, the principal eigenvalue λ₁(p, h) exhibits more special properties. Recalling that Γ₁ = {(x, y) ∈ ∂Ω : b(x, y) = 0} and H(x, y) ≥, ≢0 on Γ₁, we introduce the following assumption.

1. There does not exist an non-empty subset Ω Γ 1 ⊂ Ω , where Γ 1 ⊂ Ω Γ 1 ̄ , such that h is constant on Ω Γ 1 ̄ .

Motivated by ref. [18], Lemma 4.4], we have the following results.

For any given h ∈ L^∞(Ω), we define a function

for p ∈ [0, 1]. Based on the proof of Lemma 2.1, it is straightforward to observe that F′(p; h) ≤ 0 in (0, 1). This function is closely aligned with the formulation in ref. [18], Theorem 2.1] and underpins the core arguments of the subsequent investigation. It quantifies the difference between the resource variations in horizontal and vertical directions. For F > 0, the resource distribution exhibits greater spatial variations in the horizontal direction, while for F < 0, the variation is more pronounced in the vertical direction.

Next, we summarize several results on λ₁(p, h) obtained in prior studies, while the proofs are not repeated here.

The monotonicity of the principal eigenvalue λ 1 1 2 , h with respect to D₀ can be established through its variational formulation; see, e.g., equation (2.2). By employing techniques similar to those used in the proofs of ref. [25], Propositions 1.3.16 and 1.3.19], we can prove the asymptotic behavior of λ 1 1 2 , h in the limits D₀ → 0+ and D₀ → +∞, respectively.

Then we restate some conclusions for certain degenerate cases without providing the proof.

For the sake of simplicity, unless otherwise specified, we shall assume that the following assumptions hold throughout the remainder of this paper:

If p = 1 2 , then system (1.4) becomes

Let D 1 * ≔ D * f 1 ( z ) , where D* is defined in Lemma 2.3. We have the following results.

Recall that λ₁(p, f₁(z)) is the principal eigenvalue of problem

with eigenfunction φ₁(p, f₁(z)), and F p ; f 1 ( z ) = ∫ Ω φ 1 p , f 1 ( z ) x 2 − φ 1 p , f 1 ( z ) y 2 . It follows that once the parameters p, m₁, a₁, b(x, y), and the domain Ω are specified, the value of F(p; f₁(z)) is uniquely determined. For convenience, we write

By Lemmas 2.1 and 2.3, we obtain the following results. The proofs are similar to that of ref. [18], Theorem 5.2], thus we omit them here.

For convenience, we summarize the conclusions on the sign of λ₁(p, f₁(z)) in Table 1. Here D₀, D₁ > 0. The symbol “*” indicates that no specific requirements are imposed on the related parameters. In this table, the set [ 0,1 ] \ p 0 * , p 1 * is empty when p 0 * = 0 and p 1 * = 1 , and the set p − * , p + * is empty when p − * = p + * or when p − * , p + * do not exist. All other sets are non-empty.

Set D 2 * ≔ D * f 2 ( z ) and F 2 ≔ F 1 2 ; f 2 ( z ) . Similarly, Table 2 can be derived, which provides essential structural insights into the analysis of system (1.1)–(1.3).

3 Single-species model

To gain insight into the population dynamics in the absence of interspecific interactions, we examine the single-species models (1.4) and (1.5), aiming to determine the conditions for population persistence in the unstirred chemostat.