Analysing the systemic implications of energy efficiency and circular economy strategies in the decarbonisation context

1.   Introduction

The two-dimensional Lotka-Volterra competition model has been widely studied ^[1,2,3,4]. Murray discussed a simple model on competition between two species and explains the biological conclusion implied by the competition exclusion principle, namely, when two species compete for the same limited resources, one of the species usually becomes extinct ^[1]. Song et al. considered the Lotka-Volterra competitive system with two delays and studied the dynamical properties of the system ^[2]. Jin et al. considered a nonautonomous two-dimensional competitive Lotka-Volterra system with impulse and gave an extension of the principle of competition exclusion ^[3]. Zhang and Chen studied the effects of linear and nonlinear diffusion of the competitive Lotka-Volterra model ^[4]. Hsu and Zhao analyzed the stability of a monotone dynamical system and obtained the global dynamics of the Lotka-Volterra two-species competition model with seasonal succession ^[5]. Wang and Zhao investigated two free boundary problems of a Lotka-Volterra type competition model in a one-dimension space to understand the asymptotic behavior of the two competing species spreading through a free boundary ^[6]. Liu and Fan proposed a new definition of permanence for stochastic population models, which overcame some limitations and deficiency of the existing ones ^[7]. Ren and Liu considered the two-species chemotaxis system with Lotka-Volterra competitive kinetics in a bounded domain with smooth boundary ^[8].

In fact, all organisms are composed of multiple chemical elements such as carbon, nitrogen and phosphorus (Loladze et al. ^[9]). Ecological stoichiometry is the balance of multiple chemical substances in ecological interactions and processes, or the study of such balances (Sterner and Elser ^[10]). From the stoichiometric point of view ^[9,10], both food quantity and quality need to be explicitly modeled in producer-consumer interactions, as it has been observed that the plant quality and quantity can greatly affect the growth rate of herbivorous grazers ^{[11,12,13,14,15,16,17]}. Loladze et al. applied stoichiometric principles to construct a two-dimensional Lotka-Volterra type model (LKE model) that links energy flow with element cycling and explains the Rosenzweig's paradox of enrichment ^[9]. Wang et al. proposed a stoichiometric producer-grazer model, which extended the LKE model by taking P into account in both the prey and the media ^[18]. The results show that as the uptake rate of P by producer becomes infinite, LKE models become the limiting case of this model. Peace et al. formulated a model to explicitly track free P in the stoichiometric knife-edge model. The model shows that the fate of the grazer population can be very sensitive to excessive nutrient concentrations ^[19].

Furthermore, the growing study of ecological stoichiometry also further facilitates more insights on mechanisms of competitive systems between the two species. For example, Tilman et al. explained the mechanism of plant competition for nutrients by giving the expression of plant competition for nutrients ^[20]. Ji and Wang incorporated stoichiometry into a competition model in chemostat settings ^[21] and considered two species competing for a single nutrient, and investigated how stoichiometry, dilution rate and concentration of phosphorus input affected the result of competition between algae species. However, for the system with both interspecific and interspecific competition in ^[1], the effect of nutrients on the growth of both species was unclear and an ecological stoichiometric model was essential.

In this study, by considering the ecological stoichiometry of nutrient elements, we formulate a new stoichiometric competition model to investigate the dynamics among two plants in section two. Here, the phosphorus in the system is divided into three pools: Phosphorus in the first plant, phosphorus in the second plant and free phosphorus. We analyze the positive invariant region of the model and derive the existence and stability conditions of the equilibria. In section three, we first study the stability of boundary equilibria in Case 1, then, we provide the solution curves of the model and verify the results of theoretical analysis. Finally, we study the bifurcation diagram of the system about $K$ and $T$, and the bistability of the system with interspecific competition is employed to investigate how total phosphorus in the system and light intensity affect the competition outcomes. Section four concludes and discusses the implications of our mathematical findings.

2.   Modeling the competition between two plants

In this section, we present how we formulate the mathematical model for the competition of two plants. Next, we derive the equilibria of the model and analyze the stability.

2.1.   Model formulation

In this work, we integrate stoichiometry into the competition model of two plants. We start with a two-species Lotka-Volterra competition model ^[1]:

In model (2.1), $x_i$ ($i = 1, 2$) describe the biomasses of the plant measured in terms of carbon. $r_1$, $r_2$, $K_1$, $K_2$, $b_{12}$ and $b_{21}$ are all positive constants, $r_i$ ($i = 1, 2$) are the linear birth rates and $K_i$ ($i = 1, 2$) are the carrying capacities. We use $b_{12}$ and $b_{21}$ to measure the competitive effect of $x_1$ on $x_2$ and $x_2$ on $x_1$, respectively. Model (2.1) only considers the energy flow in a species. However, element cycling cannot be ignored in a competition model. Hence, it is realistic to take into account the stoichiometric competition model.

We assume that plants are composed of two main elements, carbon and phosphorus, and we consider two plants, plant one and plant two. Let $T$ represent the total phosphorus in the system. $p_i$ ($i = 1, 2$) denotes the total phosphorus in plant $i$, and its unit is (mg P)/L. $b_{12}$ and $b_{21}$ are the competition coefficients between two plants. Following the main assumptions used in the LKE model ^[9], we assume

