Stability Analysis of a Self-Memory Prey-Predator Diffusion Model Based on Bazykin Functional Response

Yanzhe Han, Fuqin Sun^*

School of Sciences, Tianjin University of Technology and Education, Tianjin, China.

DOI: 10.4236/jamp.2023.115090                

Abstract

To investigate the effects of self-memory diffusion on predator-prey models, we consider a predator-prey model with Bazykin functional response of self- memory diffusion. The uniqueness, boundedness, positivity, existence and stability of equilibrium point of the model are studied. In this paper, the uniqueness of the solution is discussed under the non-negative initial function and Neumann boundary conditions satisfying a specific space. The boundness of the solution is proved by the comparison principle of parabolic equations, and the positivity of the solution is proved by the strong maximum principle of parabolic equations. Hurwitz criterion and Lyapunov function construction are used to analyze the local stability and global stability of feasible equilibrium points. The results show that the system solution is unique non-negative and bounded. The model is unstable at the trivial equilibrium point E0 and the boundary equilibrium point E1, and the condition of whether the positive equilibrium point E2 is stable under certain conditions is given.

Keywords

Bazykin Functional Response, Lyapunov Function, Boundedness, Uniqueness, Stability

1. Introduction

Predator-prey interactions in nature are the main cause of the rich diversity of complex ecosystems. Therefore, the qualitative and quantitative analysis of predator-prey relationships is of great theoretical and practical significance, and is an important field of population biology. Originally, Lotka and Volterra [1] proposed the same model describing predator-prey interactions, which was named the Lotka Volterra model. Later, Holling [2] [3] proposed different functional response functions to improve the model and describe the dynamic behavior of predator and prey more accurately. However, the spatial spread of predators and captive food species may further complicate the spatio-temporal dynamics. So now in the spatial predator-prey model, it’s usually assumed that the predator and the prey plant move randomly through their habitat, which is simulated by the diffusion equation. Relevant studies can be seen [4] [5] , etc. In addition, some intelligent predators also have memory effect and cognitive behavior. For example, blue whales migrate by memory. In reference [6] , a model was proposed to describe the motion of a single population with spatial memory, and the research results indicate that memory based diffusion may have a significant impact on the distribution of the population. Subsequently, reference [7] incorporated memory based diffusion into the classical diffusion model to describe the interaction between two species. At this point, there are usually two types of memory and cognition: diffusion based on cross memory and diffusion based on self memory. Reference [8] studied a predator-prey model based on cross memory diffusion, but there is currently little research on self-memory diffusion predator-prey models.

In the existing studies on biological predation diffusion system, Holling I, Holling II and Holling III functional responses have been widely studied, while Bazykin functional responses have been less studied. Bazykin functional response can describe the destabilizing force of predator saturation and the stabilizing force of prey competition, so the study of Bazykin functional response has more practical significance. Therefore, a self-memory diffusion model of predation based on Bazykin functional response is established in this paper to study the properties of its solutions and the existence of equilibrium points, and to analyze the local stability of equilibrium points.

2. Model Composition

Consider the predator’s diffusion-reaction-diffusion model based on self-memory [9] as follows:

