Bifurcationing Analysis of Predator-Prey Diffusive System Based on Bazykin Functional Response

Mingyang Zhao, Fuqin Sun

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

DOI: 10.4236/jamp.2022.1012254                

Abstract

A predator-prey diffusion system with a Bazykin functional response is studied. The existence of equilibrium points, the stability of normal number equilibrium points and the existence of Hopf bifurcationes are investigated for the proposed system, the existence of positive solutions in the system is discussed under Neumann boundary conditions, and the stability of constant equilibrium points is focused on under the condition of Hurwitz criterion. The results show that there exist positive equilibrium points in the system under Neumann boundary conditions, and the normal number equilibrium points are stable when specific conditions are satisfied, and the bifurcation points of Hopf bifurcationes and their orders are given.

Keywords

Bazykin Functional Response, Diffusion System, Existence, Stability, Hopf Bifurcation

1. Introduction

In nature, there is an interdependent and mutually constraining way of survival between different populations: population A thrives on abundant natural resources, while population B feeds on population A, such as fish and sharks, American rabbits and bobcats, larch and aphids, etc. In ecology, population A is called predator and population B is called predator, and together they form a predator-prey system. The predator-prey system has a long and distinguished history in applied mathematics as a classical mathematical system for the study of predator-prey systems, which is originated from the exploration of chemical reactions [1] and biological interactions [2] and has been widely studied for its rich kinetic behavior [3] [4] [5] [6].

At present, Holling-type and Beddington-type functional response functions have been widely used in the study of predator-prey systems, but Bazykin functional response functions have been less studied. Bazykin functional response can describe the unstable force of predator saturation and the stable force of prey competition, which is more practical for understanding the role between populations. Moreover, in the real world, prey and predators are always in motion, and to accurately system the dynamic characteristics of prey and predators, the spatial spread of populations needs to be considered in the system [7] [8] [9]. Therefore, a predator-prey diffusive system with Bazykin functional response was constructed in this paper.

2. System Composition

We considered a predator-prey diffusion system with a Bazykin functional response, as shown in Equation (2.1).