(ⅰ) The total phosphorus in the system ($T$) in the entire system is fixed.

(ⅱ) Phosphorus to carbon ratio (P:C) varies in plants, but it never falls below a minimum ratio for plant one and plant two, denoted by $q_{1}$ and $q_{2}$, respectively.

Recall that P:C in the plants should be at least $q_i$ ($i = 1, 2$) ((mg P)/(mg C)), and one obtains that the plant density cannot exceed $p_i/q_i$ ((mg P)/L). We consider the competition of plants for carbon (light) and the restriction of phosphorus on the growth of each plant. The internal phosphorus dependent growth function follows the empirically well tested Droop form, $1-\frac{q_i}{p_i/x_i}\; (i = 1, 2)$ ^[22]. Then, the growth is restricted by the light intensity ($K$) and the total phosphorus in the system ($p_1/q_1, p_2/q_2$). The model takes the form of

Here, $T-p_{1}-p_{2}$ is the free phosphorus in the environment. Let $d_{i}$ ($i = 1, 2$) be the loss rate of phosphorus in the plants, then the phosphorus element metabolized by plants is $d_{i}p_{i}$ ($i = 1, 2$). Hence, the biological significance of other parameters is explained below. Parameter values estimated by Wang et al. (2008) ^[18] are listed in Table 1. The replenishment rate functions $g_i(y)$ are increasing, nonnegative for nonnegative variables. In general, $g_i(y)$ is a bounded smooth function that satisfies the following assumptions:

To prevent the orbits from entering the origin where the system is undefined, we devise a transformation to generate a new system in ^[15]. The transformation converts variables plant C and plant P to variables plant C and plant P:C ratio. We introduce $Q_1 = p_1/x_1$ and $Q_2 = p_2/x_2$, which are the cell quota (intracellular P:C ratio) in two plants. We set the transformation

which converts (2.2) in $\Omega$ into the following system

Here, $\varPhi(\Omega)$ = $\{ (x_1, x_2, Q_1, Q_2)$: $0 < x_{1} < \min\{K, T/q_{1}\}$, $0 < x_{2} < \min\{K, T/q_{2}\}$, $Q_1 > q_1$, $Q_2 > q_2$, $Q_1x_1+Q_2x_2 < T \}$.

More details of mathematical analysis can be found in Appendix A.

3.   Numerical dynamics and the implications

In this section, we first analyze the stability of the two boundary equilibria in Case 1, and then we analyze the bifurcation of (A.1). Finally, we study the stability of the system with interspecific competition. The initial conditions in our simulations are chosen inside the biologically meaningful region.

3.1.   Stability analysis of the boundary equilibria in Case 1

In this subsection, we analyze the stability of the two boundary equilibria in Case 1. The model (A.2) has two boundary equilibria

Based on the study of Theorem 2, the first two equations of (A.2) do not contain variables $Q_1$ and $Q_2$. Here, we consider the submodel of (A.2)

We present the phase diagram of (3.1). In Figure 1, $r_1 = 0.93$, $r_2 = 1.2$ and $K = 10$. If the initial value starting from the inside will eventually tend to the straight line composed of the internal equilibrium family of (3.1), then the boundary equilibria $E_1$ and $E_2$ are unstable. However, when the initial value starts from the boundary, the solution starting from the boundary eventually tends to the boundary equilibria. Since the forward invariant domain of (3.1) does not contain boundaries, the boundary equilibria of (3.1) are unstable. Therefore, the boundary equilibria of (A.2) are unstable in $\varPhi(\Omega)$.

3.2.   Solution curve of the model

In this subsection, we validate the results of the theoretical analysis by plotting the solution curve of (A.1) and marking the corresponding situations for Case 1 to Case 4 in the solution curve. We select the initial value $x_1(0) = 0.005$, $x_2(0) = 0.005$, $Q_1(0) = 0.001$ and $Q_2(0) = 0.001$. In Figures 2–4, the blue dot represents the solution of (A.1) in Case 1, the red dot represents the solution of (A.1) in Case 2, the green dot represents the solution of (A.1) in Case 3 and the black dot represents the solution of (A.1) in Case 4.

Figure 2 shows that, initially, the biomasses of the two plants increased. When $t \in (0, 11)$, the solution of (A.1) corresponds to the solution in Case 4. When $t = 11$, the solution of (A.1) changes from Case 4 to Case 2. When $t = 12$, the solution of (A.1) changes from Case 2 to Case 1. Finally, the solution of (A.1) converges to a stable state. This result confirms the conclusion that the internal equilibrium $E_*$ in Case 1 is globally asymptotically stable.

Figure 3 shows that the biomasses of two plants increasing over time. When $t \in (0, 11)$, the solution of (A.1) corresponds to the solution in Case 4. When $t = 10$, the solution of (A.1) transitions from Case 4 to Case 3. When $t = 11$, the solution of (A.1) transitions from Case 3 to Case 1. When $t = 49$, the solution of (A.1) transitions from Case 1 to Case 2. Finally, the solution of (A.1) converges to a stable state. This result verifies the conclusion that the internal equilibrium $E_*^{1}$ in Case 2 is locally asymptotically stable.

