
We are interested in modelling the interaction of bacteria and certain nutrient concentration within a porous medium admitting viscous flow. The governing equations in primal-mixed form consist of an advection-reaction-diffusion system representing the bacteria-chemical mass exchange, coupled to the Brinkman problem written in terms of fluid vorticity, velocity and pressure, and describing the flow patterns driven by an external source depending on the local distribution of the chemical species. A priori stability bounds are derived for the uncoupled problems, and the solvability of the full system is analysed using a fixed-point approach. We introduce a primal-mixed finite element method to numerically solve the model equations, employing a primal scheme with piecewise linear approximation of the reaction-diffusion unknowns, while the discrete flow problem uses a mixed approach based on Raviart-Thomas elements for velocity, Nédélec elements for vorticity, and piecewise constant pressure approximations. In particular, this choice produces exactly divergence-free velocity approximations. We establish existence of discrete solutions and show their convergence to the weak solution of the continuous coupled problem. Finally, we report several numerical experiments illustrating the behaviour of the proposed scheme.
Citation: Verónica Anaya, Mostafa Bendahmane, David Mora, Ricardo Ruiz Baier. 2018: On a vorticity-based formulation for reaction-diffusion-Brinkman systems, Networks and Heterogeneous Media, 13(1): 69-94. doi: 10.3934/nhm.2018004
[1] | Ilse Domínguez-Alemán, Itzel Domínguez-Alemán, Juan Carlos Hernández-Gómez, Francisco J. Ariza-Hernández . A predator-prey fractional model with disease in the prey species. Mathematical Biosciences and Engineering, 2024, 21(3): 3713-3741. doi: 10.3934/mbe.2024164 |
[2] | Saheb Pal, Nikhil Pal, Sudip Samanta, Joydev Chattopadhyay . Fear effect in prey and hunting cooperation among predators in a Leslie-Gower model. Mathematical Biosciences and Engineering, 2019, 16(5): 5146-5179. doi: 10.3934/mbe.2019258 |
[3] | Lazarus Kalvein Beay, Agus Suryanto, Isnani Darti, Trisilowati . Hopf bifurcation and stability analysis of the Rosenzweig-MacArthur predator-prey model with stage-structure in prey. Mathematical Biosciences and Engineering, 2020, 17(4): 4080-4097. doi: 10.3934/mbe.2020226 |
[4] | Xin-You Meng, Yu-Qian Wu . Bifurcation analysis in a singular Beddington-DeAngelis predator-prey model with two delays and nonlinear predator harvesting. Mathematical Biosciences and Engineering, 2019, 16(4): 2668-2696. doi: 10.3934/mbe.2019133 |
[5] | Xiaoyuan Chang, Junjie Wei . Stability and Hopf bifurcation in a diffusivepredator-prey system incorporating a prey refuge. Mathematical Biosciences and Engineering, 2013, 10(4): 979-996. doi: 10.3934/mbe.2013.10.979 |
[6] | Ceyu Lei, Xiaoling Han, Weiming Wang . Bifurcation analysis and chaos control of a discrete-time prey-predator model with fear factor. Mathematical Biosciences and Engineering, 2022, 19(7): 6659-6679. doi: 10.3934/mbe.2022313 |
[7] | Shuo Yao, Jingen Yang, Sanling Yuan . Bifurcation analysis in a modified Leslie-Gower predator-prey model with fear effect and multiple delays. Mathematical Biosciences and Engineering, 2024, 21(4): 5658-5685. doi: 10.3934/mbe.2024249 |
[8] | Gunog Seo, Mark Kot . The dynamics of a simple Laissez-Faire model with two predators. Mathematical Biosciences and Engineering, 2009, 6(1): 145-172. doi: 10.3934/mbe.2009.6.145 |
[9] | Christian Cortés García, Jasmidt Vera Cuenca . Impact of alternative food on predator diet in a Leslie-Gower model with prey refuge and Holling Ⅱ functional response. Mathematical Biosciences and Engineering, 2023, 20(8): 13681-13703. doi: 10.3934/mbe.2023610 |
[10] | Kawkab Al Amri, Qamar J. A Khan, David Greenhalgh . Combined impact of fear and Allee effect in predator-prey interaction models on their growth. Mathematical Biosciences and Engineering, 2024, 21(10): 7211-7252. doi: 10.3934/mbe.2024319 |
We are interested in modelling the interaction of bacteria and certain nutrient concentration within a porous medium admitting viscous flow. The governing equations in primal-mixed form consist of an advection-reaction-diffusion system representing the bacteria-chemical mass exchange, coupled to the Brinkman problem written in terms of fluid vorticity, velocity and pressure, and describing the flow patterns driven by an external source depending on the local distribution of the chemical species. A priori stability bounds are derived for the uncoupled problems, and the solvability of the full system is analysed using a fixed-point approach. We introduce a primal-mixed finite element method to numerically solve the model equations, employing a primal scheme with piecewise linear approximation of the reaction-diffusion unknowns, while the discrete flow problem uses a mixed approach based on Raviart-Thomas elements for velocity, Nédélec elements for vorticity, and piecewise constant pressure approximations. In particular, this choice produces exactly divergence-free velocity approximations. We establish existence of discrete solutions and show their convergence to the weak solution of the continuous coupled problem. Finally, we report several numerical experiments illustrating the behaviour of the proposed scheme.
Since Lotka and Volterra [1,2] proposed the classical Lotka-Volterra predator-prey model, it has received much attention in studying the relationships between the predator and prey [3,4,5,6,7]. In the meaning of biology, the population cannot grow unrestrictedly. However, the prey grows infinitely in the absence of the predator for the original Lotka-Volterra predator-prey model. To correct this unrealistic defect, the modified Lotka-Volterra predator-prey model with the finite environment carrying capacity for prey is considered and studied (refer to May [8] for a detailed model construction), which is in the form
$ dxdτ=Rx(1−xK)−hxy,dydτ=hxy−dy, $
|
(1.1) |
where $ R $ and $ K $ are the per capita growth rate and the carrying capacity of the prey $ x $, $ h $ and $ d $ are the hunting ability and the death rate of the predator $ y $, respectively. As part of the subsystem of the proposed model in this article, more specific analysis and results for model (1.1) will be given in Section 3.
Recently, Mobilia et al. [9] expanded Lotka-Volterra predator-prey model and studied the effects of the spatial constraints and stochastic noise on the properties of the predator-prey systems. Considering the spatial constraints and stochastic noise, Mobilia et al. [9] formulated a predator-prey model similar to model (1.1), and a stochastic lattice Lotka-Volterra model (SLLVM). They obtained a lot of important and interesting results differing from those of the classical (unrestricted) deterministic Lotka-Volterra predator-prey model by a suitable mean-field approach, field-theoretic arguments and Monte Carlo simulations. In the topical review [10], from the historical overview of population dynamics to the stochastic lattice Lotka-Volterra predator-prey models, Dobramysl et al. focused on spatially extended population dynamics models and the role of fluctuations and correlations in biological systems induced by the demographic noise. Also, they demonstrated the cyclic dominance of three-species populations and multiple species competition networks. The models in the above two articles are the improvements of the Lotka-Volterra predator-prey by considering the finite carrying capacity and incorporating the stochastic noise, spatial constraints and spatial extension.
The Lotka-Volterra predator-prey models have been explored and expanded widely in the aspects of the diffusion, time delay, harvesting, switching and herd behavior [11,12,13,14,15,16,17,18,19,20,21], etc. However, considering the effect of the population interactions for different genotypes on the ecosystem is a research hotspot that has only been developed in recent years.
Taking into account the fact that the prey is genetically distinguishable and assuming that the prey is divided into two subpopulations with two different genotypes $ x $ and $ y $, Venturino [22] proposed an ecogenetic model
$ dxdτ=(Rp−ax)(x+y)−hxz,dydτ=(Rq−by)(x+y)−gyz,dzdτ=z[e(hx+gy)−m], $
|
(1.2) |
where $ z $ represents the predator, $ R $ is the reproduction rate of the prey, $ p $ and $ q = 1-p $ are the fractions of the newborns for the two genotypes $ x $ and $ y $, $ a $ and $ b $ are the population pressures of the two genotypes of the prey, $ h $ and $ g $ are the predators' different hunting ability to prey on the two genotypes of food, respectively; $ m $ is the mortality rate of the predator, $ 0 < e < 1 $ is the conversion rates of the prey into the predator's newborns. The author showed that if the predator invasion number is greater than $ 1 $ then the predator is permanent and the model is persistent. Furthermore, model (1.2) could be persistent oscillatory under some conditions. Reference [23] also proposed an eco-genetic model with distinguishable genotypes happened in the prey. Supposed that the genetically distinguishable species are the predator with two different genotypes $ y $ and $ z $, Viberti and Venturino [24] introduced a model
$ dxdτ=R(1−xK)x−hxy−gxz,dydτ=pe(hy+gz)x−my,dzdτ=qe(hy+gz)x−nz, $
|
(1.3) |
where $ R $ and $ K $ are the reproduction rate and the environment's carrying capacity of the prey $ x $, respectively, and all other parameters have the similar meaning to those in model (1.2). Then authors investigated the extinction of the predator, the permanence of the prey and the coexistence of the two species. However, model (1.3) does not allow persistent oscillation, which is the main difference between models (1.2) and (1.3). Similar to model (1.3), Viberti and Venturino [25] proposed a predator-prey model with Holling Ⅱ response function in describing the dynamics of two different genotypes predator.
Motivated by [22,24,25], we propose a predator-prey model in which the genetic differentiation both happened in the predator and prey. The model is presented in the form as follows
$ dxdτ=(Rp1−ax)(x+y)−mxu−nxv,dydτ=(Rq1−by)(x+y)−myu−nyv,dudτ=p2e(mxu+myu+nxv+nyv)−ku,dvdτ=q2e(mxu+myu+nxv+nyv)−lv, $
|
(1.4) |
where the two genotypes population $ x $ and $ y $ of the prey have the same reproduction rate $ R $, the two genotypes population $ u $ and $ v $ of the predator have the same conversion rates $ 0 < e < 1 $; $ p_{1} $ and $ q_{1} = 1-p_1 $ are the fractions of the newborns of the prey $ x $ and $ y $ respectively, $ p_{2} $ and $ q_{2} = 1-p_2 $ are the fractions of the newborns of the predator $ u $ and $ v $ respectively, $ a $ and $ b $ are the different pressures felt by the two types of prey population, $ k $ and $ l $ are the mortality rates of $ u $ and $ v $ respectively; we also assume that the predator $ u $ has the same hunting ability to prey $ x $ and prey $ y $, which is denoted by $ m $; Similarly, the hunting ability of predator $ v $ to prey $ x $ and prey $ y $ is denoted by $ n $. All coefficients in model (1.4) are positive. Assuming that both the two types of predators have different hunting abilities to the two types of the prey, Castellino et al. [26] gave a more general predator-prey model with both species genetically distinguishable. However, they only roughly analyzed the stability of the boundary equilibria and showed some results by numerical simulations. In this paper, we shall give a detailed analysis of the model (1.4).
In order to simplify calculations, we first nondimensionalize model (1.4) by letting $ t = e\tau $, $ x(\tau) = \alpha X(t) $, $ y(\tau) = \alpha Y(t) $, $ u(\tau) = \gamma U(t) $, $ v(\tau) = \gamma V(t) $. Then it becomes
$ dXdt=1αe(Rp1−aαX)(αX+αY)−1αemαγXU−1αenαγXV,dYdt=1αe(Rq1−bαY)(αX+αY)−1αemαγYU−1αenαγYV,dUdt=1γep2e(mαγXU+mαγYU+nαγXV+nαγYV)−1γekγU,dVdt=1γeq2e(mαγXU+mαγYU+nαγXV+nαγYV)−1γelγV. $
|
Let $ \alpha = \frac{Rp_{1}}{a} $, $ \gamma = \frac{e}{m} $, and define new parameters $ r = \frac{Rp_{1}}{e} $, $ c = \frac{n}{m} $, $ w = \frac{q_{1}}{p_{1}} $, $ s = \frac{b}{a} $, $ p = \frac{p_{2}mRp_{1}}{a} $, $ q = \frac{q_{2}mRp_{1}}{a} $, $ g = \frac{k}{e} $, $ d = \frac{l}{e} $. Then the above model can be further simplified into
$ dXdt=r(1−X)(X+Y)−XU−cXV,dYdt=r(w−sY)(X+Y)−YU−cYV,dUdt=p[(X+Y)U+c(X+Y)V]−gU,dVdt=q[(X+Y)U+c(X+Y)V]−dV. $
|
(1.5) |
It is clear that all of models (1.2), (1.3), (1.4), and (1.5) have more than one predator or one prey, and are generalizations of Lotka-Volterra predator-prey model apparently.
Azzali et al. [27] studied a competitive model with genetically distinguishable species. In fact, model (1.4) is similar to the ecoepidemiology model, i.e., the population model with diseases spreading in the species. Earlier researches on the ecoepidemiology models can be found from [28,29], and influences of diseases on the predator-prey models and competitive models can refer to the articles [30,31,32,33,34,35,36,37,38,39,40].
The rest of this paper is organised as follows. In Sections 2 and 3, we shall analyze model (1.5) without predators and model (1.5) with one genotype in both the predator and prey, respectively. The full model of (1.5) will be discussed in Section 4. Some examples to illustrate our main results will be given in Section 5. In Section 6, we shall give our conclusions and discussions.
Without the predators, model (1.5) becomes
$ dXdt=r(1−X)(X+Y),dYdt=r(w−sY)(X+Y). $
|
(2.1) |
It is easy to know that $ (0, 0) $ and $ (1, \frac{w}{s}) $ are two equilibria of the model. And we can obtain the following two lemmas.
Lemma 2.1. For model (2.1), the boundary equilibrium $ (0, 0) $ is unstable, and the positive equilibrium $ \left(1, \frac{w}{s}\right) $ is locally stable.
Proof. The Jacobian matrix of (2.1) is
$ J_1 = \left[ r−2rX−rYr−rXrw−rsYrw−rsX−2rsY \right]. $
|
Then the characteristic polynomial of $ J_1 $ at $ (0, 0) $ is $ \lambda^2-(r+wr)\lambda $, which has two eigenvalues $ \lambda_{1} = 0 $, $ \lambda_{2} = r(1+w) > 0 $. So, the origin $ (0, 0) $ is unstable.
Similarly, at $ (1, \frac{w}{s}) $ the characteristic polynomial is $ \left(\lambda+r+\frac{wr}{s}\right)(\lambda+wr+sr) $, resulting in two negative eigenvalues $ \lambda_{1} = -r\left(1+\frac{w}{s}\right) $, $ \lambda_{2} = -r(s+w) $. Therefore, equilibrium $ \left(1, \frac{w}{s}\right) $ is locally stable.
Lemma 2.2. Model (2.1) does not admit any periodic solution.
Proof. Let
$ B_1(X,Y) = \frac{1}{r(X+Y)},\ W_{1}(X,Y) = r(1-X)(X+Y),\ W_{2}(X,Y) = r(w-sY)(X+Y), $ |
it is easy to know that $ B_1(X, Y) $ is a continuously differentiable function in the first quadrant. We get that
$ \frac {\partial [W_1(X,Y)B_1(X,Y)]}{\partial{X}}+\frac {\partial[W_2(X,Y)B_1(X,Y)]}{\partial{Y}} = -1-s \lt 0. $ |
Then the Bendixson-Dulac criterion ([41], Proposition 1.195) implies our result here.
Then, Lemmas 2.1 and 2.2 yield the following global behavior of model (2.1).
Theorem 2.1. Positive equilibrium $ \left(1, \frac{w}{s}\right) $ of model (2.1) is globally asymptotically stable.
Remark 2.1. In the biological sense, Theorem 2.1 means that two genotypes of the prey will always coexist without the predators. This is consistent with the Hardy-Weinberg law [42], which shows that the gene and genotype frequencies in a species will remain constant from the first daughter generation onwards.
When both the predator and prey have only one genotype, model (1.5) becomes the nondimensionalization of the Lotka-Volterra predator-prey model (1.1), which is in the following form
$ dXdt=r(1−X)X−XU,dUdt=pXU−gU. $
|
(3.1) |
It is easy to obtain that the model always has two boundary equilibria $ (0, 0) $ and $ (1, 0) $, and has a positive equilibrium $ \left(\frac{g}{p}, \frac{r(p-g)}{p}\right) $ if $ p > g $. The positive equilibrium $ \left(\frac{g}{p}, \frac{r(p-g)}{p}\right) $ approaches to the boundary equilibrium $ (1, 0) $ as $ p\to g $.
Next, by using Jacobin matrix and center manifold theory, we have Lemma 3.1 and 3.2 as shown in the following.
Lemma 3.1. For model (3.1), boundary equilibrium $ (0, 0) $ is always unstable, and the other boundary equilibrium $ (1, 0) $
(a) is stable when $ p < g $;
(b) is unstable when $ p > g $;
(c) has a local center manifold
$ W_{c}^{loc} = \left\{(x_{2},u_{2})\in{R^2}:u_{2}\in{K},x_{2} = -\frac{p}{r}u_{2}^2+\left(\frac{3p^2}{r^2}+\frac{p}{r}\right)u_{2}^3+O(u_{2}^4)\right\} $ |
when $ p = g $, where $ K $ is a small neighborhood of the origin, and
$ [x2u2]=P−1[X−1U], P=[110−r]. $
|
Proof. The Jacobian matrix of (3.1) is
$ J_2 = \left[ r−2rX−U−XpUpX−g \right]. $
|
The characteristic polynomial of matrix $ J_2 $ at $ (0, 0) $ is $ (\lambda-r)(\lambda+g) $, which has two eigenvalues $ \lambda_{1} = r $, $ \lambda_{2} = -g $. Therefore, equilibrium $ (0, 0) $ is a saddle point.
The characteristic polynomial of matrix $ J_2 $ at $ (1, 0) $ is $ (\lambda+r)(\lambda+g-p) $, which has two eigenvalues $ \lambda_{1} = -r < 0, $ $ \lambda_{2} = p-g $. We have the following three cases.
(a) If $ p < g $, then $ \lambda_{2} < 0 $, implying equilibrium $ (1, 0) $ is locally stable.
(b) If $ p > g $, then $ \lambda_{2} > 0 $. Thus, equilibrium $ (1, 0) $ is a saddle point. Hence, unstable accordingly.
(c) If $ p = g $, then $ \lambda_{2} = 0 $. By the center manifold theory [43,44], there is a local center manifold at equilibrium $ (1, 0) $, which is calculated as follows. Let $ x_{1} = X-1 $, $ u_{1} = U $. Model (3.1) becomes
$ [x′1u′1]=A[x1u1]+[−rx21−x1u1px1u1], $
|
(3.2) |
where
$ A = \left[ −r−100 \right]. $
|
Equilibrium $ (1, 0) $ of model (1.5) becomes equilibrium $ (0, 0) $ of (3.2). The eigenvalues of $ A $ are the same as those of matrix $ J_2 $ at$ (1, 0) $, which are $ \lambda_{1} = -r, $ $ \lambda_{2} = 0 $. By the center manifold theory, equilibrium $ (0, 0) $ has a one-dimensional center subspace and a one-dimensional stable subspace. The eigenvectors of $ A $ with respect to $ \lambda_{1} $ and $ \lambda_{2} $ are
$ e_{1} = \left[ 10 \right],\ e_{2} = \left[ 1−r \right], $
|
respectively. Let $ P = (e_{1}, e_{2}) $, then $ P^{-1}AP = B $, where $ B = {diag}(\lambda_{1}, \lambda_{2}) $. Therefore, Let
$ \left[ x1u1 \right] = P\left[ x2u2 \right] = \left[ x2+u2−ru2 \right], $
|
then
$ [x′2u′2]=P−1AP[x2u2]+P−1[−rx2(x2+u2)−rpu2(x2+u2)]=B[x2u2]+[−(x2+u2)(rx2+pu2)p(x2+u2)u2]. $
|
The local center manifold at $ (0, 0) $ is described by an approximation function $ \psi:span{\{e_{2}\}}\rightarrow{span{\{e_{1}\}}} $. Since $ \psi(0) = D_{\psi}(0) = 0 $, we set the Taylor expansion near $ (0, 0) $ as
$ x_{2} = \psi(u_{2}) = \beta_{2}u_{2}^2+\beta_{3}u_{3}^3+\ldots. $ |
By the center manifold approximation theorem, we have
$ -r\psi(u_{2})-(\psi(u_{2})+u_{2})(r\psi(x_{2})+pu_{2}) = (2\beta_{2}u_{2}+3\beta_{3}u_{3}^2+\ldots)p(\psi(u_{2})+u_{2})u_{2}. $ |
Simplifying the above equation and comparing the coefficients of the two sides, we get
$ \beta_{2} = -\frac{p}{r},\ \beta_{3} = \frac{3p^2}{r^2}+\frac{p}{r}. $ |
Thus, the center manifold at $ (0, 0) $ is
$ W_{c}^{loc} = \left\{(x_{2},u_{2})\in{R^2}:u_{2}\in{K},x_{2} = -\frac{p}{r}u_{2}^2+\left(\frac{3p^2}{r^2}+\frac{p}{r}\right)u_{2}^3+O(u_{2}^4)\right\}, $ |
where $ K $ is a small neighborhood of the origin, and
$ [x2u2]=P−1[X−1U]. $
|
Lemma 3.2. If $ p > g $, then model (3.1) has a stable positive equilibrium $ \left(\frac{g}{p}, \frac{r(p-g)}{p}\right) $.
Proof. From the above discussion, the positive equilibrium $ \left(\frac{g}{p}, \frac{r(p-g)}{p}\right) $ exists if and only if $ p > g $. For model (3.1), the characteristic polynomial of Jacobian matrix at $ \left(\frac{g}{p}, \frac{r(p-g)}{p}\right) $ is $ f(\lambda) = \lambda^2+\theta_1 \lambda+\sigma_1 $, where $ \theta_1 = r\frac{g}{p} $, $ \sigma_1 = rg\frac{p-g}{p} $. It is clear that $ \sigma_1 > 0 $. Let $ \Delta_1 = \theta_1^{2}-4\sigma_1 $, then $ \Delta_1 < \theta_1^{2} $. The roots of $ f(\lambda) $ are $ \lambda_{1, 2} = \frac{-\theta_1\pm\sqrt{\Delta_1}}{2} $. Then, we only need to discuss the real parts of $ \lambda_{1, 2} $.
1. If $ \Delta_1\geq 0 $, then both $ \lambda_{1, 2} $ are negative, and positive equilibrium $ \left(\frac{g}{p}, \frac{r(p-g)}{p}\right) $ is stable accordingly.
2. If $ \Delta_1 < 0 $, then $ \lambda_{1, 2} $ are two conjugate complex numbers with negative real part, and positive equilibrium $ \left(\frac{g}{p}, \frac{r(p-g)}{p}\right) $ is stable accordingly.
Above discussion indicates the local behavior of equilibria in model (3.1), next we illustrate the global dynamical behaviors of these equilibria. First, we need to discuss the existence of periodic solution in model (3.1).
Lemma 3.3. For model (3.1), there is no periodic solution.
Proof. Let
$ B_2(X,U) = \frac{1}{XU},\ H_{1}(X,U) = r(1-X)X-XU,\ H_{2}(X,U) = pXU-gU. $ |
It is easy to know that $ B_2(X, U) $ is a continuous differentiable function in the first quadrant. Then, in the first quadrant
$ \frac{\partial[H_{1}(X,U)B_2(X,U)]}{\partial{X}}+\frac{\partial[H_{2}(X,U)B_2(X,U)]}{\partial{X}} = -\frac{r}{U} \lt 0. $ |
By the Bendixson-Dulac criterion, model (3.1) has no periodic solution.
From Lemmas 3.1, 3.2 and 3.3, we know the global behavior of model (3.1).
Theorem 3.1. For model (3.1), boundary equilibrium $ (0, 0) $ is always unstable, and
(a) if $ p > g $, then boundary equilibrium $ (1, 0) $ is unstable and positive equilibrium $ \left(\frac{g}{p}, \frac{r(p-g)}{p}\right) $ is globally asymptotically stable;
(b) if $ p < g $, then boundary equilibrium $ (1, 0) $ is globally asymptotically stable and there is no positive equilibrium.
From Theorem 3.1, we know that if the death rate of predator is bigger than its conversion rate, then the predator will be extinct while the prey will be permanent; if the death rate of predators is smaller than the conversion rate, then the predator and the prey will coexist permanently.
It is easy to find that the origin $ E_{0}(0, 0, 0, 0) $ and $ E_{1}(1, \frac{w}{s}, 0, 0) $ are two boundary equilibria of model (1.5). If
$ r(1−X)(X+Y)−X(U+cV)=0,r(w−sY)(X+Y)−Y(U+cV)=0,p(X+Y)(U+cV)−gU=0,q(X+Y)(U+cV)−dV=0, $
|
(4.1) |
has a positive solution, then it is the positive equilibrium of model (1.5). Let $ (X_{*}, Y_{*}, U_{*}, V_{*}) $ is a solution of Eq (4.1), then
$ 1−X∗X∗=w−sY∗Y∗, X∗+Y∗=h, U∗=ph2rg(1−X∗)X∗, V∗=qh2rd(1−X∗)X∗, $
|
(4.2) |
where $ h = \frac{gd}{pd+cqg} $.
If $ 0 < X_{*} < 1 $, $ 0 < Y_{*} < \frac{w}{s} $, $ U_{*} > 0 $, $ V_{*} > 0 $, then $ (X_{*}, Y_{*}, U_{*}, V_{*}) $ is a positive solution of (4.1), which means it is the positive equilibrium of model (1.5).
From the first two equations of (4.2), we get
$ (1−s)Y2∗+[w+1−(1−s)h]Y∗−wh=0. $
|
(4.3) |
Let $ \Delta_{2} = [w+1-(1-s)h]^2+4wh(1-s), $ then $ \Delta_{2} = [w-1+(1-s)h]^2+4w > 0 $. Let
$ f(Y)=(1−s)Y2+[w+1−(1−s)h]Y−wh, $
|
(4.4) |
then $ f\left(\frac{w}{s}\right) > 0 $ if $ h < 1+\frac{w}{s} $, and $ f\left(\frac{w}{s}\right)\leq 0 $ if $ h\geq 1+\frac{w}{s} $.
Next, we shall prove that $ h < 1+\frac{w}{s} $ is the necessary and sufficient condition for the existence of positive equilibrium of model (1.5). First, we shall prove that if $ h < 1+\frac{w}{s} $ then $ 0 < X_{\ast} < 1 $, $ 0 < Y_{\ast} < \frac{w}{s} $, $ U_{*} > 0 $, $ V_{*} > 0 $, and model (1.5) has a positive equilibrium $ E_{*}(X_{*}, Y_{*}, U_{*}, V_{*}) $ accordingly. Three cases are considered.
Case 1. When $ s = 1 $, we get $ Y_{*} = \frac{wh}{w+1} $ from (4.3). So, by (4.2),
$ X_{*} = \frac{h}{w+1},\ U_{*} = \frac{phr}{g}(w+1-h),\ V_{*} = \frac{qhr}{d}(w+1-h). $ |
If $ h < 1+w $, then $ 0 < X_{*} < 1 $, $ 0 < Y_{*} = \frac{wh}{w+1} < \frac{w}{1+w}(1+w) = w $, $ U_{*} > 0 $, $ V_{*} > 0 $.
Case 2. When $ s < 1 $, i.e., $ 1-s > 0 $, it is easy to find that the parabola $ f(Y) $ is opening up and the curve intersects vertical axis at point $ (0, -wh) $. Thus, (4.3) has one positive root
$ Y_{*} = \frac{h}{2}+\frac{-(w+1)+\sqrt{\Delta_{2}}}{2(1-s)}. $ |
By (4.2) and the above equation, we get
$ X_{*} = \frac{h}{2}+\frac{w+1-\sqrt{\Delta_{2}}}{2(1-s)},\ U_{*} = \frac{ph^2r}{g}\left[\frac{2(1-s)}{w+1+(1-s)h-\sqrt{\Delta_{2}}}-1\right], $ |
$ V_{*} = \frac{qh^2r}{d}\left[\frac{2(1-s)}{w+1+(1-s)h-\sqrt{\Delta_{2}}}-1\right]. $ |
If $ h < 1+\frac{w}{s} $, then $ f\left(\frac{w}{s}\right) > 0 $. From the image of parabola $ f(Y) $, we know that $ f(Y_*) = 0 $ and $ f(Y) < 0 $ as $ Y\in(0, Y_{*}) $. So $ 0 < Y_{*} < \frac{w}{s} $. By the first equation of (4.2), we have $ 0 < X_{*} < 1 $. Accordingly, $ U_{*} > 0 $ and $ V_{*} > 0 $.
Case 3. When $ s > 1 $, i.e., $ 1-s < 0 $, the parabola $ f(Y) $ is opening down and the curve intersects vertical axis at point $ (0, -wh) $. Since $ 1-s < 0 $ and $ w+1-(1-s)h > 0 $, we know the center shaft of the parabola is positive. So, $ f(Y) $ has two positive roots
$ Y_{*{1,2}} = \frac{h}{2}+\frac{-(w+1)\pm\sqrt{\Delta_{2}}}{2(1-s)}, $ |
which means
$ X_{*{1,2}} = \frac{h}{2}+\frac{(w+1)\mp\sqrt{\Delta_{2}}}{2(1-s)},\ U_{*{1,2}} = \frac{ph^2r}{g}\left[\frac{2(1-s)}{w+1+(1-s)h\mp\sqrt{\Delta_{2}}}-1\right], $ |
$ V_{*{1,2}} = \frac{qh^2r}{d}\left[\frac{2(1-s)}{w+1+(1-s)h\mp\sqrt{\Delta_{2}}}-1\right] $ |
from Eq (4.2). By $ f\left(\frac{w}{s}\right) > 0 $ as $ h < 1+\frac{w}{s} $, $ f(Y) > 0 $ as $ Y\in (Y_{*1}, Y_{*2}) $ and $ f(Y) < 0 $ as $ Y\notin [Y_{*1}, Y_{*2}] $, we get that $ \frac{w}{s}\in(Y_{*{1}}, Y_{*{2}}) $, which means $ 0 < Y_{*{1}} < \frac{w}{s} < Y_{*{2}} $. The same reason as that of Case 2 gives $ 0 < X_{*{1}} < 1 $. Accordingly, $ U_{*{1}} > 0, $ $ V_{*{1}} > 0 $.
By the same analysis as the above three cases, we have that if $ h\geq 1+\frac{w}{s} $ then $ Y_{*}\ge \frac{w}{s} $ for Case 1 and 2, and both $ Y_{*{2}} > Y_{*{1}} > \frac{w}{s} $ for Case 3. From the second equation of (4.1), we know that Eqs (4.1) have no positive solutions. Therefore, model (1.5) has no positive equilibria.
From the above discussion, we get the following results.
Theorem 4.1. Model (1.5) always has two boundary equilibria $ E_{0}(0, 0, 0, 0) $ and $ E_{1}(1, \frac{w}{s}, 0, 0) $. Furthermore, $ h < 1+\frac{w}{s} $ is the necessary and sufficient condition for the existence of an unique positive equilibrium $ E_{*}(X_{*}, Y_{*}, U_{*}, V_{*}) $ of model (1.5), where $ 0 < X_{*} < 1 $, $ 0 < Y_{*} < \frac{w}{s} $, $ U_{*} > 0 $ and $ V_{*} > 0 $.
The Jacobian matrix of model (1.5) at $ E_{0} $ is
$ J(E_0) = \left[ rr00rwrw0000−g0000−d \right]. $
|
The characteristic polynomial of $ J(E_0) $ is $ \lambda(\lambda-r-wr)(\lambda+g)(\lambda+d) $, which has four eigenvalues $ \lambda_{1} = 0 $, $ \lambda_{2} = r+rw $, $ \lambda_{3} = -g $ and $ \lambda_{4} = -d $. Since $ \lambda_{2} > 0 $, $ E_{0} $ is unstable.
Theorem 4.2. The boundary equilibrium $ E_{0} $ of (1.5) is unstable.
The Jacobian matrix of model (1.5) at $ E_{1} $ is
$ J(E_1) = \left[ −r−rws0−1−c0−rw−rs−ws−cws00p(1+ws)−gcp(1+ws)00q(1+ws)cq(1+ws)−d \right]. $
|
The characteristic polynomial of $ J(E_1) $ is
$ \left(\lambda+r+r\frac{w}{s}\right)(\lambda+rw+rs)\left\{\lambda^2+\left[g+d-(p+cq)\left(1+\frac{w}{s}\right)\right]\lambda +gd-(pd+cqg)\left(1+\frac{w}{s}\right)\right\}. $ |
$ J(E_1) $ has two negative eigenvalues $ \lambda_{1} = -r-r\frac{w}{s} $ and $ \lambda_{2} = -rw-rs $. The other two eigenvalues of $ J(E_1) $ are the roots of polynomial $ f_1(\lambda) = \lambda^2-\theta\lambda+\sigma $, where
$ \theta = (p+cq)\left(1+\frac{w}{s}\right)-(g+d),\ \sigma = gd-(dp+cqg)\left(1+\frac{w}{s}\right). $ |
The two roots of $ f_1(\lambda) $ are
$ λ3,4=θ±√θ2−4σ2. $
|
(4.5) |
Let $ \Delta_{3} = \theta^2-4\sigma, $ then
$ \Delta_{3} = \left[g-d+(p-cq)\left(1+\frac{w}{s}\right)\right]^2+4cqp\left(1+\frac{w}{s}\right)^2 \gt 0. $ |
It is easy to find that
$ \frac{g+d}{p+cq}-h = \frac{cqg^2+d^2p}{(p+cq)(dp+cqg)} \gt 0, $ |
so,
$ g+dp+cq>h. $
|
(4.6) |
Now, we discuss the stability of $ E_{1} $ from three cases.
Case 1. If $ h > 1+\frac{w}{s} $, then
$ σ=gd−(dp+cqg)(1+ws)=(dp+cqg)[gddp+cqg−(1+ws)]=(dp+cqg)[h−(1+ws)]>0, $
|
$ θ=(p+cq)(1+ws)−(g+d)=[(1+ws)−g+dp+cq](p+cq)<[h−g+dp+cq](p+cq). $
|
By (4.6) and $ \sigma > 0 $, we have $ \theta < 0 $ and $ \Delta_{3} = \theta^2-4\sigma < \theta^2 $, respectively. Therefore, both $ \lambda_{{3}, {4}} $ have negative real parts, $ E_{1} $ is stable accordingly.
Case 2. If $ h < 1+\frac{w}{s} $, then
$ \sigma = (dp+cqg)\left[h-\left(1+\frac{w}{s}\right)\right] \lt 0, $ |
and $ \Delta_{3} = \theta^2-4\sigma > \theta^2 $ accordingly. For any $ \theta $, (4.5) shows that one of $ \lambda_{{3}, {4}} $ is negative and the other is positive. Therefore, $ E_{1} $ is unstable.
Case 3. If $ h = 1+\frac{w}{s} $, then
$ \sigma = (dp+cqg)\left[h-\left(1+\frac{w}{s}\right)\right] = 0,\ \theta = \left[h-\frac{g+d}{p+cq}\right](p+cq). $ |
So, $ \Delta_{3} = \theta^2-4\sigma = \theta^2 $. Because (4.6) means $ \theta < 0 $, (4.5) shows that one of $ \lambda_{{3}, {4}} $ is $ \theta $ and the other is $ 0 $. Thus, the Jacobian matrix $ J(E_1) $ has one zero eigenvalue and three negative eigenvalues. By the center manifold theory, there exist a one-dimensional center manifold and a three-dimensional stable subspace at $ E_{1} $.
From the above discussion, we get the following results.
Theorem 4.3. Equilibrium $ E_{1} $ of model (1.5) is stable if $ h > 1+\frac{w}{s} $, and unstable if $ h < 1+\frac{w}{s} $.
By Theorems 4.2 and 4.3, we know that if $ h > 1+\frac{w}{s} $ then model (1.5) has no positive equilibrium and $ E_{1} $ is stable. Since the axes are invariant sets of model (1.5), boundary equilibrium $ E_1 $ is globally stable when $ h > 1+\frac{w}{s} $.
Theorem 4.4. If $ h > 1+\frac{w}{s} $, then the boundary equilibrium $ E_{1} $ is globally stable.
Theorem 4.5 gives a local center manifold at $ E_{1} $ for Case 3 of the above discussion.
Theorem 4.5. For boundary equilibrium $ E_{1} $, if $ h = 1+\frac{w}{s} $, then there is a local center manifold
$ W_{c}^{loc}(E_{1}) = \{(x_{2},y_{2},u_{2},v_{2})\in{R^4}:v_{2}\in{K_1},x_{2} = \beta_{1}v_{2}^2+0(v_{2}^3),y_{2} = \beta_{2}v_{2}^2+0(v_{2}^3),u_{2} = \beta_{3}v_{2}^2+0(v_{2}^3)\}, $ |
where $ K_1 $ is a small neighborhood of the origin,
$ [x2y2u2v2]=P−1[X−1Y−wsUV], P=[10i1i301i2i400p(d+θ)pd00q(g+θ)qg], $
|
and
$ \beta_{1} = \frac{s(pd+cqg)}{r(s+w)}\left[\left(\frac{s}{s+w}+\frac{i_{3}-i_{1}}{\theta}\right)\left(1+\frac{w}{s^2}\right)-1\right]i_{3}, $ |
$ \beta_{2} = \frac{pd+cqg}{r(s+w)}\left[\left(\frac{w}{s+w}+\frac{wi_{3}-s^2i_{2}}{s^2\theta}\right)\left(1+\frac{w}{s^2}\right)-\frac{w}{s^2}\right]i_{3}, $ |
$ \beta_{3} = -\frac{1}{\theta^2}\left(1+\frac{w}{s^2}\right)i_{3}(pd+cqg), $ |
with
$ i_{1} = -\frac{p(d+\theta)+cq(g+\theta)}{\theta+r\left(1+\frac{w}{s}\right)},\ i_{2} = -\frac{w}{s}\frac{p(d+\theta)+cq(g+\theta)}{\theta+rs\left(1+\frac{w}{s}\right)},\ i_{3} = -\frac{pd+cqg}{r\left(1+\frac{w}{s}\right)}. $ |
Proof. First, we take a matrix transformation. Let $ x_{1} = X-1 $, $ y_{1} = Y-\frac{w}{s} $, $ u_{1} = U $ and $ v_{1} = V $, then model (1.5) becomes
$ \left[ x′1y′1u′1v′1 \right] = J(E_1)\left[ x1y1u1v1 \right]+\left[ f1(x1,y1,u1,v1)f2(x1,y1,u1,v1)f3(x1,y1,u1,v1)f4(x1,y1,u1,v1) \right], $
|
where
$ f_{1}(x_{1},y_{1},u_{1},v_{1}) = x_{1}(-rx_{1}-ry_{1}-u_{1}-cv_{1}),\ f_{2}(x_{1},y_{1},u_{1},v_{1}) = y_{1}(-srx_{1}-sry_{1}-u_{1}-cv_{1}), $ |
$ f_{3}(x_{1},y_{1},u_{1},v_{1}) = p(x_{1}+y_{1})(u_{1}+cv_{x}),\ f_{4}(x_{1},y_{1},u_{1},v_{1}) = q(x_{1}+y_{1})(u_{1}+cv_{x}). $ |
The eigenvalues of $ J(E_1) $ are $ \lambda_{1} = -r\left(1+\frac{w}{s}\right) $, $ \lambda_{2} = -rs\left(1+\frac{w}{s}\right) $, $ \lambda_{3} = \theta $, $ \lambda_{4} = 0 $. By the center manifold theory, $ (0, 0, 0, 0) $ has a one-dimensional center subspace and a three-dimensional stable subspace. The eigenvectors of $ J(E_1) $ with respect to $ \lambda_{1} $, $ \lambda_{2} $, $ \lambda_{3} $ and $ \lambda_{4} $ are
$ e_{1} = \left[ 1000 \right],\ e_{2} = \left[ 0100 \right],\ e_{3} = \left[ i1i2p(d+θ)q(g+θ) \right],\ e_{4} = \left[ i3i4pdqg \right], $
|
respectively, where
$ i_{1} = -\frac{p(d+\theta)+cq(g+\theta)}{\theta+r\left(1+\frac{w}{s}\right)},\ i_{2} = -\frac{w}{s}\frac{p(d+\theta)+cq(g+\theta)}{\theta+rs\left(1+\frac{w}{s}\right)},\ i_{3} = -\frac{pd+cqg}{r\left(1+\frac{w}{s}\right)},\ i_{4} = \frac{w}{s^2}i_{3}. $ |
Let $ P = (e_{1}, e_{2}, e_{3}, e_{4}) $, then $ P^{-1}J(E_1)P = B $, where $ B = diag(\lambda_{1}, \lambda_{2}, \lambda_{3}, \lambda_{4}) $. Let
$ \left[ x1y1u1v1 \right] = P\left[ x2y2u2v2 \right] = \left[ x2+i1u2+i3v2y2+i2u2+i4v2p(d+θ)u2+pdv2q(g+θ)u2+qgv2 \right], $
|
then,
$ [x′2y′2u′2v′2]=P−1J(E1)P[x2y2u2v2]+P−1[f1(x2,y2,u2,v2)f2(x2,y2,u2,v2)f3(x2,y2,u2,v2)f4(x2,y2,u2,v2)]=B[x2y2u2v2]+[g1(x2,y2,u2,v2)g2(x2,y2,u2,v2)g3(x2,y2,u2,v2)g4(x2,y2,u2,v2)], $
|
where
$ g_{1}(x_{2},y_{2},u_{2},v_{2}) = f_{1}(x_{2},y_{2},u_{2},v_{2})+\frac{i_{3}-i_{1}}{p\theta}f_{3}(x_{2},y_{2},u_{2},v_{2}), $ |
$ g_{2}(x_{2},y_{2},u_{2},v_{2}) = f_{2}(x_{2},y_{2},u_{2},v_{2})+\frac{i_{4}-i_{2}}{p\theta}f_{3}(x_{2},y_{2},u_{2},v_{2}), $ |
$ g_{3}(x_{2},y_{2},u_{2},v_{2}) = \frac{1}{p\theta}f_{3}(x_{2},y_{2},u_{2},v_{2}),\ g_{4}(x_{2},y_{2},u_{2},v_{2}) = -\frac{1}{p\theta}f_{3}(x_{2},y_{2},u_{2},v_{2}), $ |
and
$ f_{1}(x_{2},y_{2},u_{2},v_{2}) = (x_{2}+i_{1}u_{2}+i_{3}v_{2})[-rx_{2}-ry_{2}-(rm_{1}+m_{3}+m_{4})u_{2}-(rm_{2}+m_{3})v_{2}], $ |
$ f_{2}(x_{2},y_{2},u_{2},v_{2}) = (y_{2}+i_{2}u_{2}+i_{4}v_{2})[-srx_{2}-sry_{2}-(srm_{1}+m_{3}+m_{4})u_{2}-(srm_{2}+m_{3})v_{2}], $ |
$ f_{3}(x_{2},y_{2},u_{2},v_{2}) = p\left[x_{2}+y_{2}+(i_{1}+i_{2})u_{2}+(i_{3}+i_{4})v_{2}\right]\left[(m_{3}+m_{4})u_{2}+m_{3}v_{2}\right], $ |
$ f_{4}(x_{2},y_{2},u_{2},v_{2}) = q\left[x_{2}+y_{2}+(i_{1}+i_{2})u_{2}+(i_{3}+i_{4})v_{2}\right]\left[(m_{3}+m_{4})u_{2}+m_{3}v_{2}\right], $ |
with
$ m_{1} = i_{1}+i_{2},\ m_{2} = i_{3}+i_{4},\ m_{3} = pd+cqg,\ m_{4} = d\theta+cq\theta. $ |
The local center manifold is described by a approximation function $ \eta:span{\{e_{4}\}}\rightarrow{span{\{e_{1}, e_{2}, e_{3}\}}}, $ which we set
$ \eta(v_{2})e_{4} = \eta_{1}(v_{2})e_{1}+\eta_{2}(v_{2})e_{2}+\eta_{3}(v_{2})e_{3}. $ |
Since $ \eta(0) = D_{\eta}(0) = 0 $, the Taylor expansion at the origin is
$ x_{2} = \eta_{1}(v_{2}) = \beta_{1}v_{2}^2+\cdots,\ y_{2} = \eta_{2}(v_{2}) = \beta_{2}v_{2}^2+\cdots,\ u_{2} = \eta_{3}(v_{2}) = \beta_{3}v_{2}^2+\cdots. $ |
By the center manifold approximation theory, we have
$ -r(1+\frac{w}{s})x_{2}+g_{1} = (2\beta_{1}v_{2}+\ldots)g_{4}, $ |
$ -rs(1+\frac{w}{s})y_{2}+g_{2} = (2\beta_{2}v_{2}+\ldots)g_{4}, $ |
$ \theta{u_{2}}+g_{3} = (2\beta_{3}v_{2}+\ldots)g_{4}. $ |
Simplifying the above equations and comparing the coefficients of the the same order terms in both sides, we get
$ \beta_{1} = \frac{s(pd+cqg)}{r(s+w)}\left[\left(\frac{s}{s+w}+\frac{i_{3}-i_{1}}{\theta}\right)\left(1+\frac{w}{s^2}\right)-1\right]i_{3}, $ |
$ \beta_{2} = \frac{pd+cqg}{r(s+w)}\left[\left(\frac{w}{s+w}+\frac{wi_{3}-s^2i_{2}}{s^2\theta}\right)\left(1+\frac{w}{s^2}\right)-\frac{w}{s^2}\right]i_{3}, $ |
$ \beta_{3} = -\frac{1}{\theta^2}(i_{3}+i_{4})(pd+cqg). $ |
Thus, we get the center manifold as follow
$ W_{c}^{loc}(E_{1}) = \{(x_{2},y_{2},u_{2},v_{2})\in{R^4}:v_{2}\in{K_1},x_{2} = \beta_{1}v_{2}^2+O(v_{2}^3),y_{2} = \beta_{2}v_{2}^2+O(v_{2}^3),u_{2} = \beta_{3}v_{2}^2+O(v_{2}^3)\}, $ |
where $ K_1 $ is a small neighborhood of the origin, and
$ [x2y2u2v2]=P−1[X−1Y−wsUV]. $
|
In this section, we consider the stability of $ E_{*} $.
Theorem 4.6. If $ h < 1+\frac{w}{s} $, then the unique positive equilibrium $ E_{*} $ of model (1.5) is stable.
Proof. The Jacobian matrix of model (1.5) at $ E_{*} $ is
$ J_{*} = \left[ r−2rX∗−rY∗−U∗−cV∗r−rX∗−X∗−cX∗rw−rsY∗rw−rsX∗−2rsY∗−U∗−cV∗−Y∗−cY∗p(U∗+cV∗)p(U∗+cV∗)p(X∗+Y∗)−gcp(X∗+Y∗)q(U∗+cV∗)q(U∗+cV∗)q(X∗+Y∗)cq(X∗+Y∗)−d \right]. $
|
Let $ j = \frac{cqg^2+pd^2}{pd+cqg} $, the characteristic polynomial of $ J_{*} $ is
$ f∗(λ)=λ4+a3λ3+a2λ2+a1λ+a0, $
|
(4.7) |
where
$ a_{3} = r(s+1)(X_{*}+Y_{*})+(U_{*}+cV_{*})+j, $
$ a_{2} = r^2s(X_{*}+Y_{*})^2+r^2(X_{*}+Y_{*})[s(w-sY_{*})+(1-X_{*})]+(d+g)(U_{*}+cV_{*})+r(s+1)j(X_{*}+Y_{*}), $
$ a1=r(p+cq)(sX∗+Y∗)(X∗+Y∗)(U+cV∗)+r2(X∗+Y∗)j[s(w−sY∗+(1−X∗))]+gd(U+cV∗)+(p+cq)(X∗+Y∗)(U+cV∗)2+r2sj(X∗+Y∗)2,
$ a_{0} = gd[(U+cV_{*})+r(sX_{*}+Y_{*})](U_{*}+cV_{*}). $
For $ 0 < X_{*} < 1, $ $ 0 < Y_{*} < \frac{w}{s} $, we know $ a_{3} > 0 $, $ a_{2} > 0 $, $ a_{1} > 0 $ and $ a_{0} > 0 $. According to the Routh-Hurwitz stability criterion, if
$ a2a3−a1>0 $
|
(4.8) |
also holds, then positive equilibrium $ E_{*} $ is stable.
Next, we prove that (4.8) is satisfied. Let $ g(r) = a_{2}a_{3}-a_{1} $. Substituting the values of $ U_{*} $ and $ V_{*} $ into $ g(r) $, we get that $ g(r) = b_{2}r^3+b_{1}r^2+b_{0}r $, where
$ b_{2} = (X_{*}+Y_{*})(sX_{*}+Y_{*}+w+1)\left[s(w-sY_{*})+(1-X_{*})+s(X_{*}+Y_{*})\right], $
$ b_{1} = (s+1)^2(X_{*}+Y_{*})^2j+(w-sY_{*}+1-X_{*})(2sX_{*}+2Y_{*}+1+w)j, $
$ b_{0} = (s+1)(X_{*}+Y_{*})j^2+(w-sY_{*}+1-X_{*})[(d+g)j-dg]. $
From (4.2), we know that $ X_{*} $ and $ Y_{*} $ do not depend on $ r $, and $ b_{2} > 0 $, $ b_{1} > 0 $. Since
$ (d+g)j-dg = \frac{pd^3+cqg^3}{pd+cqg} \gt 0, $ |
we have $ b_{0} > 0 $. So, $ g(r) > 0 $, i.e., (4.8) holds. Therefore, (4.8) and Theorem 4.1 mean that $ E_* $ is stable.
Remark 4.1. Theorem 4.1 and 4.6 mean $ h < 1+\frac{w}{s} $ is the necessary and sufficient condition for the existence and stability of the unique positive equilibrium $ E_{*} $ of model (1.5).
Remark 4.2. By Theorems 4.1, 4.3 and 4.6, we know that if $ h > 1+\frac{w}{s} $ then model (1.5) has no positive equilibrium and $ E_{1} $ is stable, and if $ h < 1+\frac{w}{s} $ then the unique positive equilibrium of model (1.5) is stable and $ E_{1} $ is unstable. Therefore, there is a transcritical bifurcation in model (1.5).
In this section, some examples are given to illustrate our main results. Furthermore, we get some strategies for species protection by insight into the examples.
From the discussion in Section 4, we know that the boundary equilibrium $ E_{1} $ is stable when $ h > 1+\frac{w}{s} $, and the positive equilibrium $ E_{*} $ appears and becomes stable when $ h < 1+\frac{w}{s} $. In order to testify our analyzed results, we take some examples by simulation. For convenience, fix
$ r = 1, c = 1/4, w = 1/2, p = 1/3, q = 2/7, g = 3/4, d = 1/4, $ |
and the initial condition $ (0.3, 0.5, 0.4, 0.8). $ Next, we observe the time series diagrams of the four populations with four different values of $ s $ (see Figure 1). In this figure, the red dotted line represents the size of prey $ X $, the blue dotted line represents the size of prey $ Y $, the red solid line represents the size of predator $ U $, the blue solid line represents the size of predator $ V $.
In Figure 1(a), we choose $ s = 27/17 $ which satisfies the equality $ h > 1+\frac{w}{s} $. According to Theorem 4.3 and Remark 4.1, in this case the boundary equilibrium $ E_{1}(1, \frac{17}{54}, 0, 0) $ is stable while the positive equilibrium doesn't exist. In the figure, the size of two prey populations $ X $ and $ Y $ approach to 1 and $ 0.31\approx\frac{17}{54} $ respectively as the time increases, while the two predator populations $ U $ and $ V $ both tend to zero. This shows that $ E_{1} $ is stable, which verify our results.
In Figure 1(b), we choose $ s = 1.1 $ which satisfies the equality $ h < 1+\frac{w}{s} $. According to Theorem 4.3 and Remark 4.1, in this case the boundary equilibrium $ E_{1}(1, \frac{5}{11}, 0, 0) $ is unstable while the positive equilibrium $ E_{*}(0.96, 0.41, 0.04, 0.09) $ is stable. In the figure, the size of two prey populations $ X $ and $ Y $ approach to 0.96 and 0.41, and the two predator populations $ U $ and $ V $ both tend to 0.04 and 0.09 respectively as the time increases. This shows that $ E_{*} $ is stable, which verify our results.
In Figure 1(c), we choose $ s = 1 $ which satisfies the equality $ h < 1+\frac{w}{s} $. According to Theorem 4.3 and Remark 4.1, in this case the boundary equilibrium $ E_{1}(1, \frac{1}{2}, 0, 0) $ is unstable while the positive equilibrium $ E_{*}(0.91, 0.46, 0.08, 0.20) $ is stable. In the figure, the size of two prey populations $ X $ and $ Y $ approach to 0.91 and 0.46, and the two predator populations $ U $ and $ V $ both tend to 0.08 and 0.20 respectively as the time increases. This verifies our results.
In Figure 1(d), we choose $ s = 0.8 $ which satisfies the equality $ h < 1+\frac{w}{s} $. According to Theorem 4.3 and Remark 4.1, in this case the boundary equilibrium $ E_{1}(1, \frac{5}{8}, 0, 0) $ is unstable while the positive equilibrium $ E_{*}(0.86, 0.51, 0.14, 0.35) $ is stable. In the figure, the size of two prey populations $ X $ and $ Y $ approach to 0.86 and 0.51, and the two predator populations $ U $ and $ V $ both tend to 0.14 and 0.35 respectively as the time increases. This verifies our results.
From the biological point of view, the decrease of $ s $ can be regarded as the decrease in the ratio of the survival pressure of population $ Y $ to the survival pressure of population $ X $. In the process of decreasing, the population $ Y $ has more advantages in survival than population $ X $, which leads to an increase in the population density. Conversely, the density of population $ X $ decreases due to the stronger stress in the same environment. At the same time, both the densities of predator $ U $ and $ V $ are increased. It is worth noting that, with the change of $ s $ in the above examples, the total amount of two prey populations remains unchanged, while the total density of two predator populations expands. This result gives us two implications for model (1.5): First, the balance of two prey groups can be adjusted by controlling the environmental pressure experienced by one of the prey population; Second, endangered predator populations can be saved by controlling the survival pressure of the prey population.
In this paper, a predator-prey model with genetic differentiation both happened in the predator and prey is studied, and some interesting results are obtained.
For model (1.5) without the predators, by the qualitative analysis and the Bendixson-Dulac criterion, we get that the boundary equilibrium $ (0, 0) $ is unstable, and the positive equilibrium $ (1, \frac{w}{s}) $ is global asymptotical stability. These indicate that the two genotypes of the prey will be permanent.
Model (1.5) with one genotype of the predator and prey has two boundary equilibria $ (0, 0) $ and $ (1, 0) $, and a positive equilibrium $ \left(\frac{g}{p}, \frac{r(p-g)}{p}\right) $ when $ p > g $. The boundary equilibrium $ (0, 0) $ is always unstable. The boundary equilibrium $ (1, 0) $ is stable when $ p < g $, unstable when $ p > g $, and there is a stable local center manifold at $ (1, 0) $ when $ p = g $. The positive equilibrium $ \left(\frac{g}{p}, \frac{r(p-g)}{p}\right) $ is always stable if it exists. We also show that there is no periodic solution in the model by using Bendixson-Dulac criterion. Therefore, the global dynamics of the model is clear as follows: boundary equilibrium $ (1, 0) $ is globally stable when $ p < g $, positive equilibrium $ \left(\frac{g}{p}, \frac{r(p-g)}{p}\right) $ is globally stable when $ p > g $, and there is a transcritical bifurcation in the model when $ p = g $.
Model (1.5) has two boundary equilibria $ E_{0}(0, 0, 0, 0), $ $ E_{1}(1, \frac{w}{s}, 0, 0) $, and a positive equilibrium $ E_{*} $ if and only if $ h < 1+\frac{w}{s} $. The boundary equilibrium $ E_{0} $ is always unstable. Another boundary equilibrium $ E_{1} $ is stable when $ h > 1+\frac{w}{s} $, unstable when $ h < 1+\frac{w}{s} $, and there is a stable local center manifold at $ (1, 0) $ when $ h = 1+\frac{w}{s} $. The positive equilibrium $ E_{*} $ is stable if it exists. Furthermore, there is a transcritical bifurcation in the model when $ h = 1+\frac{w}{s} $.
If $ h > 1+\frac{w}{s} $, then model (1.5) has no positive equilibrium, an unstable $ E_0 $ and a stable $ E_1 $. Since the axes are invariant with respect to model (1.5), boundary equilibrium $ E_1 $ is globally stable when $ h > 1+\frac{w}{s} $. This means that the prey will be permanent and the predators will be extinct, which is similar to the principle of competitive exclusion for two competitive species [45]. However, if $ h < 1+\frac{w}{s} $ then model (1.5) has a local positive equilibrium and two unstable boundary equilibria. The global stability of the positive equilibrium needs to be further analyzed.
The authors would like to thank the referees for their valuable comments and helpful suggestions.
The research was supported by the Visiting Scholar Program of Ludong University. The corresponding author would like to thank Department of Mathematics and Statistics at Memorial University of Newfoundland for kindly hosting during his visiting.
The authors declare that they have no conflict of interest.
1. | Jiandong Zhao, Tonghua Zhang, Zhixia Han, Dynamic behavior of a stochastic SIRS model with two viruses, 2020, 0, 2191-0294, 10.1515/ijnsns-2019-0208 | |
2. | Eric Planas, Jordi Mill, Andy L. Olivares, Xabier Morales, Maria Isabel Pons, Xavier Iriart, Hubert Cochet, Oscar Camara, 2022, Chapter 18, 978-3-030-93721-8, 160, 10.1007/978-3-030-93722-5_18 | |
3. | G. P. Neverova, O. L. Zhdanova, E. Ya. Frisman, Evolutionary dynamics of predator in a community of interacting species, 2022, 108, 0924-090X, 4557, 10.1007/s11071-022-07372-z | |
4. | O. L. Zhdanova, G. P. Neverova, E. Ya. Frisman, Predator Evolution in a Model of Interacting Species: To the Question about Maintaining Polymorphism by Litter Size in Natural Populations of Arctic Fox, 2022, 58, 1022-7954, 94, 10.1134/S1022795422010136 | |
5. | Li Li, Dongshen Fang, Qiyao Ye, Tan Hu, Shaobo Shi, An ensemble framework for risk prediction of left atrial thrombus based on undersampling with replacement, 2024, 36, 0941-0643, 18613, 10.1007/s00521-024-10166-6 | |
6. | G. Santhosh Kumar, C. Gunasundari, Salah Mahmoud Boulaaras, M. Aakash, N.B. Sharmila, Genetic differentiation of prey predator interaction model along with an Holling type-II functional response, 2024, 9, 26668181, 100649, 10.1016/j.padiff.2024.100649 |
Example 1. Convergence tests for the spatial (left) and temporal (right) discretisation via mixed
Example 2: snapshots at
Example 3A: snapshots of FitzHugh-Nagumo dynamics on a porous mixture at early (left) and advanced (right) times. Computed solutions from top to bottom: membrane voltage, vorticity, and velocity.
Example 3A: Number of inner Newton steps and outer Picard steps needed to reach residual convergence to a tolerance of 1e-6.
Approximate membrane voltage, velocity, and pressure for the FitzHugh-Nagumo dynamics on a porous mixture at early (top), moderate (middle row), and advanced (bottom panels) times.
Example 4: Computed solutions (cytosolic calcium, sarcoplasmic calcium, vorticity, velocity, and pressure) for the intracellular calcium dynamics at early (left) and advanced (right) times.