Transmission dynamics of breast cancer through Caputo Fabrizio fractional derivative operator with real data

1.   Introduction

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

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

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

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

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

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

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.

2.   Stability of model (1.5) without the predators

Without the predators, model (1.5) becomes

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

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

it is easy to know that $B_1(X, Y)$ is a continuously differentiable function in the first quadrant. We get that

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.

3.   Stability of model (1.5) with one genotype in both the predator and prey

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

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

when $p = g$, where $K$ is a small neighborhood of the origin, and

Proof. The Jacobian matrix of (3.1) is

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

where

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

respectively. Let $P = (e_{1}, e_{2})$, then $P^{-1}AP = B$, where $B = {diag}(\lambda_{1}, \lambda_{2})$. Therefore, Let

then

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

By the center manifold approximation theorem, we have

Simplifying the above equation and comparing the coefficients of the two sides, we get

Thus, the center manifold at $(0, 0)$ is

where $K$ is a small neighborhood of the origin, and

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

It is easy to know that $B_2(X, U)$ is a continuous differentiable function in the first quadrant. Then, in the first quadrant

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.

4.   The dynamics of model (1.5)

4.1.   Existence of the equilibria

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

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

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

Let $\Delta_{2} = [w+1-(1-s)h]^2+4wh(1-s),$ then $\Delta_{2} = [w-1+(1-s)h]^2+4w > 0$. Let

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),

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

By (4.2) and the above equation, we get

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

which means

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$.

4.2.   Stability of the boundary equilibria

The Jacobian matrix of model (1.5) at $E_{0}$ is

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

The characteristic polynomial of $J(E_1)$ is

$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

The two roots of $f_1(\lambda)$ are

Let $\Delta_{3} = \theta^2-4\sigma,$ then

It is easy to find that

so,

Now, we discuss the stability of $E_{1}$ from three cases.

Case 1. If $h > 1+\frac{w}{s}$, then

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

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

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

where $K_1$ is a small neighborhood of the origin,

and

with

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

where

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

respectively, where

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

then,

where

and

with

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

Since $\eta(0) = D_{\eta}(0) = 0$, the Taylor expansion at the origin is

By the center manifold approximation theory, we have

Simplifying the above equations and comparing the coefficients of the the same order terms in both sides, we get

Thus, we get the center manifold as follow

where $K_1$ is a small neighborhood of the origin, and

4.3.   Stability of the positive equilibrium

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

Let $j = \frac{cqg^2+pd^2}{pd+cqg}$, the characteristic polynomial of $J_{*}$ is

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_{*}),$

$\begin{aligned} a_{1} = &r(p+cq)(sX_{*}+Y_{*})(X_{*}+Y_{*})(U+cV_{*})+r^2(X_{*}+Y_{*})j[s(w-sY_{*}+(1-X_{*}))]\\ & +gd(U+cV_{*})+(p+cq)(X_{*}+Y_{*})(U+cV_{*})^2+r^2sj(X_{*}+Y_{*})^2, \end{aligned}$

$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

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

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).

5.   Examples

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

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.

6.   Conclusions and discussions

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.

Acknowledgments

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.

Conflict of interest

The authors declare that they have no conflict of interest.