Figure 4 shows that the biomasses of two plants increases over time. When $t \in (0, 10)$, the solution of (A.1) corresponds to the solution in Case 4. When $t = 10$, the solution of (A.1) changes from Case 4 to Case 2. When $t = 11$, the solution of (A.1) changes from Case 2 to Case 1. When $t = 49$, the solution of (A.1) changes from Case 1 to Case 3. Finally, the solution of (A.1) converges to a stable state. This result verifies the conclusion that the internal equilibrium $E_*^{2}$ in Case 3 is locally asymptotically stable.

According to Figures 2–4, it can be observed that the solution of (A.1) will appear in three stable states (Cases 1–3). However, no stable solution appears in Case 4, which is consistent with the theoretical analysis above.

3.3.   Impact of light intensity and total phosphorus in the system on competition dynamics

In this subsection, we investigate how the light intensity ($K$) and the total phosphorus in the system ($T$) affect the competition outcomes by plotting bifurcation diagrams of (A.1), with respect to $K$ and $T$. The phosphorus uptake rates are chosen as $g_i(y) = c_i y$ ($i = 1, 2$), where $y = T-p_{1}-p_{2}$. Figure 5 presents the bifurcation diagrams where the densities of both plants are plotted along the gradient of $K$. In the bifurcation diagrams, the solid lines represent stable equilibria and the dashed lines signify unstable equilibria. $BP$ means a bifurcation point, at which the existence or stability of the equilibria will be transformed. From Figure 5, one observes that when $T = 0.03$, $K \in (0, 8.61)$, the internal equilibrium $E_{*}$ is stable and other equilibria ($E_0, E_1, E_2$) are unstable. This means that these two plants can coexist, even if the light is low ($K \in (0, 8.61)$), when the total phosphorus in the system is adequate. Then, at $K = 8.61$ (point $BP$), the transcritical bifurcation occurs; that is, $E_{*}$ changes its stability from being stable to unstable, while the boundary equilibrium $E_2$ changes from being unstable to stable. When $K$ further increases ($K > 8.61$), there is only one stable equilibrium $E_2$. It indicates that only plant one can survive. The reason for this phenomenon is that the biomasses of both plants increase with the increase of light intensity ($K$), and there is not enough phosphorus in the ecosystem to support the growth of both plants. Hence, the plant with the lower phosphorus loss rate can adapt to the low total phosphorus in the system conditions and win the competition, and, thus eventually survives.

In particular, through Figure 6, we find that when the parameters of both plants are equal the model has no bifurcation and there is a unique stable equilibrium $E_*$. This means that two plants will coexist with the same phosphorus loss rate. The biomasses of plants will increase with the increase of light. Due to the total phosphorus constraints in the system, the biomasses of both plants will not continue increasing with the increase of light.

From Figure 7, we can observe that when $K = 0.25$, $T \in (0, 0.735 \times 10^{-3})$, the boundary equilibrium $E_0$ is stable and other equilibria ($E_*, E_1, E_2$) are unstable. Hence, both plants will go extinct under low total phosphorus in the system ($T \in (0, 0.735 \times 10^{-3})$). Then, at $T = 0.735 \times 10^{-3}$ (point $BP$), the transcritical bifurcation occurs, where $E_{0}$ changes from being stable to unstable and $E_2$ changes from being unstable to stable. It indicates that only plant one can survive. This is because the biomasses of both plants increase with the increase of total phosphorus in the system ($T$), and the plant with the lower phosphorus loss rates has a competitive advantage and eventually survives. When $T\approx 0.98 \times 10^{-3}$ (point $BP$), the transcritical bifurcation occurs ($E_1$ becomes unstable). When $T$ continues increasing until $T\approx 1.86 \times 10^{-3}$ (point $BP$), a transcritical bifurcation occurs, where $E_{2}$ changes its stability from being stable to unstable and $E_*$ changes from being unstable to stable. When $T$ further increases ($T > 1.86 \times 10^{-3}$), there is only one stable equilibrium $E_*$. This indicates that the two plants will coexist as adequate phosphorus. To conclude, when the light intensity is at a low level, low total phosphorus in the system ($T \in (0, 0.735 \times 10^{-3})$) will lead to the extinction of both plants, and when the total phosphorus in the system is at a moderate level ($0\leq 0.735 \times 10^{-3} \leq 0.98 \times 10^{-3}$), only the plant with the lower phosphorus loss rate will survive. When the total phosphorus in the system is sufficient, both plants can survive and coexist ($0.98 \times 10^{-3} \leq T \leq 1.86 \times 10^{-3}$). However, the biomasses of both plants no longer increases because of the limitation of light intensity.

3.4.   Impact of interspecific competition on competition dynamics

In this subsection, we consider the interspecific competition between the two plants and discuss the effect of competition coefficients ($b_{12}$, $b_{21}$) on the dynamics of (2.2), and we expound the details for the high, intermediate and low level of total phosphorus in the system with $T = 0.03$, $T = 0.001$ and $T = 0.0009$, respectively. We choose $g_i(y)$ as Michaelis-Menten (Monod) function; that is, $g_i(y) = \frac{c_i y}{a_i+y}$ ($i = 1, 2$) in (2.2), where $y = T-p_{1}-p_{2}$.