$\{\begin{array}{l}{u}_t=d_1{u}_{xx}+ru\left(1-u/K\right)-puv,\\ v_t=d_2{v}_{xx}+d_3{\left(v{v}_x\left(x,t-\tau \right)\right)}_x+cpuv-\theta v.\end{array}$ (1)

where, $u=u\left(x,t\right)$ represents the population density of prey at time; $v=v\left(x,t\right)$ Represents the population density of predators at time; r and $\theta$ represent the prey birth rate and predator death rate, respectively; K represents environmental bearing rate; The conversion rate $c\left(0<c<1\right)$ of prey to predator; $d_1,d_2,d_3$ represents the spatial diffusion term, in particular, $d_3>0$ represents the movement of a predator from a high density location to a low density location.

The predator’s functional response follows the Holling I functional response $puv$ . $1/\left(1+\alpha u\right)\left(1+\beta v\right)$ is called Bazykin-type functional response function, where $\alpha ,\beta$ are two normal numbers. Bazykin functional response can describe the destabilizing force of predator saturation and the stabilizing force of prey competition, so the study of Bazykin functional response has more practical significance. Therefore, in this paper, Holling I functional response is replaced by Bazykin functional response to study the one-dimensional spatial domain $\Omega =\left(0,l\pi \right),l\in R^{+}$ of the system. Considering the above assumptions, and after a certain scale transformation, the following model can be obtained

$\{\begin{array}{ll}{u}_t=d_1{u}_{xx}+u\left(1-u\right)-\frac{uv}{\left(1+\alpha u\right)\left(1+\beta v\right)},\hfill & x\in \left(0,l\pi \right),t>0,\hfill \\ v_t=d_2{v}_{xx}+d_3{\left(v{v}_x\left(x,t-\tau \right)\right)}_x-\frac{c{v}^2}{\left(1+\alpha u\right)\left(1+\beta v\right)}+\gamma v,\hfill & x\in \left(0,l\pi \right),t>0,\hfill \\ u_x\left(x,t\right)=v_x\left(x,t\right)=0,\hfill & x=0,l\pi ,t>0,\hfill \\ u\left(x,0\right)=u\left(x\right)\ge 0,v\left(x,t\right)=\phi \left(x,t\right)\ge 0,\hfill & x\in \left[0,l\pi \right],-\tau \le t\le 0.\hfill \end{array}$ (2)

where x denotes the position of the predator or prey at the moment t. None of the above parameters are negative.

3. Correlation Suitability of Global Solutions of System (2) When $\tau >0$

When $\tau >0$ , Let the initial value function $u_0\left(x\right),\phi \left(x,t\right)$ satisfy

$u_0\left(x\right)\in C^{2+\alpha }\left(\left[0,l\pi \right]\right),\phi \left(x,t\right)\in C^{2+\alpha ,1+\alpha /2}\left(\left[0,l\pi \right]\times \left[-\tau ,0\right]\right).$ (3)

Theorem 1 assumes that all parameters of (2) are positive and satisfies (2) for any non-negative initial value function $u_0\left(x\right)\ge 0,v_0\left(x\right)\ge 0$ . For (2) there is a unique global solution $\left(u\left(x,t\right),v\left(x,t\right)\right)$ , and $u\left(x,t\right)\ge 0,v\left(x,t\right)\ge 0$ ( $\cancel{\equiv }0$ ). If satisfied $\beta \left(1+\alpha \right)<c/\gamma$ , then $u\left(x,t\right),v\left(x,t\right)$ satisfied

$\underset{t\to \infty }{\lim }\sup \underset{\bar{\Omega }}{\max }u\left(x,t\right)\le 1,\underset{t\to \infty }{\lim }\sup \underset{\bar{\Omega }}{\max }v\left(x,t\right)\le \frac{\gamma \left(1+\alpha \right)}{c-\beta \gamma \left(1+\alpha \right)}.$

Proof: (I) For $t\in \left[0,\tau \right]$ , The initial value function for predator density $v\left(x,t\right)$ can be expressed as $v_{\tau }=\phi \left(x,t-\tau \right)$ .

First consider the following initial value boundary problem

$\{\begin{array}{ll}{u}_t=d_1{u}_{xx}+u\left(1-u\right)-\frac{uv}{\left(1+\alpha u\right)\left(1+\beta v\right)},\hfill & x\in \left(0,l\pi \right),t>0,\hfill \\ v_t=d_2{v}_{xx}+d_3{\left(v{\phi }_x\left(x,t-\tau \right)\right)}_x-c\frac{v^2}{\left(1+\alpha u\right)\left(1+\beta v\right)}+\gamma v,\hfill & x\in \left(0,l\pi \right),t>0,\hfill \\ u_x\left(x,t\right)=0,v_x\left(x,t\right)=0,\hfill & x=0,l\pi ,t>0,\hfill \\ u\left(x,t\right)=u_0\left(x\right),v\left(x,t\right)=\phi \left(x,0\right),\hfill & x\in \left[0,l\pi \right].\hfill \end{array}$ (4)

where $u\left(x,t\right)$ satisfy

$\{\begin{array}{ll}{u}_t=d_1{u}_{xx}+u\left(1-u\right)-\frac{uv}{\left(1+\alpha u\right)\left(1+\beta v\right)}\le d_1{u}_{xx}+u\left(1-u\right),\hfill & x\in \left(0,l\pi \right),t>0,\hfill \\ u_x\left(x,t\right)=0,\hfill & x=0,l\pi ,t>0.\hfill \end{array}$

Let $u^1\left(t\right)$ be the solution of the following system (5)

$\{\begin{array}{ll}{u}_t=d_1{u}_{xx}+u\left(1-u\right),\hfill & x\in \left(0,l\pi \right),t>0，\hfill \\ u_x\left(x,t\right)=0,\hfill & x=0,l\pi ,t>0.\hfill \end{array}$ (5)

Then $t\to \infty$ , $u^1\left(t\right)\to 1$ . Any given $\epsilon$ satisfaction the following formula can be obtained from the comparison principle of parabolic equations

$0<\epsilon <\frac{c-\gamma \beta \left(1+\alpha \right)}{\gamma \alpha \beta },$

exist $T_0=T_0\left(\epsilon \right)\gg 1$ , s.t

$u\left(x,t\right)\le 1+\epsilon ,\left(x,t\right)\in \bar{\Omega }\times \left[T_0,\infty \right).$

i.e.

$\underset{t\to \infty }{\lim }\sup \underset{\bar{\Omega }}{\max }u\left(x,t\right)\le 1.$

Let $v^1\left(t\right)$ be the solution of the following system (6)

$\{\begin{array}{ll}{v}_t=d_2{v}_{xx}+d_3{\left(v{\phi }_x\left(x,t-\tau \right)\right)}_x+v\frac{\gamma \left(1+\alpha \left(1+\epsilon \right)\right)-\left[c-\beta \gamma \left(1+\alpha \left(1+\epsilon \right)\right)\right]v}{\left[1+\alpha \left(1+\epsilon \right)\right]\left(1+\beta v\right)},\hfill & x\in \left(0,l\pi \right),t>0,\hfill \\ v_x\left(x,t\right)=0.\hfill & x=0,l\pi ,t>0.\hfill \end{array}$ (6)

when $t\to +\infty$ , there is

$v^1\left(t\right)\to \frac{\gamma \left[1+\alpha \left(1+\epsilon \right)\right]}{c-\beta \gamma \left(1+\alpha \left(1+\epsilon \right)\right)}.$

Due to

$\begin{array}{l}\gamma v-c\frac{v^2}{\left(1+\alpha u\right)\left(1+\beta v\right)}\\ \le v\frac{\gamma \left(1+\alpha \left(1+\epsilon \right)\right)-\left[c-\beta \gamma \left(1+\alpha \left(1+\epsilon \right)\right)\right]v}{\left[1+\alpha \left(1+\epsilon \right)\right]\left(1+\beta v\right)},\left(x,t\right)\in \bar{\Omega }\times \left[T_0,\infty \right).\end{array}$

According to the principle of comparison,

$v\left(t\right)\le \frac{\gamma \left[1+\alpha \left(1+\epsilon \right)\right]}{c-\beta \gamma \left(1+\alpha \left(1+\epsilon \right)\right)}+\epsilon ,\left(x,t\right)\in \bar{\Omega }\times \left[T_0,\infty \right).$

I.e.

$\underset{t\to \infty }{\lim }\sup \underset{\bar{\Omega }}{\max }v\left(x,t\right)\le \frac{\gamma \left[1+\alpha \left(1+\epsilon \right)\right]}{c-\beta \gamma \left(1+\alpha \left(1+\epsilon \right)\right)}.$

By the arbitrariness of $\epsilon$ , there is

$\underset{t\to \infty }{\lim }\sup \underset{\bar{\Omega }}{\max }v\left(x,t\right)\le \frac{\gamma \left(1+\alpha \right)}{c-\beta \gamma \left(1+\alpha \right)}.$

That is, $\left(u,v\right)$ is $L^1$ bounded on $\left[0,\tau \right]$ . According to Theorem 3 of literature [10] , we can know that the system has a unique solution on $\left[0,\tau \right]$ . According to the strong maximum principle of parabolic equation: $u\left(x,t\right)>0$ , $v\left(x,t\right)>0$ , $\left(x,t\right)\in \left[0,l\pi \right]\times \left(0,\tau \right)$ . In the same way, for the analysis of $t\in \left[\tau ,2\tau \right),t\in \left[2\tau ,3\tau \right),\cdots$ is consistent with the results of $t\in \left(0,\tau \right)$ . Finally

$\left(u\left(x,t\right),v\left(x,t\right)\right)>0,\left(x,t\right)\in \left[0,l\pi \right]\times \left(0,\infty \right).$

4. Analysis of System (2) about the Existence and Stability of Steady-State Solutions

4.1. Existence of Equilibrium Point

The equilibrium point of the system can be obtained by solving $u_t=0,v_t=0$ . There are obviously three equilibrium points: trivial equilibrium point $E_0=\left(0,0\right)$ , boundary equilibrium point $E_1=\left(1,0\right)$ and positive equilibrium point $E_2=\left(u_{*},v_{*}\right)$ . The simultaneous equations can be solved

$u_{*}=1-\frac{\gamma }{c},v_{*}=\frac{{\alpha }^{*}}{c^2-{\alpha }^{*}\beta }.$

where ${\alpha }^{*}=\left[\left(c-\gamma \right)\alpha +c\right]\gamma$ .

If $E_2$ exist, then satisfy

$\{\begin{array}{l}1-\gamma /c>0,\\ {\alpha }^{*}/\left(c^2-{\alpha }^{*}\beta \right)>0.\end{array}$

Solve the set of inequalities can get

(H₀) $0<\gamma <c,0<\beta <\frac{4\alpha }{{\left(\alpha +1\right)}^2}.$

4.2. Stability Analysis of E₀, E₁

Define 1. Real valued Sobolev space

${\rm X}=\left\{\left(u,v\right)\in H^2\left(0,l\pi \right)\times H^2\left(0,l\pi \right):{\left(u_x,v_x\right)|}_{x=0,l\pi }=0\right\}.$

For arbitrary $U_1,U_2\in {\rm X}$ , Define the inner product on space X to be

$\langle U_1,U_2\rangle ={\displaystyle {\int }_0^{l\pi }{U}_1^{\text{T}}{U}_2\text{d}x},U_1,U_2\in {\rm X}.$

Its eigenvalue problem is

$\{\begin{array}{l}-u^{″}=\mu u,x\in \left(0,l\pi \right),\hfill \\ u^{\prime }\left(0\right)=u^{\prime }\left(l\pi \right)=0.\hfill \end{array}$ (7)

The eigenvalues and eigenfunctions corresponding to the Equations (7) are

${\mu }_n={\left(n/l\right)}^2,{\eta }_n\left(x\right)=\cos \left(nx/l\right),n=0,1,\cdots .$

Theorem 2. The following conclusions for system (2) are valid: $E_0=\left(0,0\right)$ is unstable; $E_1=\left(1,0\right)$ is unstable.

Proof: Let

$f\left(u,v\right)=u\left(1-u\right)-\frac{uv}{\left(1+\alpha u\right)\left(1+\beta v\right)},g\left(u,v\right)=-c\frac{v^2}{\left(1+\alpha u\right)\left(1+\beta v\right)}+\gamma v.$

The Jacobian matrix of system (2) can be obtained at any point $J=\left(a_{ij}\right)\in R^{2\times 2}$ .

Where

$\begin{array}{l}{a}_{11}=f_u\left(u,v\right)=1-2u-\frac{v}{\left(1+\beta v\right){\left(1+\alpha u\right)}^2},\\ a_{12}=f_v\left(u,v\right)=-\frac{u}{\left(1+\alpha u\right){\left(1+\beta v\right)}^2},\\ a_{21}=g_u\left(u,v\right)=\frac{c\alpha v^2}{\left(1+\beta v\right){\left(1+\alpha u\right)}^2},\\ a_{22}=g_v\left(u,v\right)=\frac{-c(2v+\beta v^2)}{\left(1+\alpha u\right){\left(1+\beta v\right)}^2}+\gamma .\end{array}$

Linearize system (2) at $E_0$

$\left(\begin{array}{c}{u}_t\\ v_t\end{array}\right)=d\Delta \left(\begin{array}{c}u\\ v\end{array}\right)+A\left(\begin{array}{c}u\\ v\end{array}\right),$

where $d\Delta =\left(\begin{array}{cc}{d}_1\Delta & 0\\ 0& d_2\Delta \end{array}\right)$ , $A=J\left(0,0\right)=\left(\begin{array}{cc}1& 0\\ 0& \gamma \end{array}\right)$ .

It is known that the corresponding eigenvalue is ${\mu }_n={\left(n/l\right)}^2$ for $-\Delta$ , so the corresponding eigenvalue is ${\sigma }_n=-d_1{\left(n/l\right)}^2$ for $d_1\Delta$ , the corresponding eigenvalue is ${\epsilon }_n=-d_2{\left(n/l\right)}^2$ for $d_2\Delta$ .

Let

$M_n=\left(\begin{array}{cc}{\sigma }_n& 0\\ 0& {\epsilon }_n\end{array}\right)$

Then $\lambda$ satisfy $\det \left(\lambda I-M_n-A\right)=0$ .

$\det \left[\left(\begin{array}{cc}\lambda & 0\\ 0& \lambda \end{array}\right)-\left(\begin{array}{cc}-d_1{\left(n/l\right)}^2& 0\\ 0& -d_2{\left(n/l\right)}^2\end{array}\right)-\left(\begin{array}{cc}1& 0\\ 0& \gamma \end{array}\right)\right]=0.$

Reduce to

$\left(\lambda +d_1{\left(n/l\right)}^2-1\right)\left(\lambda +d_2{\left(n/l\right)}^2-\gamma \right)=0.$ (8)

The eigenvalue of the characteristic Equation (8) is ${\lambda }_1=-d_1{\left(n/l\right)}^2+1$ , ${\lambda }_2=-d_2{\left(n/l\right)}^2+\gamma$ .

Let $n_0=0$ , then ${\lambda }_1=1>0$ . It can be known from [11] that it is an unstable point.

By the same token, the characteristic equation of system (2) is linearized at $E_1$

$\left(\lambda +d_1{\left(n/l\right)}^2+1\right)\left(\lambda +d_2{\left(n/l\right)}^2-\gamma \right)=0.$ (9)

The eigenvalue of the characteristic Equation (9) is ${\lambda }_1=-d_1{\left(n/l\right)}^2-1$ , ${\lambda }_2=-d_2{\left(n/l\right)}^2+\gamma$ .

Let $n_0=0$ , then ${\lambda }_2=\gamma >0$ . So $E_1$ is an unstable point.

4.3. Stability Analysis of Positive Equilibrium Point E₂

If the hypothesis (H₀) is true, let

$U={\left(u\left(x,t\right),v\left(x,t\right)\right)}^{\text{T}},U_{\tau }={\left(0,v\left(x,t-\tau \right)\right)}^{\text{T}}.$

System (2) is linearized at a positive equilibrium point $E_2=\left(u_{*},v_{*}\right)$

$U_t=D_0\Delta U+D_1\Delta U_{\tau }+J_{1}U.$ (10)

where

$D_0=\left(\begin{array}{cc}{d}_1& 0\\ 0& d_2\end{array}\right),D_1=\left(\begin{array}{cc}0& 0\\ 0& d_3{v}_{*}\end{array}\right),J_1=J\left(u_{*},v_{*}\right)=\left(\begin{array}{cc}{a}_{11}^{*}& a_{12}^{*}\\ a_{21}^{*}& a_{22}^{*}\end{array}\right).$

and

$\begin{array}{l}{a}_{11}^{*}=\frac{\left(2\gamma -c\right){\alpha }^{*}-c{\gamma }^2}{c{\alpha }^{*}}=\frac{2\gamma }{c}-\frac{{\gamma }^2}{{\alpha }^{*}}-1,a_{12}^{*}=\frac{\gamma \left(c-\gamma \right){\left(c^2-{\alpha }^{*}\beta \right)}^2}{c^4{\alpha }^{*}},\\ a_{21}^{*}=\frac{c\alpha {\gamma }^2}{c^2-{\alpha }^{*}\beta },a_{22}^{*}=-\gamma +\frac{\beta \gamma {\alpha }^{*}}{c^2}.\end{array}$

The characteristic equation of system (10) is given

$E\left(n,\lambda ,\tau \right)=\det \left(\lambda I-{M^{\prime }}_n-J_1\right)=0.$ (11)

where ${M^{\prime }}_n=-{\mu }_n{D}_0-{\mu }_n{D}_1{\text{e}}^{-\lambda \tau }$ . Simplify (11) can get

$E\left(n,\lambda ,\tau \right)={\lambda }^2+\left(a_n+b_n{\text{e}}^{-\lambda \tau }\right)\lambda +\left(c_n+d_n{\text{e}}^{-\lambda \tau }\right)=0.$ (12)

where

$\begin{array}{l}{a}_n={\left(\frac{n}{l}\right)}^2\left(d_1+d_2\right)-a_{11}^{*}-a_{22}^{*},b_n={\left(\frac{n}{l}\right)}^2{d}_3{v}_{*},\\ c_n={\left(\frac{n}{l}\right)}^2{d}_2\left({\left(\frac{n}{l}\right)}^2{d}_1-a_{11}^{*}\right)-{\left(\frac{n}{l}\right)}^2{d}_1{a}_{22}^{*}+a_{11}^{*}{a}_{22}^{*}-a_{12}^{*}{a}_{21}^{*},\\ d_n={\left(\frac{n}{l}\right)}^2{d}_3{v}_{*}\left({\left(\frac{n}{l}\right)}^2{d}_1-a_{11}^{*}\right).\end{array}$

4.3.1. Stability of Positive Equilibrium Point in System (2) When $\tau =0$

When $\tau =0$ , the characteristic Equation (12) becomes

$E\left(n,\lambda ,0\right)={\lambda }^2+\left(a_n+b_n\right)\lambda +\left(c_n+d_n\right)=0.$ (13)

From (H₀) $0<\theta <c,0<\beta <4\alpha /{\left(\alpha +1\right)}^2$ , We can get

$a_{12}^{*}>0,a_{21}^{*}>0,a_{22}^{*}<0.$

If the conditions are met (H₁) ${\alpha }^{*}\le {\gamma }^2$ , then have $a_{11}^{*}\le 0$ , $a_{11}^{*}{a}_{22}^{*}-a_{12}^{*}{a}_{21}^{*}\ge 0$ , i.e.

$a_n+b_n={\left(\frac{n}{l}\right)}^2\left(d_1+d_2+d_3{v}_{*}\right)-a_{11}^{*}-a_{22}^{*}>0,$

$c_n+d_n=\left({\left(\frac{n}{l}\right)}^2{d}_1-a_{11}^{*}\right)\left({\left(\frac{n}{l}\right)}^2\left(d_2+d_3{v}_{*}\right)\right)-{\left(\frac{n}{l}\right)}^2{d}_1{a}_{22}^{*}+a_{11}^{*}{a}_{22}^{*}-a_{12}^{*}{a}_{21}^{*}>0.$

According to the Routh-Hurwitz discriminant theorem, all eigenroots of (13) have negative real parts, so the following theorem holds.

Theorem 3. Assumes that the condition (H₀), (H₁) is true, then the positive equilibrium point of system (2) is locally asymptotically stable when $\tau =0$ .

Theorem 4. When $\tau =0$ , the positive equilibrium point of system (2) is globally stable if the conditions (H₂) $\alpha <\beta ,\beta \left(1+\alpha \right)<c/\gamma$ are met.

Proof: Let $F_1\left(u,v\right)=1-u-v/\left(1+\alpha u\right)\left(1+\beta v\right)$ , $F_2\left(u,v\right)=\gamma -cv/\left(1+\alpha u\right)\left(1+\beta v\right)$ . Then

$\begin{array}{l}\frac{\partial F_1}{\partial u}=A_{11}=\frac{\alpha v}{{\left(1+\alpha u\right)}^2\left(1+\beta v\right)}-1,\frac{\partial F_1}{\partial v}=A_{12}=-\frac{1}{\left(1+\alpha u\right){\left(1+\beta v\right)}^2},\\ \frac{\partial F_2}{\partial u}=A_{21}=\frac{c\alpha v}{{\left(1+\alpha u\right)}^2\left(1+\beta v\right)},\frac{\partial F_2}{\partial v}=A_{22}=-\frac{c}{\left(1+\alpha u\right){\left(1+\beta v\right)}^2}.\end{array}$

Let

$s=\frac{\gamma \left(1+\alpha \right)}{c-\beta \gamma \left(1+\alpha \right)},$

We can get

$\begin{array}{l}{A}_{11}\le -H_{11}=\frac{\alpha }{\beta }-1,A_{22}\le -H_{22}=-\frac{c}{\left(1+\alpha \right){\left(1+\beta s\right)}^2},\\ \left|A_{12}\right|\le H_{12}=1,\left|A_{21}\right|\le H_{21}=\frac{c\alpha }{\beta }.\end{array}$ (14)

According to Taylor’s theorem, model (2) can be written as

$\{\begin{array}{ll}{u}_t=d_1{u}_{xx}+u\left(A_{21}\left(u-u_{*}\right)+A_{22}\left(v-v_{*}\right)\right),\hfill & x\in \left(0,l\pi \right),t>0,\hfill \\ v_t=d_2{v}_{xx}+d_3{\left(v{v}_x\right)}_x+v\left(A_{21}\left(u-u_{*}\right)+A_{22}\left(v-v_{*}\right)\right),\hfill & x\in \left(0,l\pi \right),t>0.\hfill \end{array}$ (15)

The Lyapunov function is constructed as follows:

$V\left(u,v\right)=c_1{\displaystyle {\int }_0^{l\pi }\left(u-u_{*}-u_{*}\ln \frac{u}{u_{*}}\right)\text{d}x}+c_2{\displaystyle {\int }_0^{l\pi }\left(v-v_{*}-v_{*}\ln \frac{v}{v^{*}}\right)\text{d}x}.$

Where $c_1,c_2$ is the undetermined normal number. The solutions along the above system are

$\frac{\text{d}V}{\text{d}t}=c_1{\displaystyle {\int }_0^{l\pi }\frac{u-u_{*}}{u}{u}_t\text{d}x}+c_2{\displaystyle {\int }_0^{l\pi }\frac{v-v_{*}}{v}{v}_t\text{d}x}:=c_1{I}_1+c_2{I}_2.$ (16)

where

$\begin{array}{c}{I}_1={\displaystyle {\int }_0^{l\pi }\frac{u-u_{*}}{u}{u}_t}\text{d}x={\displaystyle {\int }_0^{l\pi }\frac{u-u_{*}}{u}\left[d_1{u}_{xx}+u\left(A_{11}\left(u-u_{*}\right)+A_{12}\left(v-v_{*}\right)\right)\right]\text{d}x}\\ =-d_1{u}_{*}{\displaystyle {\int }_0^{l\pi }{\left(\frac{u_x}{u}\right)}^2\text{d}x}+{\displaystyle {\int }_0^{l\pi }\left(u-u_{*}\right)\left(A_{11}\left(u-u_{*}\right)+A_{12}\left(v-v_{*}\right)\right)\text{d}x}\\ \le {\displaystyle {\int }_0^{l\pi }\left(A_{11}{\left(u-u_{*}\right)}^2+\left|A_{12}\right|\left|u-u_{*}\right|\left|v-v_{*}\right|\right)\text{d}x}-d_1{u}_{*}{\displaystyle {\int }_0^{l\pi }{\left(\frac{u_x}{u}\right)}^2\text{d}x},\\ I_2={\displaystyle {\int }_0^{l\pi }\frac{v-v_{*}}{v}{v}_t}\text{d}x={\displaystyle {\int }_0^{l\pi }\frac{v-v_{*}}{v}\left[d_2{v}_{xx}+d_3{\left(v{v}_x\right)}_x+v\left[A_{21}\left(u-u_{*}\right)+A_{22}\left(v-v_{*}\right)\right]\right]\text{d}x}\\ =-d_2{v}_{*}{\displaystyle {\int }_0^{l\pi }{\left(\frac{v_x}{v}\right)}^2\text{d}x}-d_3{v}_{*}{\displaystyle {\int }_0^{l\pi }\frac{v_x^2}{v}\text{d}x}+\left(v-v_{*}\right)\left[A_{21}\left(u-u_{*}\right)+A_{22}\left(v-v_{*}\right)\right]\text{d}x\\ \le {\displaystyle {\int }_0^{l\pi }\left(A_{22}{\left(v-v_{*}\right)}^2+\left|A_{21}\right|\left|u-u_{*}\right|\left|v-v_{*}\right|\right)\text{d}x}-d_2{v}_{*}{\displaystyle {\int }_0^{l\pi }{\left(\frac{v_x}{v}\right)}^2\text{d}x}-d_3{v}_{*}{\displaystyle {\int }_0^{l\pi }\frac{v_x^2}{v}\text{d}x}.\end{array}$ (17)

Obtained by (14), (16), (17)

$\begin{array}{c}\frac{\text{d}V}{\text{d}t}=c_1{\displaystyle {\int }_0^{l\pi }\frac{u-u_{*}}{u}{u}_t}\text{d}x+c_2{\displaystyle {\int }_0^{l\pi }\frac{v-v_{*}}{v}{v}_t}\text{d}x\\ \le c_1{\displaystyle {\int }_0^{l\pi }\left(A_{11}{\left(u-u_{*}\right)}^2+\left|A_{12}\right|\left|u-u_{*}\right|\left|v-v_{*}\right|\right)\text{d}x}-c_1{d}_1{u}_{*}{\displaystyle {\int }_0^{l\pi }{\left(\frac{u_x}{u}\right)}^2\text{d}x}\\ +c_2{\displaystyle {\int }_0^{l\pi }\left(A_{22}{\left(v-v_{*}\right)}^2+\left|A_{21}\right|\left|u-u_{*}\right|\left|v-v_{*}\right|\right)\text{d}x}-d_2{v}_{*}{\displaystyle {\int }_0^{l\pi }{\left(\frac{v_x}{v}\right)}^2\text{d}x}-d_3{v}_{*}{\displaystyle {\int }_0^{l\pi }\frac{v_x^2}{v}\text{d}x}\\ \le {\displaystyle {\int }_0^{l\pi }\left(c_1{A}_{11}{\left(u-u_{*}\right)}^2+c_2{A}_{22}{\left(v-v_{*}\right)}^2+c_1\left|A_{12}\right|\left|u-u_{*}\right|\left|v-v_{*}\right|+c_2\left|A_{21}\right|\left|u-u_{*}\right|\left|v-v_{*}\right|\right)\text{d}x}\\ \le {\displaystyle {\int }_0^{l\pi }\left(-c_1{H}_{11}{\left(u-u_{*}\right)}^2-c_2{H}_{22}{\left(v-v_{*}\right)}^2+c_1{H}_{12}\left|u-u_{*}\right|\left|v-v_{*}\right|+c_2{H}_{21}\left|u-u_{*}\right|\left|v-v_{*}\right|\right)\text{d}x.}\end{array}$

Let $V_1=\left|u-u_{*}\right|,V_2=\left|v-v_{*}\right|$ , then

$\frac{\text{d}V}{\text{d}t}=-\frac{1}{2}V\left(CH+H^{\text{T}}C\right)V^{\text{T}},$

where $V=\left(V_1,V_2\right),C=diag\left(c_1,c_2\right)$ .

Let $c_1=1,c_2=H_{12}/H_{21}>0$ , then

$2c_1{H}_{11}>0,2c_2{H}_{22}>0,4c_1{c}_2{H}_{11}{H}_{22}-\left(c_1{H}_{12}+c_2{H}_{21}\right)>0.$

i.e. $CH+H^{\text{T}}C$ is positive definite, and

$\frac{\text{d}V}{\text{d}t}=-\frac{1}{2}V\left(CH+H^{\text{T}}C\right)V^{\text{T}}\le 0,$

is constant, if and only if $\left(u,v\right)=\left(u_{*},v_{*}\right)$ , the equal sign is true. Therefore, when $\tau =0$ , model (2) is globally stable at the positive equilibrium point.

4.3.2. Stability of Positive Equilibrium Point in System (2)

To discuss the stability of the positive equilibrium point of the system with time delay $\tau >0$ , we first consider the case that the characteristic Equation (12) has pure imaginary roots:

Let ${\lambda }_n=i{\omega }_n\left({\omega }_n>0\right)$ is a pure imaginary root of Equation (12), then plug it into (12) can get

$\{\begin{array}{l}{b}_n{\omega }_n\sin \left({\omega }_n\tau \right)+d_n\cos \left({\omega }_n\tau \right)+c_n-{\omega }_n^2=0,\\ b_n{\omega }_n\cos \left({\omega }_n\tau \right)-d_n\sin \left({\omega }_n\tau \right)+a_n{\omega }_n=0.\end{array}$ (18)

By sorting out the Equations (18), we can obtain

${\omega }_n^4+\left(a_n^2-b_n^2-2c_n\right){\omega }_n^2+c_n^2-d_n^2=0.$ (19)

Let $y_n={\omega }_n^2$ , then (19) becomes

$y_n^2+\left(a_n^2-b_n^2-2c_n\right)y_n+c_n^2-d_n^2=0.$ (20)

Solve the Equation (20) to find the root $y_n$

$y_n=\frac{1}{2}\left(b_n^2+2c_n-a_n^2\right)\pm \sqrt{{\left(b_n^2+2c_n-a_n^2\right)}^2-4\left(c_n^2-d_n^2\right)}.$ (21)

Due to

$\begin{array}{c}{a}_n^2-b_n^2-2c_n={\left(n/l\right)}^4\left(d_1^2+d_2^2-d_3^2{v}_{*}^2\right)-2{\left(n/l\right)}^2\left(a_{11}^{*}{d}_1+a_{22}^{*}{d}_2\right)\\ +{\left(a_{11}^{*}+a_{22}^{*}\right)}^2+2a_{12}^{*}{a}_{21}^{*}-2a_{11}^{*}{a}_{22}^{*},\end{array}$

$\begin{array}{c}{c}_n^2-d_n^2=[{\left(n/l\right)}^2\left(d_2+d_3{v}_{*}\right)\left({\left(n/l\right)}^2{d}_1-a_{11}^{*}\right)-{\left(n/l\right)}^2{d}_1{a}_{22}^{*}\\ {}_{}{}^{{}_{}{}^{}}+a_{11}^{*}{a}_{22}^{*}-a_{12}^{*}{a}_{21}^{*}][{\left(n/l\right)}^2\left(d_2-d_3{v}_{*}\right)\left({\left(n/l\right)}^2{d}_1-a_{11}^{*}\right)\\ -{\left(n/l\right)}^2{d}_1{a}_{22}^{*}+a_{11}^{*}{a}_{22}^{*}-a_{12}^{*}{a}_{21}^{*}].\end{array}$

Let $a_{11}^{*}{a}_{22}^{*}-a_{12}^{*}{a}_{21}^{*}=0$ , then

$\alpha =\frac{c^2-\left(2-c\right)\gamma }{2c\gamma -\left(c^2+{\gamma }^2\right)},$

And need to satisfy $\alpha =\frac{c^2-\left(2-c\right)\gamma }{2c\gamma -\left(c^2+{\gamma }^2\right)}>0$ , i.e. $c^2-\left(2-c\right)\gamma <0$ , by $0<\gamma <c$ , there is

$\frac{c^2}{2-c}<\gamma <c,0<c<1.$

Assume

(H₃) $0<d_3<\frac{d_2}{v^{*}},\alpha =\frac{c^2-\left(2-c\right)\gamma }{2c\gamma -\left(c^2+{\gamma }^2\right)},\frac{c^2}{2-c}<\gamma <c,0<c<1.$

From the above discussion, we can see that if the conditions (H₀), (H₁), (H₃) are true, then

$a_n^2-b_n^2-2c_n>0,c_n^2-d_n^2>0$

True, for any $n\in {\mathbb{N}}_0$ . Equation (20) has no positive roots, and its roots are distributed on the negative half plane, i.e.

$\mathrm{Re}{\lambda }_n<0,n\in {\mathbb{N}}_0.$

Then the following theorem holds:

Theorem 5. If the conditions (H₀), (H₁), (H₃) are true, then for all $\tau \ge 0$ , the positive equilibrium point $E_2$ of systems (2) is locally asymptotically stable.

Proof: When the condition (H₁) is true, let $\tau =0$ , all eigenroots of the characteristic Equation (13) have negative real parts; When both the (H₁), (H₃) are true, let $\tau >0$ , all eigenroots of the characteristic Equation (20) have negative real parts. I.e. when condition (H₀), (H₁), (H₃) are true, for all $\tau \ge 0$ , the positive equilibrium point of system (2) is locally asymptotically stable.

Theorem 6. On the premise of $\tau >0$ , let ${\lambda }_0$ is the root of the characteristic Equation (12), and let the set

${\sigma }_n\left(\tau \right)=\left\{\lambda \in ℂ:E\left(n,\tau ,\lambda \right)=0\right\}.$

Then exist $N\in \mathbb{N}$ and ${\left\{{\lambda }_n\right\}}_{n\ge N}\in ℂ$ , s.t. The solution element with a positive real part in ${\sigma }_n\left(\tau \right)$ is completely determined by $E_{\infty }\left(\tau ,\lambda \right)=0$ .

Proof: let $\stackrel{⌢}{E}\left(n,\tau ,\lambda \right)=E\left(n,\tau ,\lambda \right)/{\mu }_n^2$ , i.e.

$\begin{array}{c}\stackrel{⌢}{E}\left(n,\tau ,\lambda \right)=\frac{{\lambda }^2}{{\mu }_n^2}+\left(\frac{d_1+d_2}{{\mu }_n}-\frac{a_{11}^{*}+a_{22}^{*}}{{\mu }_n^2}+\frac{d_3{v}_{*}{\text{e}}^{-\lambda \tau }}{{\mu }_n}\right)\lambda +d_1{d}_2\\ -\frac{d_2{a}_{11}^{*}+d_1{a}_{22}^{*}}{{\mu }_n}-\frac{a_{12}^{*}{a}_{21}^{*}+a_{11}^{*}{a}_{22}^{*}}{{\mu }_n^2}+\left(d_1{d}_3{v}_{*}-\frac{a_{11}^{*}{d}_3{v}_{*}}{{\mu }_n}\right){\text{e}}^{-\lambda \tau }\\ =0.\end{array}$

Let ${\mu }_n^{-1}=\gamma$ , then

$\begin{array}{c}F\left(\gamma ,\lambda \right)=\stackrel{⌢}{E}\left(n,\tau ,\lambda \right)\\ ={\gamma }^2{\lambda }^2+\left[\left(d_1+d_2\right)\gamma -\left(a_{11}^{*}+a_{22}^{*}\right){\gamma }^2\right]\lambda \\ +\left[d_1{d}_2-\left(d_2{a}_{11}^{*}+d_1{a}_{22}^{*}\right)\gamma +\left(a_{12}^{*}{a}_{21}^{*}+a_{11}^{*}{a}_{22}^{*}\right){\gamma }^2\right]\\ +\left(d_3{v}_{*}\gamma \lambda +d_1{d}_3{v}_{*}-a_{11}^{*}{d}_3{v}_{*}\gamma \right){\text{e}}^{-\lambda \tau }\\ =0.\end{array}$

Plug $\gamma =0,\lambda ={\lambda }_0$ into the above equation can get

$F\left(0,{\lambda }_0\right)=d_1{d}_2+d_1{d}_3{v}_{*}{\text{e}}^{-\lambda \tau }=0.$

Suppose that $F:\left(-\delta ,\delta \right)\times ℂ\to ℂ$ , theorem 2.6 in reference [12] shows that there must be a ${\delta }_1\in \left(0,\delta \right)$ and analytic function $\lambda :\left(-{\delta }_1,{\delta }_1\right)\to ℂ$ , s.t. $F\left(\gamma ,\lambda \left(\gamma \right)\right)=0$ and $\lambda \left(0\right)={\lambda }_0$ .

Let

$E_{\infty }\left(\tau ,\lambda \right)=d_2+d_3{v}_{*}{\text{e}}^{-\lambda \tau }=0,$ (22)

Then $\lambda ={\lambda }_0$ . According to Corollary 2.7 of reference [9] , exist $N\in \mathbb{N}$ and ${\left\{{\lambda }_n\right\}}_{n\ge N}\in ℂ$ , makes

$E\left(n,\tau ,{\lambda }_n\right)=0,\underset{n\to \infty }{\lim }{\lambda }_n={\lambda }_0.$

Let ${\sigma }_n\left(\tau \right)=\left\{\lambda \in ℂ:E\left(n,\tau ,\lambda \right)=0\right\}$ , then the solution element with a positive real part on ${\sigma }_n\left(\tau \right)$ is completely determined $E_{\infty }\left(\tau ,\lambda \right)=0$ .

If the conditions (H₄) $d_3>d_2/v_{*}$ are true, theorem 6 shows that when n sufficiently large, the distribution estimation of the eigenroots of the characteristic Equation (12) can be converted to the distribution estimation of the eigenroots of the limit Equation (22).

Theorem 7. If the conditions (H₀), (H₁), (H₄) are true, when $\tau =0$ , the positive equilibrium point $E_2$ of system (2) is locally asymptotically stable, and when $\tau >0$ , the positive equilibrium point $E_2$ of system (2) is unstable.

Proof: According to Theorem 3, when $\tau =0$ , the positive equilibrium point of system (2) is locally asymptotically stable.

When $\tau >0$ , the characteristic root of the limit Equation (22) is assumed to be $\lambda =a\pm ib$ , plug it in the equation

$\{\begin{array}{l}\sin \left(b\tau \right)=0,\\ d_2{\text{e}}^{a\tau }+d_3{v}_{*}\cos \left(b\tau \right)=0.\end{array}$ (23)

Solving the first equation of the system of Equations (23) is obtained

$b=\ln \left(d_3{v}_{*}/d_2\right)/\tau +i\left(\left(2k+1\right)\pi /\tau \right),k\in \mathbb{Z}.$

From $\sin \left(b\tau \right)=0$ , we can get $\cos \left(b\tau \right)=\pm 1$ .

For the solution of a, discuss by case:

When $\cos \left(b\tau \right)=1$ , the second equation of the system (23) is reduced to

$d_2{\text{e}}^{a\tau }+d_3{v}_{*}=0$ (24)

It is known that $d_2,d_3>0$ , so Equation (24) has no solution.

When $\cos \left(b\tau \right)=-1$ , the second equation of the system (23) is reduced to

$d_2{\text{e}}^{a\tau }-d_3{v}_{*}=0.$ (25)

The solution to Equation (25) is

$a=\ln \left(d_3{v}_{*}/d_2\right)/\tau .$

According to the conditions (H₃) $d_3>d_2/v_{*}$ , it is known that $\ln \left(d_3{v}_{*}/d_2\right)>\ln 1>0$ .

I.e. the characteristic roots $\lambda$ of the limit Equation (22) have positive real parts. In other words, the positive real part exists when n with the eigenroots of the characteristic Equation (12) is sufficiently large, then the positive equilibrium point $E_2$ of system (2) is unstable.

5. Conclusion

In this paper, we consider a Bazykin type functional reactive predator prey model based on self memory diffusion with non-negative initial value functions and homogeneous Neumann boundary conditions. Firstly, the uniqueness, boundedness and positivity of the solution of the system with non-negative initial value function at that time delay are proved, which is consistent with the biological meaning. Secondly, the stability at the equilibrium point of the system was studied. We have identified three equilibrium points $E_0,E_1,E_2$ of the system, analyzed the instability of $E_0$ , $E_1$ , and analyzed the stability of the positive equilibrium point $E_2$ under certain conditions $\tau =0$ and $\tau >0$ . The results indicate that it is locally asymptotically stable and globally stable at $\tau =0$ , and unstable at $\tau >0$ , at which point branching may occur. By analyzing the stability of the equilibrium points of the self memory reaction diffusion model under the Bazykin type functional response, it is helpful to further study the various branching problems of this type of model. For example, the Hopf branch. Study the existence of branches, and if there is a Hopf branch, use $\tau$ as a branch parameter to study the direction and stability of the Hopf branch at the positive equilibrium point of the system through the central manifold theorem and normal form theory.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Din, Q. (2013) Dynamics of a Discrete Lotka-Volterra Model. Advances in Difference Equations, 2013, 1-13. https://doi.org/10.1186/1687-1847-2013-95
[2]	Dawes, J.H.P. and Souza, M.O. (2013) A Derivation of Holling’s Type I, II and III Functional Responses in Predator-Prey Systems. Journal of Theoretical Biology, 327, 11-22. https://doi.org/10.1016/j.jtbi.2013.02.017
[3]	Denny, M. (2014) Buzz Holling and the Functional Response. Bulletin of the Ecological Society of America, 95, 200-203. https://doi.org/10.1890/0012-9623-95.3.200
[4]	Chang, X. and Wei, J. (2013) Stability and Hopf Bifurcation in a Diffusive Predator-Prey System Incorporating a Prey Refuge. Mathematical Biosciences and Engineering, 10, 979-996. https://doi.org/10.3934/mbe.2013.10.979
[5]	Li, X., Jiang, W. and Shi, J. (2013) Hopf Bifurcation and Turing Instability in the Reaction-Diffusion Holling-Tanner Predator-Prey Model. IMA Journal of Applied Mathematics, 78, 287-306. https://doi.org/10.1093/imamat/hxr050
[6]	Shi, J., Wang, C., Wang, H., et al. (2020) Diffusive Spatial Movement with Memory. Journal of Dynamics and Differential Equations, 32, 979-1002. https://doi.org/10.1007/s10884-019-09757-y
[7]	Shi, J., Wang, C. and Wang, H. (2021) Spatial Movement with Diffusion and Memory-Based Self-Diffusion and Cross-Diffusion. Journal of Differential Equations, 305, 242-269. https://doi.org/10.1016/j.jde.2021.10.021
[8]	Song, Y., Peng, Y. and Zhang, T. (2021) The Spatially Inhomogeneous Hopf Bifurcation Induced by Memory Delay in a Memory-Based Diffusion System. Journal of Differential Equations, 300, 597-624. https://doi.org/10.1016/j.jde.2021.08.010
[9]	Li, S., Li, Z. and Dai, B. (2022) Stability and Hopf Bifurcation in a Prey-Predator Model with Memory-Based Diffusion. Discrete and Continuous Dynamical Systems-B, 27, 6885-6906. https://doi.org/10.3934/dcdsb.2022025
[10]	Dung, L. (1998) Global L∞ Estimates for a Class of Reaction-Diffusion Systems. Journal of Mathematical Analysis and Applications, 217, 72-94. https://doi.org/10.1006/jmaa.1997.5702
[11]	Wu, J. (1996) Theory and Applications of Partial Functional Differential Equations. Springer Science & Business Media, Berlin. https://doi.org/10.1007/978-1-4612-4050-1
[12]	Shi, J., Wang, C. and Wang, H. (2021) Spatial Movement with Diffusion and Memory-Based Self-Diffusion and Cross-Diffusion. Journal of Differential Equations, 305, 242-269. https://doi.org/10.1016/j.jde.2021.10.021