$\{\begin{array}{l}\frac{\text{d}\tilde{u}}{\text{d}t}=\tilde{u}\left(\tilde{a}-\tilde{u}-\tilde{b}\tilde{v}\tilde{f}\left(\tilde{u},\tilde{v}\right)\right)\hfill \\ \frac{\text{d}\tilde{v}}{\text{d}t}=\tilde{v}\left(\tilde{c}-\tilde{d}\tilde{v}\tilde{f}\left(\tilde{u},\tilde{v}\right)\right)\hfill \end{array}$ (2.1)

In the real world, prey and predators are always in motion, and this phenomenon can be simulated using self-diffusion. Assume that the two populations under consideration remain in motion and each population follows a path, the length of which is denoted by x. Considering the above assumptions, the system (2.1) can be written in the following form

$\{\begin{array}{l}\begin{array}{cc}\frac{\text{d}\tilde{u}}{\text{d}t}-{\tilde{d}}_1{u}_{xx}=\tilde{u}\left(\tilde{a}-\tilde{u}-\tilde{b}\tilde{v}\tilde{f}\left(\tilde{u},\tilde{v}\right)\right),& x\in \left(0,l\pi \right),t>0\end{array}\hfill \\ \begin{array}{cc}\frac{\text{d}\tilde{v}}{\text{d}t}-{\tilde{d}}_2{v}_{xx}=\tilde{v}\left(\tilde{c}-\tilde{d}\tilde{v}\tilde{f}\left(\tilde{u},\tilde{v}\right)\right),& x\in \left(0,l\pi \right),t>0\end{array}\hfill \end{array}$ (2.2)

where x is the location of the prey or the predator at the time t, and $l\pi$ is the domain size, $u\left(t\right)$ is the number of prey populations at time t, $v\left(t\right)$ is the number of predator populations at time t, $d_1,d_2$ is the diffusion rate of prey and predators, respectively, $\tilde{a},\tilde{b},\tilde{c},\tilde{d}$ is a positive parameter, $f\left(u,v\right)=1/\left(1+\tilde{\alpha }u\right)\left(1+\tilde{\beta }v\right)$ is a Bazykin functional response function, $\tilde{\alpha },\tilde{\beta }$ is two positive constants, and $\tilde{\alpha }>\tilde{\beta }$, it is used to describe the destabilizing force of predator saturation and the stabilizing force of prey contention.

In order to study the pattern formation of (2.2), we consider the system (2.3) with the homogeneous Neumann boundary conditions.

Where $u{=}_1\tilde{u}$, $v=\tilde{d}\tilde{v}/\tilde{c}$, $a=\tilde{a}$, $b=\tilde{b}\tilde{c}/\tilde{d}$, $\lambda =\tilde{c}$, $d_1={\tilde{d}}_1$, $d_2={\tilde{d}}_2$, $\alpha =\tilde{\alpha }$, $\beta =\tilde{\beta }\tilde{c}/\tilde{d}$, $f\left(u,v\right)=1/\left(1+\alpha u\right)\left(1+\beta v\right)$.

The problem (2.3) with homogeneous Neumann boundary conditions and non-negative initial values is shown below:

$\{\begin{array}{l}\begin{array}{cc}\frac{\text{d}u}{\text{d}t}-d_1{u}_{xx}=u\left(a-u-bvf\left(u,v\right)\right),& x\in \left(0,l\pi \right),t>0\end{array}\hfill \\ \begin{array}{cc}\frac{\text{d}v}{\text{d}t}-d_2{v}_{xx}=\lambda v\left(1-vf\left(u,v\right)\right),& x\in \left(0,l\pi \right),t>0\end{array}\hfill \\ \begin{array}{cc}{u}_x\left(0,t\right)=u_x\left(l\pi ,t\right)=v_x\left(0,t\right)=v_x\left(l\pi ,t\right)=0,& \forall t\ge 0\end{array}\hfill \\ \begin{array}{cc}\begin{array}{cc}u\left(x,0\right)=\phi \left(x\right)\ge 0,& v\left(x,0\right)=\psi \left(x\right)\ge 0,\end{array}& x\in \left(0,l\pi \right)\end{array}\hfill \end{array}$ (2.3)

where x denotes the location of the predator or prey at time t, $l\pi$ denotes the extent of the region, and the Neumann boundary condition indicates that the prey and predator move between 0 and $l\pi$. All of the above parameters are non-negative.

3. Preliminary Results

In this section, the local stability of the positive equilibrium $E_3$ has been studied in the presence of diffusion and also Hopf bifurcation has been analyzed (calculating the Hopf bifurcation points and their order).

3.1. Existence of Equilibrium Points

The equilibrium point of system (2.3) is obtained by solving for $\text{d}u/\text{d}t=0$, $\text{d}v/\text{d}t=0$. Assuming that $a,b,\alpha ,\beta$ satisfies $0<\alpha \beta \left(a-b\right)<1-\beta$, we obtain that system (2.3) has one ordinary constant equilibrium point $E_0\left(0,0\right)$, two semi-ordinary constant equilibrium points $E_1\left(a,0\right),E_2\left(0,1/\left(1-\beta \right)\right)$, and a unique normal number equilibrium point $E_3\left(u^{*},v^{*}\right)$. $u^{*},v^{*}$ can be solved from the system of Equations (3.1) to yield:

$\{\begin{array}{l}{u}^{*}\left(a-u^{*}-b{v}^{*}f\left(u^{*},v^{*}\right)\right)=0\hfill \\ \lambda v^{*}\left(1-v^{*}f\left(u^{*},v^{*}\right)\right)=0\hfill \end{array}$ (3.1)

From the system of Equations (3.1), it is obtained that:

$u^{\ast }=a-b$, $v^{*}=\frac{1+\alpha \left(a-b\right)}{1-\beta \left[1+\alpha \left(a-b\right)\right]}$ (3.2)

3.2. Stability of the Equilibrium Point E*

The Jacobi matrix of the system (2.3) at the positive equilibrium point $\left(u^{*},v^{*}\right)$ is:

$J=\left(\begin{array}{cc}A& B\\ C& D\end{array}\right)$ (3.3)

$A=a-2u^{*}-\frac{b{v}^{*}}{{\left(1+\alpha u^{*}\right)}^2\left(1+\beta v^{*}\right)}$ (3.4)

$B=\frac{-b{u}^{*}}{\left(1+\alpha u^{*}\right){\left(1+\beta v^{*}\right)}^2}$ (3.5)

$C=\frac{\lambda v^{*}{}^2\alpha }{{\left(1+\alpha u^{*}\right)}^2\left(1+\beta v^{*}\right)}$ (3.6)

$D=\lambda -\frac{\lambda \left(2v^{*}+\beta v^{*}{}^2\right)}{\left(1+\alpha u^{*}\right){\left(1+\beta v^{*}\right)}^2}$ (3.7)

Consider the following problem:

$\{\begin{array}{l}\begin{array}{cc}-{\varphi }^{″}=k\varphi & x\in \left(0,l\pi \right)\end{array}\\ {\varphi }^{\prime }\left(0\right)={\varphi }^{\prime }\left(l\pi \right)=0\end{array}$ (3.8)

the eigenvalue of the problem (3.3) are $k_n={\left(n/l\right)}^2,n=0,1,2,\cdots$ and the corresponding eigenfunction are

${\varphi }_n\left(x\right)=\{\begin{array}{cc}\sqrt{\frac{1}{l\pi }},& n=0\\ \sqrt{\frac{2}{l\pi }},& n=1,2,3,\cdots \end{array}$ (3.9)

${\left\{{\varphi }_n\left(x\right)\right\}}_1^{\infty }$ describes the orthogonal basis of $L^2\left(0,l\pi \right)$.

$D=\text{diag}\left(d_1,d_2\right)$ and its corresponding real-valued Sobolev space with Neumann boundary conditions can be defined as

$\chi =\left\{\begin{array}{cc}U={\left(u,v\right)}^{\text{T}}\in W^{2,2}|u_x=v_x=0& x=0,l\pi \end{array}\right\}$ (3.10)

$\Omega =L^2\left(0,l\pi \right)\times L^2\left(0,l\pi \right)$ is the Hibert space whose inner product is

$\langle U_1,U_2\rangle ={\displaystyle {\int }_0^{l\pi }\left(u_1{u}_2+v_1{v}_2\right)\text{d}x}$ (3.11)

The Hilbertian parametrization of the associated $\chi$ is denoted as ${\Vert \cdot \Vert }_{2.2}$. Define the following mapping $F:\left(0,\infty \right)\times \chi \to \Omega$

$F\left(\lambda ,U\right)=\left(\begin{array}{c}{d}_1{u}_{xx}+u\left(a-u-bvf\left(u,v\right)\right)\\ d_2{v}_{xx}+\lambda v\left(1-vf\left(u,v\right)\right)\end{array}\right)$ (3.12)

where $U={\left(u,v\right)}^{\text{T}}$. Then for any ${\left(u,v\right)}^{\text{T}}\in \chi$, $U={\left(u,v\right)}^{\text{T}}$ is the solution of (2.3) $\iff F\left(\lambda ,U\right)=0$.

The Frechet derivative of $F\left(\lambda ,U\right)$ with respect to U is:

$\begin{array}{c}L\left(\lambda \right)=\text{diag}\left(d_1\Delta ,d_2\Delta \right)+J\left(u^{*},v^{*}\right)\\ =\left(\begin{array}{cc}2b-a-\frac{b}{c}-d_1{\left(\frac{n}{l}\right)}^2& -\frac{b\left(c-1\right){\left(1-\beta c\right)}^2}{\alpha c}\\ \frac{\lambda \alpha }{1-\beta c}& -\lambda \left(1-\beta c\right)-d_2{\left(\frac{n}{l}\right)}^2\end{array}\right)\end{array}$ (3.13)

The characteristic equation of $L\left(\lambda \right)$ is $L\left(\lambda \right)\left(u,v\right)=\eta \left(u,v\right)$, and $\eta$ is the eigenvalue, we derive:

$\{\begin{array}{ll}\left[2b-a-\frac{b}{c}-d_1{\left(\frac{n}{l}\right)}^2\right]u-\left[\frac{b\left(c-1\right){\left(1-\beta c\right)}^2}{\alpha c}\right]v=\eta u\hfill & x\in \left(0,l\pi \right)\hfill \\ \left[\frac{\lambda \alpha }{1-\beta c}\right]u-\left[\lambda \left(1-\beta c\right)+d_2{\left(\frac{n}{l}\right)}^2\right]v=\eta v\hfill & x\in \left(0,l\pi \right)\hfill \\ u_x=v_x=0\hfill & x=0,l\pi \hfill \end{array}$ (3.14)

Letting $u={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}{a}_n{\varphi }_n},v={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}{b}_n{\varphi }_n}$, the characteristic equation becomes of the following form.

${\displaystyle \underset{n=0}{\overset{\infty }{\sum }}\left(\text{diag}\left(-d_1{k}_n,-d_2{k}_n\right)+J\left(E^{*}\right)\right)\left(\begin{array}{c}{a}_n\\ b_n\end{array}\right){\varphi }_n}=0$ (3.15)

Letting $\left|\text{diag}\left(-d_1{k}_n,-d_2{k}_n\right)+J\left(E^{*}\right)\right|=0,n=0,1,2,\cdots$, we obtain

${\eta }^2-T_n\left(\lambda \right)\eta +D_n\left(\lambda \right)=0$ (3.16)

$T_n\left(\lambda \right)=2b-a-\frac{b}{c}-\lambda \left(1-\beta c\right)-\left(d_1+d_2\right){\left(\frac{n}{l}\right)}^2$ (3.17)

$\begin{array}{c}{D}_n\left(\lambda \right)=\lambda \left(1-\beta c\right)\left(a-b\right)-\left(2b-a-\frac{b}{c}\right)d_2{\left(\frac{n}{l}\right)}^2\\ +\lambda \left(1-\beta c\right)d_1{\left(\frac{n}{l}\right)}^2+\left(d_1+d_2\right){\left(\frac{n}{l}\right)}^2\end{array}$ (3.18)

The linearized system at equilibrium point $E_3\left(u^{*},v^{*}\right)$ is $U_t=D\Delta U+J\left(u^{*},v^{*}\right)U$. Let $J_2$ be the Jacobian matrix of the system (2.3) at equilibrium point $E_3$. Then

$J_2=-{\left(n/l\right)}^{2}D+J\left(u^{*},v^{*}\right)=\left(\begin{array}{cc}2b-a-\frac{b}{c}-d_1{\left(\frac{n}{l}\right)}^2& -\frac{b\left(c-1\right){\left(1-\beta c\right)}^2}{\alpha c}\\ \frac{\lambda \alpha }{1-\beta c}& -\lambda \left(1-\beta c\right)-d_2{\left(\frac{n}{l}\right)}^2\end{array}\right)$ (3.19)

Letting $c=1+\alpha \left(a-b\right)$, from $0<\alpha \beta \left(a-b\right)<1-\beta$ we get $1<c<1/\beta$.

The eigenvalues of the matrix (3.19) are the solutions of the characteristic equations given by Equation (3.16) (3.17) (3.18).

Theorem 1. When $T_n\left(\lambda \right)<0$, $E_3\left(u^{*},v^{*}\right)$ is locally asymptotically stable, and when $T_n\left(\lambda \right)<0$, $E_3\left(u^{*},v^{*}\right)$ is unstable.

3.3. Existence of Hopf Bifurcations

A Hopf bifurcation occurs when and only when $T_n\left(\lambda \right)=0,D_n\left(\lambda \right)>0$, clearly $D_0\left(\lambda \right)>0$ and $\underset{n\to +\infty }{\lim }{D}_n\left(\lambda \right)>0$. The bifurcation parameter $\lambda$ must be a solution of $\lambda$ in Equation (3.20)

$2b-a-\frac{b}{c}-\lambda \left(1-\beta c\right)-\left(d_1+d_2\right){\left(\frac{n}{l}\right)}^2=0$ (3.20)

The equivalent of

$\lambda =\lambda \left(n\right)$ (3.21)

$\lambda \left(n\right)=\frac{2b-a-\frac{b}{c}-\left(d_1+d_2\right){\left(\frac{n}{l}\right)}^2}{1-\beta c}$ (3.22)

Theorem 2.

Letting

$N_1=\left\{n\in {\rm N}/2b-a-b/c-\lambda \left(1-\beta c\right)-\left(d_1+d_2\right){\left(n/l\right)}^2\right\}$ (3.23)

The system (2.3) occurs Hopf bifurcationing at $\lambda =\lambda \left(n\right)$ and $n>N_1$ (where $\lambda \left(n\right)$ see (3.22)), and $\lambda \left(n\right)$ has the following estimate:

$\lambda \left(0\right)>\lambda \left(1\right)>\cdots >\lambda \left(n\right)>\lambda \left(n+1\right)>\cdots >\lambda \left(N_1\right)$ (3.24)

Proof: Hopf bifurcationing occurs when and only when $T_0\left(\lambda \right)=0$, equivalently:

$T_0\left(\lambda \right)=\left(d_1+d_2\right){\left(\frac{n}{l}\right)}^2$ (3.25)

${T^{\prime }}_0\left(\lambda \right)=-\left(1-\beta c\right)<0$ (3.26)

$\underset{\lambda \to +\infty }{\lim }{T}_0\left(\lambda \right)=2b-a-\frac{b}{c}-\lambda \left(1-\beta c\right)<0$ (3.27)

According to (3.26) and (3.27), $T_0\left(\lambda \right)$ is strictly decreasing with respect to $\lambda$. Equation (3.25) has a solution when and only when the positive integers n satisfy:

$\left(d_1+d_2\right){\left(\frac{n}{l}\right)}^2>2b-a-\frac{b}{c}-\lambda \left(1-\beta c\right)$ (3.28)

Equivalent to $n>N_1$.

Equation $\left(d_1+d_2\right){\left(n/l\right)}^2$ is strictly decreasing with respect to $\lambda$, which leads to the estimate in (3.23).

Now we take $\xi \left(\lambda \right)=\eta \left(\lambda \right)\pm i\omega \left(\lambda \right)$ as the solution of the characteristic equation, where:

$\eta \left(\lambda \left(n\right)\right)=0$, $\omega \left(\lambda \left(n\right)\right)=\sqrt{D\left(\lambda \left(n\right)\right)}$ (3.29)

${\eta }^{\prime }\left(\lambda \left(n\right)\right)=\frac{-\left(1-\beta c\right)}{2}$ (3.30)

Under condition (3.30), the bifurcation point and its order are given by the following theorem.

Theorem 3. If there exists $N^{*}\ge N_1$, a critical value $j_0,\cdots ,j_{N^{*}}$ such that $j_0=0>j_1>\cdots >j_{N^{*}}>N_1$, there are the following estimates:

$\lambda \left(0\right)>\lambda \left(1\right)>\cdots >\lambda \left(n\right)>\lambda \left(n+1\right)>\cdots >\lambda \left(j_{N^{*}}\right)$ (3.31)

4. Conclusion

The article considers a predator-prey diffusion system with Bazykin functional response under chi-square Neumann boundary conditions. Firstly, the existence of equilibrium points of the system is proved, and four equilibrium points $E_0,E_1,E_2,E_3$ of the system are obtained. Secondly, the local stability of $E_3\left(u^{*},v^{*}\right)$ and the existence of Hopf bifurcations under specific conditions are proved by analyzing the characteristic equations of the equilibrium point $E_3\left(u^{*},v^{*}\right)$.

Conflicts of Interest

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

References

[1]	Lotka, A.J. (1925) Elements of Physical Biology. Nature, 116, 461. https://doi.org/10.1038/116461b0
[2]	Volterra, V. (1926) Fluctuations in the Abundance of a Species Considered Mathematically. Nature, 118, 558-560. https://doi.org/10.1038/118558a0
[3]	Bazykin, A.D. (1998) Nonlinear Dynamics of Interacting Populations. World Scientific Series on Nonlinear Science Series A: Volume 11, World Scientific, Singapore, p. 216. https://doi.org/10.1142/2284
[4]	Murdoch, W., Briggs C.J. and Nisbet, R.M. (2003) Consumer-Resource Dynamics. Vol. 36, Princeton University Press, Princeton, NJ.
[5]	Turchin, P. (2003) Complex Population Dynamics: A Theoretical/Empirical Synthesis. Princeton University Press, Princeton, NJ.
[6]	Wang, S.T., Yu, H.G., Dai, C.J. and Zhao, M. (2020) The Dynamical Behavior of a Certain Predator-Prey System with Holling Type II Functional Response. Journal of Applied Mathematics and Physics, 8, 527-547. https://doi.org/10.4236/jamp.2020.83042
[7]	Song, Y. and Tang, X. (2017) Stability, Steady-State Bifurcation, and Turing Patterns in Predator-Prey System with Herd Behavior and Prey-Taxis. Studies in Applied Mathematics, 139, 371-404. https://doi.org/10.1111/sapm.12165
[8]	Yang, R., Liu, M. and Zhangt, C. (2018) A Diffusive Predator-Prey System with Additional Food and Intra-Specific Competition among Predators. International Journal of Biomathematics, 11, Article ID: 1850060. https://doi.org/10.1142/S1793524518500602
[9]	Yang, W. (2018) Analysis on Existence of Bifurcation Solutions for a Predator-Prey System with Herd Behavior. Applied Mathematical Modelling, 53, 433-446. https://doi.org/10.1016/j.apm.2017.09.020