For $T = 0.03$, we observe from Figure 8(a) that the dynamics are similar to those of Lotka-Volterra competition model ^[1]. That is, when $b_{12} < 1$ and $b_{21} < 1$, the competition between the two plants is not intense and both plants coexist ($E_{*}$ is stable); when $b_{12} < 1$, $b_{21} > 1$, $x_2$ goes extinct, $x_1$ excludes $x_2$ during competition ($E_{1}$ is stable); when $b_{12} > 1$, $b_{21} < 1$, $x_1$ goes extinct, $x_2$ excludes $x_1$ during competition ($E_{2}$ is stable); when $b_{12} > 1$, $b_{21} > 1$, the bistability occurs between two boundary equilibria ($E_1$ and $E_2$). Which plant wins the competition depends on the initial values of the two plants.

When $T$ is decreased to 0.001 (see Figure 8(b)), the dynamics in the region $0 < b_{12} < 1$ and $0 < b_{21} < 1$ are the same as those in Figure 8, implying that the coexistent regions of these two plants are the same. However, the stable regions of $E_1$ and $E_2$ are both larger than that in Figure 8(a), and the area of the bistability of $E_1$ and $E_2$ is smaller than that with the adequate total phosphorus in the system. Note that $b_{12}$ and $b_{21}$ are the competition coefficients of the two plants, which represent the ability of competing against light, and the higher the value of the two coefficients indicates greater intensity of competition. When the total phosphorus in the system is limited and the competition intensity is high, the survival of the plants depends more on their own competition ability for light than on the initial state of the plants. For example, the species with greater competitive ability ($b_{12} > b_{21} > 1$) (plant two) will always win the competition, while the other plant becomes extinct (i.e., $E_2$ is stable).

As is shown in Figure 8(c), when $T$ is further decreased to $0.0009$, the coexistent region of the two plants is larger than that when $T = 0.03$ and $T = 0.001$, while there are also two areas of bistability ($E_1$ and $E_{*}$, $E_2$ and $E_{*}$). This suggests that when the total phosphorus in the system is insufficient, less intense competition for light is more favorable for the coexistence of the two plants. Moreover, the limitation of nutrient elements also makes the system unstable, and bistability occurs. In addition, the area of bistability of $E_1$ and $E_2$ is smaller than that in Figure 6. This further confirms our findings that when the total phosphorus in the system is limited, the survival of plants is more dependent on their competition ability for light.

3.5.   Impact of light intensity and total phosphorus in the system on competition dynamics

In this subsection, we investigate how the light intensity ($K$) and the total phosphorus in the system ($T$) affect the competition outcomes by plotting bifurcation diagrams of (A.1) with respect to $K$ and $T$. The phosphorus uptake rates are chosen as $g_i(y) = \frac{c_i y}{a_i+y}$ ($i = 1, 2$), where $y = T-p_{1}-p_{2}$. Figure 9 presents the bifurcation diagrams where the densities of both plants are plotted along the gradients of $K$ and $T$.

Figure 9 depicts the scenario that the influence coefficient of plant one on plant two is less than the influence coefficient of plant two on plant one, i.e., $b_{12} < b_{21}$. In this condition, plant one will always survive, and its biomass is the highest when both light intensity and total phosphorus are high. Plant two is quite different; the biomass of plant two tends to go extinct even when light intensity is high, and its biomass reaches a maximum when light intensity is low and total phosphorus is high (see Figure 9(a), (b)). In addition, the cell quotas of both plants will tend to be zero with the increase of light intensity (see Figure 9(c), (d)). This is because plant one has a higher competition ability for light and a lower phosphorus loss rate compared to plant two ($b_{12} < b_{21}$, $d_1 < d_2$). As a result, when light intensity increases and phosphorus becomes limited, plant one gains a competitive advantage, while plant two suffers from a light and nutrient deficiency, leading to its extinction. The cell quotas of both plant one and plant two tends to be zero due to excessive light intensity and insufficient phosphorus. Hence, even though the biomass of plants remains high under phosphorus limitation, the available "nutrition'' for their predators is insufficient, which is consistent with the paradox of enrichment.

Similarly, Figure 10 illustrates the changes in plant biomasses and quotas with $K$ and $T$ when $b_{12} > b_{21}$. From Figure 10 (a), (b), it can be observed that plant one and plant two alternate in their competitive dominance. When light intensity is low, plant two is more likely to dominate, and with increasing $T$, the dominant range of $K$ expands. On the other hand, when light intensity is high, plant one is more likely to dominate. This could be because plant two has a higher competitive ability for light compared to plant one. Therefore, when light intensity is low, plant two is capable of absorbing more light and gaining a competitive advantage. However, as light intensity increases and phosphorus becomes limiting, the lower phosphorus loss rate of plant one helps it retain more nutrients and allows it to dominate in the competition.

In summary, the light intensity ($K$), total phosphorus ($T$), interspecific competition ability ($b_{12}, \; b_{21}$) and loss rate of phosphorus ($d_1, \; d_2$) all affect the coexistence and competition of two plants. Additionally, excessive light intensity may reduce the cell quotas for both plants, which is in good agreement with the findings of the stoichiometric competition model involving two predators and one prey ^[26], where both predators were limited by the prey's phosphorus content.

4.   Discussion

In this study, we established a stoichiometric competition model between two plants that explicitly incorporated the impact of the element cycle on plants. We systematically discussed the effects of light intensity and total phosphorus in the system on the dynamics of the model.

In mathematical analysis, when the competition for lights between two plants is equal with $b_{12} = b_{21}$, we give the existence and stability analysis of the equilibria of this model in four conditions. When the growth of two plants is only limited by light intensity, the two boundary equilibria are unstable and the coexistence equilibrium is globally asymptotically stable. When the growth of plant one is limited by light intensity and plant two is limited by the total phosphorus in the system, the boundary equilibrium is unstable and the coexistence equilibrium is locally asymptotically stable if $d_2 < \frac{c_2d_1}{c_1}$. When the growth of plant one is limited by the total phosphorus in the system and plant two is limited by light intensity, the boundary equilibrium is unstable and the coexistence equilibrium is locally asymptotically stable if $d_1 < \frac{c_1d_2}{c_2}$. When the growth of two plants are only limited by the total phosphorus in the system, the model does not have coexistence equilibrium and we give the conditions for the existence and stability of boundary equilibria.

When the competition for lights between two plants is equal with $b_{12} = b_{21}$, the numerical results indicate that light intensity and total phosphorus in the system play important roles in the growth and coexistence of two plants. When the total phosphorus in the system $T$ is fixed, moderate light intensity and sufficient total phosphorus in the system are beneficial for the coexistence of two plants. However, excessive light intensity can disrupt this kind of balance. Due to the limitation of the total phosphorus in the system, the loss rate of phosphorus determines the outcome of competition. Plants with a lower phosphorus loss rate will win the competition. This means that a higher light intensity can help plants with lower phosphorus loss rates win. In particular, when the parameters of both plants are equal, the two plants will coexist under sufficient light and nutrients. Moreover, as the light intensity increases, the quality of both plants will deteriorate. When the light intensity $K$ is fixed, two plants will go extinct due to the low total phosphorus in the system. With the enhancement of the total phosphorus in the system, the plant with the lower phosphorus loss rate can survive. This means that both plants cannot survive with lower total phosphorus in the system, while moderate total phosphorus in the system can help plants with lower phosphorus loss rates win. As the total phosphorus in the system further increases, the two plants will coexist. However, due to light limitation, the biomasses of the two plants will not continue increasing. At the same time, the increase in nutrients will improve the quality of two plants.

When the competition for light between two plants is not equal with $b_{12} \neq b_{21}$, the numerical simulation results also show that for a high level of total phosphorus in the system ($T = 0.03$), the dynamic properties of our model are the same as those in ^[1]. That is, when two plants have different competitive abilities the model exhibits bistability ($E_1$ and $E_2$ are stable). At this time, the competition results of the two plants depend on their initial values. For a moderate level of the total phosphorus in the system ($T = 0.001$), the bistable region of $E_1$ and $E_2$ will decrease. When two plants have a great competition ability, the competitive outcome of the two plants depends on their ability to compete for light. At a low level of total phosphorus in the system ($T = 0.0009$), the model becomes unstable, resulting in more types of bistability ($E_1$ and $E_*$ are stable, $E_2$ and $E_*$ are stable). At this time, the bistable region of $E_1$ and $E_2$ continues decreasing, and the coexistence area of the two plants becomes larger. Meanwhile, the survival of the two plants depends on their competition for light; the smaller the competition for light, the more favorable for the coexistence of the two plants. We conclude that when the total phosphorus in the system is at moderate or low level, competitive results depend on a plant's ability of competing for light rather than its initial state.

Our study still has several limitations. First, in order to facilitate the mathematical analysis of the model, we assumed that $b_{12} = b_{21}$. However, from a biological perspective, equal competition for light ($b_{12}$ equals $b_{21}$) is improbable, as no two plants possess identical parameters. Second, due to the complexity of the model, we didn't obtain the global asymptotic stability of the coexistence equilibrium. Third, for the mathematical analysis, we assumed that the phosphorus uptake rate function was linear, while for the nonlinear scenario, we used numerical simulations to investigate the dynamics of the model, which will be studied in future work.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

LX is funded by the National Natural Science Foundation of China 12171116. XR is funded by the National Natural Science Foundation of China 12101157. We thank Dr. Hao Wang for the valuable discussion.

Conflict of interest

The authors declare that there is no conflict of interest.

Appendix A.   Mathematical analysis of (2.4)

In this section, we verify the dissipativity of (2.4), then derive the equilibria and analyze their stabilities. For the convenience of subsequent analysis of the model, we let $b_{12} = b_{21} = 1$.

Theorem 1. The solutions with initial conditions in the set

remain there for all forward times.

Proof We consider a solution $X(t) = (x_{1}(t), x_{2}(t), Q_{1}(t), Q_{2}(t))$ of (2.2) with initial condition in $\varPhi(\Omega)$. Hence, $0 < x_{1}(0) < \min\{K, T/q_{1}\}, 0 < x_{2}(0) < \min\{K, T/q_{2}\}, Q_1(0) > q_1, Q_2(0) > q_2, Q_1(0)x_1(0)+Q_2(0)x_2(0) < T$. We assume that $\exists$ $t_{1} > 0$, such that $X(t)$ touches or crosses the boundary of $\bar\Omega$ (closure of $\varPhi(\Omega)$) for the first time, then $(x_{1}(t), x_{2}(t), Q_{1}(t), Q_{2}(t))$ $\in$ $\varPhi(\Omega)$ for $0\le t < t_{1}$. We discuss the following cases:

(i) $x_{1}(t_1) = 0$. $Q_{1}(t) > q_1$, $Q_{2}(t) > q_2$, $Q_{1}(t)x_1(t)+Q_{2}(t)x_2(t) < T$ for $0 \leq t \leq t_{1}$. Let $a_{1}$ = $\min$$\{ Q_{1}(t): t \in [0, t_{1}] \}$$> 0$. Then, for $0 \leq t \leq t_{1}$, we have

(ii) $x_{2}(t_{1}) = 0$. The proof is the same as that for (i).

(iii) $x_{1}(t_{1}) = \min\{K, T/q_{1}\}$. $0 < x_{2}(t) < \min\{k, T/q_{2}\}$, $Q_1(t)x_1(t)+Q_2(t)x_2(t) \leq T$ for $0 \leq t \leq t_{1}$ implies $Q_1(t)x_1(t) \le T$ for $0 \leq t \leq t_{1}$. Then, for $0 \leq t \leq t_{1}$, we have

if $\frac{dx_1}{dt} = r_{1}(1-\frac{x_{1}+x_{2}}{K})$, then $x_{1}(t)+x_{2}(t) < K$, $x_1(t) < K$. Moreover, $Q_1(t)x_1(t) \le T$, then $x_1(t) < T/q_1$. Therefore, $x_{1} < \min\{K, T/q_{1}\}$ for all $0 \leq t \leq t_{1}$, which is a contradiction.

(iv) $x_{2}(t_{1}) = \min\{K, T/q_{2}\}$. The proof is the same as that for (iii).

(v) $Q_{1}(t_{1}) = q_1$. $Q(t) = Q_1(t)-q_1$ then $Q(t_1) = 0$ and $Q(t) > 0$ for $0 \leq t < t_{1}$. For $0 \leq t \leq t_{1}$, we have

where $a$ is a constant. Thus, $Q(t) \ge Q_1(0)e^{-d_{1}-at} > 0$ for $0 \leq t \leq t_{1}$, which is a contradiction.

(vi) $Q_{2}(t_{1}) = q_2$. The proof is the same as that for (v).

(vii) $Q_1(t_1)x_1(t_1)+Q_2(t_1)x_2(t_1) = T$. $z(t) = T-Q_1(t)x_1(t)+Q_2(t)x_2(t)$ then $z(t_{1}) = 0$ and $z(t) > 0$ for $0 \leq t < t_{1}$. Then, for $0 \leq t \leq t_{1}$, we have

where $\tilde{v} > 0$ is a constant. Thus, $z(t)\ge z(0)e^{-\tilde{v}t} > 0$ for $0 \leq t \le t_{1}$, which is a contradiction. □

For simplicity, we assume that $g_i(x) = c_ix$, $i = 1, 2$. In this case, (2.4) becomes

We analyze the stability of the equilibrium of (A.1) in the following four cases.

Case 1. $\frac{ x_{1}+ x_{2}}{ K} > \frac{ q_{1}}{ Q_1}$, $\frac{ x_{1}+ x_{2}}{ K} > \frac{ q_{2}}{ Q_2}$

In Case 1, the boundary equilibria are $E_{1} = (0, K, \frac{c_{1}Td_2}{d_1(c_2K+d_2)}, \frac{c_2T}{c_2K+d_2})$ and $E_{2} = (K, 0, \frac{c_1T}{c_1K+d_1}, \frac{c_{2}Td_1}{d_{2}(c_1K+d_1)})$. The coexistence equilibrium is $E_{*} = (\tilde{x}_1, \tilde{x}_2, \tilde{Q}_1, \tilde{Q}_2)$, where

The Jacobian matrix at $(x_1, x_2, Q_1, Q_2)$ is

where

The boundary equilibria, $E_1$ and $E_2$, are unstable and their stability is complex. We will discuss their stability in the numerical simulation.

Next, we explore the globally asymptotic stability of $E_*$. The main approach involves the theory of asymptotic autonomous systems and Dulac's criterion ^[25]. We shall first summarize a few concepts and results of asymptotic autonomous system from Markus (2016) ^[12]. Consider the following differential equations

Equation (A.3) is asymptotically autonomous with limit equation (A.4) if

i.e., for $x$ in any compact subset of $R^{n}$. Assume that $f(t, x), g(x)$ are continuous functions and locally Lipschitz in $x$.

Lemma 1. (^[23]) The $\omega$-limit set $\Omega$ of a forward bounded solution $x$ to (A.3) is not empty, compact and connected. Moreover, $\Omega$ attracts $x$, i.e.,

Hence, $\Omega$ is invariant under (A.3). In particular, any point in $\Omega$ lies on an orbit of (A.4) in $\Omega$.

Lemma 2. (^[23]) Let $\hat{E}$ be a locally asymptotically stable equilibrium of (A.4) and let $\Omega$ be the $\omega$-limit set of forward bounded solution $x$ of (A.3). If $\Omega$ contains a point $y_0$, such that the solution $y$ of (A.4) with $y(0) = y_0$ converges to $\hat{E}$ when $t \to \infty$. Then, $\Omega = \{\hat{E}\}$, i.e., $x(t) \to \hat{E}$ when $t \to \infty$.

The existence and stability of coexistence equilibrium $E_{*} = (\tilde{x}_1, \tilde{x}_2, \tilde{Q}_1, \tilde{Q}_2)$ is analyzed below, where

Theorem 2. If $\tilde{Q_1} > q_1$, $\tilde{Q_2} > q_2$ and $\frac{c_2d_1\tilde{Q_1}^2-c_1d_2\tilde{Q_2}^2}{c_1c_2(\tilde{Q_2}-\tilde{Q_1})} < 0$, $E_*$ exists and it is globally asymptotically stable.

Proof The Jacobian matrix at $E_{*}$ is

The eigenvalues of $J(E_*)$ are

We divide (A.2) into two subsystems. Next, we consider the following submodel of (A.2)

Let

then

Therefore, (A.5) has no limit cycle by Dulac's criterion ^[25]. Both boundary equilibria of (A.5) are unstable. Therefore, the solution of (A.5) will tend to the globally asymptotically stable internal equilibrium, implying that the solution of (A.5) will approach the line $x_1 = K-x_2$.

Note that the first two equations of (A.2) are independent on $Q_1$ and $Q_2$. Moreover, $(\tilde{x}_1, \tilde{x}_2)$ is globally asymptotically stable. Next, we analyze the last equations of (A.2)

According to the above analysis in (A.5), the solution of the (A.5) will approach the line $x_1 = K-x_2$. Let $\lim_{t \to \infty} x_1(t) = m$, $\lim_{t \to \infty} x_2(t) = K-m = n$. Here, $m < \min\{K, T/q_1\}$ and $n < \min\{K, T/q_2\}$. Moreover, $Q_1m+Q_2n < T$. Then, the limiting system of (A.6) is

Let

then

Clearly, $nQ_2 < T$ and $mQ_1 < T$. It implies that $\frac{\partial(DP_1)}{\partial Q_1}+\frac{\partial(DP_2)}{\partial Q_2} < 0$. Therefore, it can be concluded that (A.7) has no limit cycle by applying Dulac's criterion ^[25] and the solution of (A.7) will tend to be globally asymptotically stable. This means that the solution of (A.6) will approach the line $Q_1 = \frac{c_1d_2}{c_2d_1} Q_2$. By the theory of asymptotic autonomous systems, the internal equilibrium $(\tilde{Q}_1, \tilde{Q}_2)$ of (A.6) is globally asymptotically stable. Hence, (A.2) has no nontrivial periodic solutions in $\varPhi(\Omega)$. The internal equilibrium $E_*$ of (A.2) is globally asymptotically stable.

Case 2. $\frac{ x_{1}+ x_{2}}{ K} > \frac{ q_{1}}{ Q_1}$, $\frac{ x_{1}+ x_{2}}{ K} < \frac{ q_{2}}{ Q_2}$

In Case 2, the boundary equilibrium is $E_{1}^{1} = \Big(0, \frac{T}{q_{2}}-\frac{d_{2}}{c_{2}}, \frac{c_1Kd_2q_{2}^{2}}{c_2q_2K(d_1+r_1)-r_1(c_2T-d_2q_2)}, q_{2}\Big)$ and the coexistence equilibrium is $E_{*}^{1} = (x_1^{*1}, x_2^{*1}, Q_1^{*1}, Q_2^{*1})$, where

The Jacobian matrix at $(x_1, x_2, Q_1, Q_2)$ is

Theorem 3. If $0 < \frac{T}{q_2}-\frac{d_2}{c_2} < K$ and $K < \frac{(c_2T-d_2q_2)(c_1d_1q_2+r_1q_1c_2)}{q_1q_2c_{2}^{2}(d_1+r_1)}$, $E_1^{1}$ exists and it is unstable.

Proof The Jacobian matrix at $E_{1}^{1}$ is

The eigenvalues of $J(E_1^{1})$ are

Since $\lambda_{1} > 0$, $E_{1}^{1}$ is unstable.

The existence and local stability of coexistence equilibrium $E_{*}^{1} = (x_1^{*1}, x_2^{*1}, Q_1^{*1}, Q_2^{*1})$ is analyzed below, where

Theorem 4. For the case of $d_2 > \frac{c_2d_1q_1}{c_1q_2}$, $0 < \frac{d_2q_2(d_1+c_1K)-c_2d_1T}{q_2(d_2c_1-d_1c_2)} < K$, $d_2c_1-d_1c_2 \neq 0$, $E_{*}^{1}$ exists. $E_{*}^{1}$ is locally asymptotically stable if $d_2 < \frac{c_2d_1}{c_1}$.

Proof The Jacobian matrix at $E_{*}^{1}$ is

The corresponding characteristic equation is

where

We verify that

If $d_2 < \frac{c_2d_1}{c_1}$, $A_4 > 0$. According to the Routh-Hurwitz criterion ^[24], $E_{*}^{1}$ is locally asymptotically stable.

Case 3. $\frac{ x_{1}+ x_{2}}{ K} < \frac{ q_{1}}{ Q_1}$, $\frac{ x_{1}+ x_{2}}{ K} > \frac{ q_{2}}{ Q_2}$

In case 3, the boundary equilibrium is $E_{2}^{2} = \Big(\frac{T}{q_1}-\frac{d_1}{c_1}, 0, q_1, \frac{c_2Kq_1^{2}d_1}{Kq_1c_1(d_2+r_2)-r_2(c_1T-d_1q_1)}\Big)$ and the coexistence equilibrium is $E_{*}^{2} = (x_1^{*2}, x_2^{*2}, Q_1^{*2}, Q_2^{*2})$, where

The Jacobian matrix at $(x_1, x_2, Q_1, Q_2)$ is

It is verified that $E_{0}^{2}$ and $E_1^{2}$ do not exist.

Theorem 5. If $0 < \frac{T}{q_1}-\frac{d_1}{c_1} < K$ and $K < \frac{(c_1T-d_1q_1)(c_2q_1d_1+c_1r_2q_2)}{c_1^{2}q_1q_2(d_2+r_2)}$, $E_{2}^{2}$ exists and it is unstable.

Proof The Jacobian matrix at $E_{2}^{2}$ is

The eigenvalues of $J(E_2^{2})$ are

Thus, $\lambda_{1} > 0$ and $E_{2}^{2}$ is unstable.

The existence and local stability of coexistence equilibrium $E_{*}^{2} = (x_1^{*2}, x_2^{*2}, Q_1^{*2}, Q_2^{*2})$ is analyzed below, where

Theorem 6. For the case of $d_1 > \frac{c_1d_2q_2}{c_2q_1}$, $0 < \frac{d_1q_1(d_2+c_2K)-c_1d_2T}{q_1(d_1c_2-d_2c_1)} < K$, $d_1c_2-d_2c_1 \neq 0$, $E_{*}^{2}$ exists. $E_{*}^{2}$ is locally asymptotically stable if $d_1 < \frac{c_1d_2}{c_2}$.

Proof The Jacobian matrix at $E_{*}^{2}$ is

The corresponding characteristic equation is

where

We verify that

If $d_1 < \frac{c_1d_2}{c_2}$, $B_4 > 0$. According to the Routh-Hurwitz criterion ^[24], if $d_1 < \frac{c_1d_2}{c_2}$, $E_{*}^{2}$ is locally asymptotically stable.

Case 4. $\frac{ x_{1}+ x_{2}}{ K} < \frac{ q_{1}}{ Q_1}$, $\frac{ x_{1}+ x_{2}}{ K} < \frac{ q_{2}}{ Q_2}$

In Case 4, (A.10) does not have coexistence equilibrium. The boundary equilibria are

The Jacobian matrix at $(x_1, x_2, Q_1, Q_2)$ is

In the following analysis, let $T_1 = \frac{d_1q_1}{c_1}$ and $T_2 = \frac{d_2q_2}{c_2}$.

Theorem 7. If $T \geq max\{T_1, T_2\}$, $E_{0}^{3}$ exists. If $T > T_1$ or $T > T_2$, $E_0^{3}$ is unstable in $\varPhi(\Omega)$.

Proof The Jacobian matrix at $E_{0}^{3}$ is

The eigenvalues of $J(E_0^{3})$ are

If $T > T_1$ or $T > T_2$, $\lambda_{1} > 0$ or $\lambda_{2} > 0$, $E_0^{3}$ is unstable. If $T = T_1$ or $T = T_2$, the stability cannot be determined.

Theorem 8. If $T_1 \leq T_2 < T < q_2K+T_2$ and $K > \frac{(c_2T-q_2d_2)(c_1d_2q_2+c_2r_1q_1)}{q_1q_2c_2^{2}(d_1+r_1)}$, $E_{1}^{3}$ exists. If $T_2 > T_1$, $E_{1}^{3}$ is unstable in $\varPhi(\Omega)$.

Proof The Jacobian matrix at $E_{1}^{3}$ is

The eigenvalues of $J(E_1^{3})$ are

If $T_2 > T_1$ and $\lambda_{1} > 0$, $E_{1}^{3}$ is unstable. If $T_2 = T_1$, the stability of $E_{1}^{3}$ cannot be determined.

Theorem 9. If $T_2 \leq T_1 < T < q_1K+T$ and $K > \frac{(c_1T-q_1d_1)(c_2d_1q_1+c_1r_2q_2)}{q_1q_2c_1^{2}(d_2+r_2)}$, $E_{2}^{3}$ exists. If $T_1 > T_2$, $E_{2}^{3}$ is unstable in $\varPhi(\Omega)$.

Proof The Jacobian matrix at $E_{2}^{3}$ is

The eigenvalues of $J(E_2^{3})$ are

If $T_1 > T_2$, $\lambda_{1} > 0$, $E_{2}^{3}$ is unstable. If $T_1 = T_2$, the stability of $E_{2}^{3}$ cannot be determined.