Hopf Bifurcation for a Predator-Prey System with θ Logistic Growth

Suqi Ma^1*, Meisheng Li²

¹Department of Mathematics, China Agricultural University, Beijing, China.
²Department of Mathematics, Beijing Beihang University, Beijing, China.

DOI: 10.4236/ijmnta.2025.144006                

Abstract

In predator-prey population model, the density capacity restricts the growth rate of prey population. Prey species is assumed with $\theta$ logistic growth rate, which is $\theta$ exponent dependent. The generalized Hopf point is found, hence the Bautin bifurcation with limit cycle hysteresis phenomena is observed, which brings forth the limit point cycle appearing on Bautin bifurcation line. The homoclinc bifurcation line starts from Bogdanov-Takens point and tangents to saddle-node bifurcation line. The homoclinic orbits and near dynamics are simulated, which manifests the complex dynamics of system.

Keywords

Homoclinc Solution, Bogdanov-Takens Bifurcation, Generalized Hopf Bifurcation

1. Introduction

The investigation of models for population dynamics of predator-prey system has long been continued topics after the famous Lotka-Voterra equations. The focus of population themes is ubiquitous both in ecology and mathematical ecology. With the enthusiasm in the study trends of highly nonlinearity of population dynamics, Holling’s functional response function provides the best understanding of relationship between predator and prey population. The usual Holling function is named as several types and they are often described by exciting feedback regulation function [1]-[4]. With the endeavor to describe prey refuge rate, people replace Holling function by introducing a new refuge rate in saturated prey population density [5]-[8] and Allie effects [9] [10].

Motivated by the work of Chowdhury et al. [10], the species with $\theta$ logistic growth rate is governed by the listed formula

$x^{\prime }=x\left(1-{\left(\frac{x}{K}\right)}^{\theta }\right)$ (1)

With Holling-I type functional response, predator-prey population is constructed by the following ODEs

$\begin{array}{l}\dot{x}=x\left(1-{\left(\frac{x}{K}\right)}^{\theta }\right)-\frac{\beta xy}{1+x}-bx\\ \dot{y}=\frac{xy}{1+x}-ay-\delta y^2\end{array}$ (2)

wherein $K$ represents the carrying capacity, $a$ and $b$ are the death rate of predator population and prey respectively, and $\delta$ is the perturbation parameter. Noticed in Holling function [3] [8] [10], the coefficient $\beta$ is assumed as positive constant, and we discuss both cases of $\beta >1$ and $\beta <1$ . Replace $y$ with $sy$ in Equation (2), both cases can be satisfied.

The generalized Hopf bifurcation (GH) happens in the case of the first Lyapunov exponent becoming zero as varying free parameter. As usual, the limit point cycle is observed near GH point since the stable limit cycle collides with the unstable limit cycle. Hence, the near dynamics of GH point is classified topologically and geometrically [11] [12]. The normal form computation is expanded to the fifth order and GH point is calculated along Hopf line. In view of three-dimensional free parameter space, GH point arrays are observed and limit point cycle bifurcation lines are plotted.

Bogdanov-Takens (BT) point is found as two eigen roots of the characteristic equation becoming zero. Using Matcont software, homoclinic line starts from BT point and is calculated and continued as varying free parameters [13] [14]. Matcont software is the handtools in dynamical bifurcation analysis, and has quick speed under the Windows circumstance. The dynamics bifurcation diagram is drawn as using Matcont program. The homoclinic solution bifurcates from BT point by Matcont [15] [16]. The saddle-node point line bifurcates from BT point, then cusp point is observed, and another new branch of saddle-node line expands from cusp point again. The point $P$ lying on the bifurcating branch of cusp point collides with the homoclinic curve starting from BT point. The second case of BT point is the heteroclinic orbit bifurcation. System has three equilibrium solutions. As heteroclinic orbits from unstable focus to stable node are built, the boundary is calculated by Matcont, which is the stable manifold of a saddle. Hence after, the heteroclinic bifurcation line starts from BT point and is continued as varying free parameter [13] [14]. As the third case of BT bifurcation, the homoclinic line tangents to the saddle node line branch expanding from cusp point, then two lines coalesce.

The whole paper is organized as the listed sections. In Section 2, the generalized bifurcation is analyzed and the near dynamics of GH point is analyzed. Based on the center manifold theory, the normal form is computed and the second Lyapunov exponent is calculated. In Section 3, BT points are calculated and three kinds of homclinic or heteroclinic bifurcation are simulated and analyzed. In Section 4, the dynamical bifurcation diagrams are analyzed by Matcont software. Finally, a discussion is given in Section 5.

2. The Generalized Hopf Bifurcation

System (2) has equilibrium solution $E^{\text{*}}=\left(x^{\text{*}},y^{\text{*}}\right)$ that satisfies

$y^{\text{*}}=-\frac{a{x}^{\text{*}}+a-x^{\text{*}}}{\delta \left(1+x^{\text{*}}\right)}$ (3)

and

$-x^{\text{*}}\left(-1+{\left(\frac{x^{\text{*}}}{K}\right)}^{\theta }\right)+\frac{\left(\left(a-1\right)x^{\text{*}}+a\right)\beta x^{\text{*}}}{\delta {\left(1+x^{\text{*}}\right)}^2}-b{x}^{\text{*}}=0$ (4)

The linearized system of system (2) is written as

$\left|\begin{array}{cc}-{\left(\frac{x^{\text{*}}}{K}\right)}^{\theta }\theta +\frac{\beta x^{\text{*}}{y}^{\text{*}}}{{\left(1+x^{\text{*}}\right)}^2}& -\frac{\beta x^{\text{*}}}{1+x^{\text{*}}}\\ \frac{y^{\text{*}}}{1+x^{\text{*}}}-\frac{x^{\text{*}}{y}^{\text{*}}}{{\left(1+x^{\text{*}}\right)}^2}& -\delta y^{\text{*}}\end{array}\right|$ (5)

The corresponding characteristic equation is written as

${\lambda }^2-tr\left(A\right)\lambda +det\left(A\right)=0$ (6)

with

$\begin{array}{l}\begin{array}{l}tr\left(A\right)=\delta +{\left(\frac{x^{\text{*}}}{K}\right)}^{\theta }\theta -\frac{\beta x^{\text{*}}{y}^{\text{*}}}{{\left(1+x^{\text{*}}\right)}^2}\\ det\left(A\right)=\delta \left(-{\left(\frac{x^{\text{*}}}{K}\right)}^{\theta }\theta +\frac{\beta x^{\text{*}}{y}^{\text{*}}}{{\left(1+x^{\text{*}}\right)}^2}\right)+\frac{\beta x^{\text{*}}{y}^{\text{*}}}{{\left(1+x^{\text{*}}\right)}^2}\end{array}\hfill \end{array}$ (7)

The sufficient condition for fold bifurcation is given as $det\left(A\right)=0$ .

The sufficient condition for Hopf bifurcation is listed as $tr\left(A\right)=0$ , $det\left(A\right)>0$ . The characteristic Equation (6) has imaginary roots $\lambda =\pm i\omega$ with $\omega =\sqrt{det\left(A\right)}$ . Suppose Hopf bifurctaion happens at $P=\left(K^{\text{*}},a^{\text{*}}\right)$ , assume that the characteristic equation has imaginary roots $\pm i\omega$ , and the eigenvector is $q$ given that

$Lq=i\omega q,L^{\text{*}}p=-i\omega p$ (8)

and $\langle p,q\rangle =1$ . For computing the normal form at the generalized Hopf point, we write the trunction Taylor expansion to the fifth order of multilinear form

$x^{\prime }=L\left(\alpha \right)x+B\left(x,x\right)+C\left(x,x,x\right)+D\left(x,x,x,x\right)+E\left(x,x,x,x,x\right)+o\left(x^5\right)$ (9)

Substitute $x=zq+\overline{z}\overline{q}$ into the above equation to get

$\begin{array}{c}{z}^{\prime }=\lambda \left(\alpha \right)z+\langle p\left(\alpha \right),B\left(q,q\right)\rangle z^2+2\langle p\left(\alpha \right),B\left(q,\overline{q}\right)\rangle z\overline{z}+\langle p\left(\alpha \right),B\left(\overline{q},\overline{q}\right)\rangle {\overline{z}}^2\\ +\langle p\left(\alpha \right),C\left(q,q,q\right)\rangle z^3+3\langle p\left(\alpha \right),C\left(q,\overline{q},q\right)\rangle z^2\overline{z}+3\langle p\left(\alpha \right),C\left(q,\overline{q},\overline{q}\right)\rangle z{\overline{z}}^2\end{array}$

$\begin{array}{c}+\langle p\left(\alpha \right),C\left(\overline{q},\overline{q},\overline{q}\right)\rangle {\overline{z}}^3+\langle p\left(\alpha \right),D\left(q,q,q,q\right)\rangle z^4+4\langle p\left(\alpha \right),D\left(q,\overline{q},q,q\right)\rangle z^3\overline{z}\\ +6\langle p\left(\alpha \right),D\left(q,\overline{q},\overline{q},q\right)\rangle z^2{\overline{z}}^2+4\langle p\left(\alpha \right),D\left(q,\overline{q},\overline{q},\overline{q}\right)\rangle z{\overline{z}}^3\\ +\langle p\left(\alpha \right),D\left(\overline{q},\overline{q},\overline{q},\overline{q}\right)\rangle {\overline{z}}^4+10\langle p\left(\alpha \right),E\left(q,q,\overline{q},\overline{q},q\right)\rangle z^3{\overline{z}}^2\end{array}$(10)

set

$\begin{array}{l}z=w+h_{20}{w}^2+h_{11}w\overline{w}+h_{02}{\overline{w}}^2+h_{30}{w}^3+h_{21}{w}^2\overline{w}+h_{12}w{\overline{w}}^2+h_{03}{\overline{w}}^3+\cdots \\ z^{\prime }=w^{\prime }+2h_{20}w{w}^{\prime }+h_{11}{w}^{\prime }\overline{w}+h_{11}w{\overline{w}}^{\prime }+2h_{02}\overline{w}{\overline{w}}^{\prime }+3h_{30}{w}^2{w}^{\prime }+2h_{21}w\overline{w}{w}^{\prime }\\ +h_{21}{w}^2{\overline{w}}^{\prime }+h_{12}{w}^{\prime }\overline{w}+2h_{12}w\overline{w}{\overline{w}}^{\prime }+3h_{03}{\overline{w}}^2{\overline{w}}^{\prime }+\cdots \end{array}$ (11)

Equation (16) is equivalent to its normal form style, which is universal. Therefore, the universal form at GH point is

$z^{\prime }=\lambda z+C\left(\lambda \right)z^2\overline{z}+F\left(\lambda \right)z^3{\overline{z}}^2$ (12)

Compared with Equation (11), one gets that

$\begin{array}{l}C\left(\lambda \right)=\frac{2{\overline{g}}_{02}{g}_{02}\lambda \overline{\lambda }-{\overline{g}}_{11}{g}_{11}{\lambda }^2+2{\overline{g}}_{11}{g}_{11}\lambda \overline{\lambda }-g_{11}{g}_{20}{\lambda }^2}{\lambda \left(\lambda -2\overline{\lambda }\right)\overline{\lambda }}\\ +\frac{4g_{11}{g}_{20}{\overline{\lambda }}^2+g_{21}{\lambda }^2\overline{\lambda }-2g_{21}\lambda {\overline{\lambda }}^2}{\lambda \left(\lambda -2\overline{\lambda }\right)\overline{\lambda }}\end{array}$

Do polar axis transformation, one derives from Equation (12)

${\rho }^{\prime }=\rho \left(\Re \left(\lambda \left(\alpha \right)\right)+\Re \left(C\left(\alpha \right)\right){\rho }^2+\Re \left(F\left(\alpha \right)\right){\rho }^4\right)$ (13)

To compute the saddle node bifurcation of the limit cycle, we analyze the equilibrium bifurcation of the normal form to get the condition

$\Delta =\Re \left(C\left(\alpha \right)\right)=0$ (14)

Example 1. A GH point is calculated at $K=12.30766$ , $a=0.68107986$ by Matcont as continuing Hopf bifurcation curve in $\left(K,a\right)$ plane. As shown in Figure 1, a limit point cycle (LPC) appears when we continue limit cycles with varying free parameter $K$ . In Figure 1(a), the subcritical Hopf bifurcation (the right Hopf point) occurs, which brings forth the hysteresis phenomena of limit cycles, while red color line denotes the branch of unstable limit cycles, and the blue color line represents the stable limit cycles. It is also seen that the super-critical Hopf bifurcation happens (the left Hopf point), which bifurcates the stable limit cycle. This explains that between two Hopf points, the generalized Hopf phenomena with the first Lyapunov exponent becoming zero may occur. If varying free parameter $a$ , the subcritical Hopf bifurcation occurs, which produces a limit point cycle (LPC) while continuing the limit cycles, as shown in Figure 1(b). As discussed above, Bautin bifurcation line Γ starts from GH point to be computed by Equation (14). Using Matcont, we also draw the green line of limit point cycle (LPC) bifurcation in $\left(K,a\right)$ parameter space, as shown in Figure 2. The near dynamics of GH point is classified into different regimes. In regime (1), the equilibrium solution is stable, however the equilibrium solution losses its stability as crossing Hopf line. Hence after, the equlibrium solution is unstable in regime (3)

Figure 1. The bifurcating periodical solutions from Hopf point, with chosen parameters $\theta =3$ , $\beta =0.5$ , $\delta =0.022$ . (a) The continue of limit cycles which arising from Hopf point as varying $K$ with $a=0.678$ ; both subcritical Hopf and super-critical Hopf bifurcation phenomena occur; (b) The continue of limit cycles with as varying $a$ with $K=16.3148$ .

Figure 2. The bifurcation diagram near GH point in $\left(K,a\right)$ plane, with chosen parameters $\theta =3$ , $\beta =0.5$ , $\delta =0.022$ . The near dynamics of GH point is classified into three different regimes topologically.

with a stable limit cyle observed. In regime (2), the unstable limit cycle occurs due to subcritical Hopf bifurctaion phenomena, and collides with the stable limit cycle on the curve Γ. The blue line in Figure 2 describes Hopf bifurcation phenomena.

3. TB Bifurcation Analysis

The coefficient matrix of the linearized system (5) is rewritten as

$A=\left(\begin{array}{cc}-{\left(\frac{x^{\text{*}}}{K}\right)}^{\theta }\theta +\frac{\beta x^{\text{*}}{y}^{\text{*}}}{{\left(1+x^{\text{*}}\right)}^2}& -\frac{\beta x^{\text{*}}}{1+x^{\text{*}}}\\ \frac{y^{\text{*}}}{1+x^{\text{*}}}-\frac{x^{\text{*}}{y}^{\text{*}}}{{\left(1+x^{\text{*}}\right)}^2}& -\delta y^{\text{*}}\end{array}\right)$

The sufficient condition for fold bifurcation is given as $det\left(A\right)=0$ , and TB bifurcation happens if the matrix $A$ has two zero roots. Set TB point $E=\left(K^{\text{*}},a^{\text{*}}\right)$ , suppose the right eigenvectors $q_{1,2}$ and left eigenvector $p_{1,2}$ satisfying

$A{q}_1=0,A{q}_2=q_1,A{p}_2=0,A{p}_1=p_2$ (15)

with $\langle p_1,q_1\rangle =1,\langle p_2,q_2\rangle =1$ . The multilinear form style of system (2) is written as

$X^{\prime }=A\left(\alpha \right)X+B\left(X,X\right)+C\left(X,X,X\right)+o\left({\left|X\right|}^3\right)$ (16)

Without loss generality, we set $X=q_1{z}_1+q_2{z}_2$ , with $\alpha =0$ , one easily get that

$\left(\begin{array}{c}{z^{\prime }}_1\\ {z^{\prime }}_2\end{array}\right)=\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right)\left(\begin{array}{c}{z}_1\\ z_2\end{array}\right)+\left(\begin{array}{l}+A_{02}{z}_2^2+A_{11}{z}_1{z}_2+A_{20}{z}_1^2+A_{03}{z}_2^3+A_{12}{z}_1{z}_2^2+A_{21}{z}_1^2{z}_2+A_{30}{z}_1^3\hfill \\ B_{03}{z}_2^3+B_{12}{z}_1{z}_2^2+B_{21}{z}_1^2{z}_2+B_{30}{z}_1^3+B_{02}{z}_2^2+B_{11}{z}_1{z}_2+B_{20}{z}_1^2\hfill \end{array}\right)$ (17)

with

$A_{20}=\langle p_1,B\left(q_1,q_1\right)\rangle ,A_{11}=2\langle p_1,B\left(q_1,q_2\right)\rangle ,A_{02}=\langle p_1,B\left(q_2,q_2\right)\rangle ,$

$A_{30}=\langle p_1,C\left(q_1,q_1,q_1\right)\rangle ,A_{21}=3\langle p_1,C\left(q_1,q_1,q_2\right)\rangle ,$

$A_{12}=3\langle p_1,C\left(q_1,q_2,q_2\right)\rangle ,A_{03}=\langle p_1,C\left(q_2,q_2,q_2\right)\rangle ,$

$B_{20}=\langle p_2,B\left(q_1,q_1\right)\rangle ,B_{11}=2\langle p_2,B\left(q_1,q_2\right)\rangle ,B_{02}=\langle p_2,B\left(q_2,q_2\right)\rangle ,$

$B_{30}=\langle p_2,C\left(q_1,q_1,q_1\right)\rangle ,B_{21}=3\langle p_2,C\left(q_1,q_1,q_2\right)\rangle ,$

$B_{12}=3\langle p_2,C\left(q_1,q_2,q_2\right)\rangle ,B_{03}=\langle p_2,C\left(q_2,q_2,q_2\right)\rangle .$

Setting

$\begin{array}{l}{y}_1=z_1,\\ y_2=z_2+A_{02}{z}_2^2+A_{11}{z}_1{z}_2+A_{20}{z}_1^2+A_{03}{z}_2^3+A_{12}{z}_1{z}_2^2+A_{21}{z}_1^2{z}_2+A_{30}{z}_1^3\end{array}$

one gets that

$\left(\begin{array}{c}{y^{\prime }}_1\\ {y^{\prime }}_2\end{array}\right)=\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right)\left(\begin{array}{c}{y}_1\\ y_2\end{array}\right)+\left(\begin{array}{c}0\\ a_1{y}_1^2+a_2{y}_1{y}_2+a_3{y}_2^2+a_4{y}_1^3+a_5{y}_1^2{y}_2+a_6{y}_1{y}_2^2+a_7{y}_2^3\end{array}\right)$ (18)

with

$\begin{array}{l}{a}_1=B_{20},a_2=\left(2A_{20}+B_{11}\right),a_3=\left(A_{11}+B_{02}\right),\\ a_4=\left(A_{11}{B}_{20}-A_{20}{B}_{11}+B_{30}\right),a_5=\left(2A_{02}{B}_{20}-A_{11}{A}_{20}-2A_{20}{B}_{02}+3A_{30}+B_{21}\right),\\ a_6=\left(A_{02}{B}_{11}-A_{11}^2-A_{11}{B}_{02}+2A_{21}+B_{12}\right),a_7=\left(-A_{02}{A}_{11}+A_{12}+B_{03}\right)\end{array}$

By time transform, the equivalent form of Equation (19) can be written as

$\left(\begin{array}{c}{y^{\prime }}_1\\ {y^{\prime }}_2\end{array}\right)=\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right)\left(\begin{array}{c}{y}_1\\ y_2\end{array}\right)+\left(\begin{array}{c}0\\ a_1{y}_1^2+a_2{y}_1{y}_2+\left(a_4-2a_1{a}_3\right)y_1^3\end{array}\right)$ (19)

Considering the linear part to get

$\left(\begin{array}{c}{z^{\prime }}_1\\ {z^{\prime }}_2\end{array}\right)=\left(\begin{array}{c}{x}_{00}\\ y_{00}\end{array}\right)+\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right)\left(\begin{array}{c}{z}_1\\ z_2\end{array}\right)+\left(\begin{array}{c}{\alpha }_1{z}_1+{\alpha }_2{z}_2\\ {\beta }_1{z}_1+{\beta }_2{z}_2\end{array}\right)$ (20)

with

$\begin{array}{l}{\alpha }_1=\langle p_1,A\left(\alpha \right)q_1\rangle ,{\alpha }_2=\langle p_1,A\left(\alpha \right)q_2\rangle ,\hfill \\ {\beta }_1=\langle p_2,A\left(\alpha \right)q_1\rangle ,{\beta }_2=\langle p_2,A\left(\alpha \right)q_2\rangle \hfill \end{array}$

Set

$\begin{array}{l}{y}_1=z_1,\\ y_2=x_{00}+{\alpha }_1{z}_1+{\alpha }_2{z}_2\end{array}$

one derives that

$\begin{array}{l}{y^{\prime }}_1=y_2\\ {y^{\prime }}_2={\alpha }_2{y}_{00}-{\beta }_2{x}_{00}+\left({\alpha }_1+{\beta }_2\right){\alpha }_2{y}_2+\left({\alpha }_2{\beta }_1-{\alpha }_1{\beta }_2\right)y_1\end{array}$ (21)

Adding the linear part of Equation (21) to Equation (19) to get

$\begin{array}{l}\left(\begin{array}{c}{y^{\prime }}_1\\ {y^{\prime }}_2\end{array}\right)=\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right)\left(\begin{array}{c}{y}_1\\ y_2\end{array}\right)+\left(\begin{array}{c}0\\ {\alpha }_2{y}_{00}-{\beta }_2{x}_{00}+\left({\alpha }_1+{\beta }_2\right){\alpha }_2{y}_2+\left({\alpha }_2{\beta }_1-{\alpha }_1{\beta }_2\right)y_1\end{array}\right)\\ +\left(\begin{array}{c}0\\ a_1{y}_1^2+a_2{y}_1{y}_2+\left(a_4-2a_1{a}_3\right)y_1^3\end{array}\right)\end{array}$ (22)

Doing the shift transformation in $y_1$ direction by $u_1=y_1+{\delta }_y,u_2=y_2$ to get

$\left(\begin{array}{c}{u^{\prime }}_1\\ {u^{\prime }}_2\end{array}\right)=\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right)\left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)+\left(\begin{array}{c}0\\ {\mu }_0+{\mu }_1{y}_1+a_1{y}_1^2+a_2{y}_1{y}_2+\left(a_4-2a_1{a}_3\right)y_1^3\end{array}\right)$ (23)

with

Figure 3. The homoclinic bifurcation with chosen parameters $\theta =3$ , $\beta =1.2$ , $\delta =0.2504$ . (a) The homoclinc bifurcation curve arising from BT point tangents to LP curve at point P; (b) The homoclinc solution with $K=7.0064$ , $a=0.2500$ ; (c) The homoclinc loop with a saddle $P_2$ and a node $P_3$ is plotted, with chosen $K=7.4297$ , $a=0.2405$ ; (d) The homoclinc loop near saddle node bifurcation point P is shown, wherein $K=7.1504$ , $a=0.2440$ .

Figure 4. The classification phase portraits near BT point. The near dynamics is classified into different regimes, from regime (1) to regime (4), as shown in Figure 3(a).

$\begin{array}{l}{\delta }_y=\frac{\left({\alpha }_1+{\beta }_2\right){\alpha }_2}{a_2},\\ {\mu }_0=\left({\alpha }_2{y}_{00}-{\beta }_2{x}_{00}\right),\\ {\mu }_1=\left({\alpha }_2{\beta }_1-{\alpha }_1{\beta }_2\right)-2\left({\alpha }_1+{\beta }_2\right){\alpha }_2\frac{a_1}{a_2}\end{array}$ (24)

By the universal norm form expressed by Equation (23), the near dynamics of TB point can be classified geometrically and topologically. By computing the equilibrium solution, the limit point line is calculated as

$L:{\mu }_1-4{\mu }_0{a}_1=0.$

and Hopf line is listed as

$H:{\mu }_1>0,{\mu }_0=0.$

The usual formula for the saddle homoclinic bifurcation line is

$Hom:{\mu }_0+\frac{6}{25}{\mu }_1^2=0.$

Example 1: With chosen parameters $\theta =3$ , $\beta =1.2$ , $\delta =0.2504$ , BT point is computed at $a=0.25295887$ , $K=7.7220041$ , and $E^{\text{*}}=\left(1.8093201,1.5618335\right)$ . As analyzed above, Homoclinc bifurcation starts from BT point. Using Matcont, the bifurcation diagram of system (2) is drawn as shown in Figure 3(a). We produce the homoclinic solution as shown in Figures 3(b)-(d). The homoclinic bifurcation curve starts from BT point in system (2) by Matcont, we notice that the homoclinic curve tangents to the limit curve LP at point P, as shown in Figure 3(a). A homoclinic solution with saddle node lying on loop is observed. As shown in Figure 3(a), The continuation of homoclinic solutions is drawn as varying free parameter $a$ and $K$ . The homoclinic solution with saddle is observed in Figure 3(b); whilst the homoclinc loop arises with saddle and node is manifested in Figure 3(c); It is concluded that saddle-node bifurcation of equilibrium solution happens at the homoclinic loop as varying free parameter along the homoclinc curve at point P. We also draw homoclinc solution near the point P, as shown in Figure 3(d).

Figure 5. The dynamical bifurcation near BT point with chosen parameters $\theta =3$ , $\beta =1.2$ , $\delta =0.12$ . (a) The bifurcation diagram near BT point; (b) The unstable periodical solution observed in regime (3) as shown in Figure 5(a); (c) The homoclinc loop with a saddle and a node on the loop, chosen $K=10.2859$ , $a=0.4894$ ; (d) The homoclinc loop with a saddle as $K=10.4953$ , $a=0.4884$ .

The classification diagram near BT point is drawn on $\left(K,a\right)$ parameter plane, as shown in Figure 3(a), which leads to the near dynamics of BT point being analyzed topologically and geometrically. The classification regime near BT point is denoted by regimes (1) - (4), and the phase portraits in different regimes are drawn in Figure 4.

Example 2. We calculate that BT bifurcation happens at $K=10.310284$ , $a=0.48918857$ with $\theta =3$ , $\beta =1.2$ , $\delta =0.12$ , and system (2) has equilibrium solution $x^{\text{*}}=4.0613503$ , $y^{\text{*}}=2.6102975$ . Using Matcont, we compute the homoclinic curve from BT point, and the near dynamics of BT point is analyzed. As discussed above, Hopf bifurcation curve $H$ tangents to the limit point curve LP at BT point. The homoclinc curve starts from BT point and the homoclinc solutions are simulated. With free parameter values $K=10.2859$ , $a=0.4894$ , the homoclinc solution with a saddle and a node on the loop is simulated, as shown in Figure 5(c); Whilst the homoclinc solution with a saddle is simulated in Figure 5(d). It is seen a saddle node bifurcation of equilibrium solution happens at the tangent point while the homoclinic curve tangents to the curve LP. In Figure 5(b), an unstable periodical solution is observed in regime (3). The near dynamics near BT point is classified into different regimes, and the phase portraits of regimes (1) - (5) are classified topologically and geometrically, as shown in Figure 6.

Figure 6. The phase portrait is classified near BT point topologically and geometrically. The near dynamics of BT point is classified different regimes in Figure 5(a).

4. Bifurcation Diagram

With fixed parameter $\theta =3,\beta =1.2$ , we understand that Hopf bifurcation often happens in system (2), therefore we drawn Hopf surface in three-dimensional parameter space $\left(K,a,\delta \right)$ , as shown in Figure 7. The sufficient condition for Hopf bifurcation is $Tr\left(A\right)=0$ , which is listed in Equation (7), and the sufficient condition for fold bifurcation is $det\left(A\right)=0$ as given in Equation (6). Furthermore, BT bifurcation occurs as $Tr\left(A\right)=0$ , $det\left(A\right)=0$ . Hence, we also drew BT bifurcation curve onto Hopf surface, noted as H-surface in Figure 7. As discussed in section 2, the generalized Hopf bifurcation occurs if the first Lyapunov coefficient becomes zero, which means GH point adjoint the super-critical Hopf curve to the sub-critical Hopf curve, and we also calculate GH bifurcation curve on H-surface. It is seen that when varying free parameter $\delta$ , GH bifurcation behavior arises as $\delta$ below its super limit value. Hence, BT and GH intersects at a point. As discussed in Section 3, the two examples conclude that the mechanism that leads to homoclinic solution arising from BT point are different. The homoclinic solution appears in Example 1 and the stable periodical solution is broken as free parameter crossing over the homoclinic curve. However, the homoclinc solution appears in Example 2 and the unstable limit cycle bifurcates from the stable equilibrium solution as free parameter across the homoclinic curve. Therefore, we deduce that the above two examples show that GH bifurcation changes the stability of homoclinc solution. The Homoclinc solution is stable if the inside equilibrium solution is unstable. However, the homoclinc solution is unstable as GH happens near BT point.

Figure 7. The Hopf surface in parameter space $\left(K,a,\delta \right)$ is drawn with fixed parameters $\theta =3$ , $\beta =1.2$ .

5. Discussion

In two-species population model, prey species was assumed with $\theta$ logistic growth rate and the population dynamics was complex. We discussed that Hopf bifurcation often appears with free parameters varying in $\left(K,a,\delta \right)$ space. Hence, Hopf surface was drawn. It was seen that in addition to BT bifurcation, which arises at the tangent point while Hopf curve tangents to the limit point line LP, GH Hopf bifurcation point was found on Hopf surface as the first Lyapunov coefficient became zero. Chosen fixed $\delta$ , the bifurcation diagram was drawn in $\left(K,a\right)$ parameter plane and the homoclinc bifurcation was discussed.

Conflicts of Interest

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

References

[1]	Cai, L., Chen, G. and Xiao, D. (2012) Multiparametric Bifurcations of an Epidemiological Model with Strong Allee Effect. Journal of Mathematical Biology, 67, 185-215.[CrossRef] [PubMed]
[2]	Hsu, S.B. and Huang, T.W. (1999) Hopf Bifurcation Analysis for a Predator-Prey System of Holling and Leslie Type. Taiwanese Journal of Mathematics, 3, 35-53.[CrossRef]
[3]	Gasull, A., Kooij, R.E. and Torregrosa, J. (1997) Limit Cycles in the Holling-Tanner Model. Publicacions Matemàtiques, 41, 149-167.[CrossRef]
[4]	Lamontagne, Y., Coutu, C. and Rousseau, C. (2008) Bifurcation Analysis of a Predator-Prey System with Generalised Holling Type III Functional Response. Journal of Dynamics and Differential Equations, 20, 535-571.[CrossRef]
[5]	Xiang, C., Huang, J. and Wang, H. (2022) Linking Bifurcation Analysis of Holling-Tanner Model with Generalist Predator to a Changing Environment. Studies in Applied Mathematics, 149, 124-163.[CrossRef]
[6]	Xiang, C., Huang, J., Ruan, S. and Xiao, D. (2020) Bifurcation Analysis in a Host-Generalist Parasitoid Model with Holling II Functional Response. Journal of Differential Equations, 268, 4618-4662.[CrossRef]
[7]	Beddington, J.R. and May, R.M. (1982) The Harvesting of Interacting Species in a Natural Ecosystem. Scientific American, 247, 62-69.[CrossRef]
[8]	Brauer, F. and Soudack, A.C. (1979) Stability Regions in Predator-Prey Systems with Constant-Rate Prey Harvesting. Journal of Mathematical Biology, 8, 55-71.[CrossRef]
[9]	Brauer, F. and Soudack, A.C. (1982) Coexistence Properties of Some Predator-Prey Systems under Constant Rate Harvesting and Stocking. Journal of Mathematical Biology, 12, 101-114.[CrossRef]
[10]	Chowdhury, J., Al Basir, F., Mukherjee, A. and Roy, P.K. (2025) A Theta Logistic Model for the Dynamics of Whitefly Borne Mosaic Disease in Cassava: Impact of Roguing and Insecticide Spraying. Journal of Applied Mathematics and Computing, 71, 4897-4914.[CrossRef]
[11]	Chow, S.N., Li, C. and Wang, D. (1994) Normal Forms and Bifurcation of Planar Vector Fields. Cambridge University Press.[CrossRef]
[12]	Kuznetsov, Yu.A. (1998) Elements of Applied Bifrucation Theory. 2nd Edition, Springer-Verlag.
[13]	Kuznetsov, Yu.A. (2011) Practical Computation of Normal Forms on Center Manifolds at Degenerate Bogdanov-Takens Bifurcations. International Journal of Bifurcation & Chaos, 15, 3535-3546.
[14]	Dhooge, A., Govaerts, W. and Kuznetsov, Y.A. (2003) MATCONT: A MATLAB Package for Numerical Bifrucation Analysis of ODEs. ACM Transactions on Mathematical Software, 29, 141-164.[CrossRef]
[15]	Doedel, E.J., Champneys, A.R., Fairgrieve, T.F., Kuznetsov, Yu.A., Sandstede, B. and Wang, X.J. (2000) Auto97-Auto2000: Continuation and Bifurcation Software for Ordinary Differential Equations (with HomCont), User’s Guide. Concordia University. http://indy.cs.concordia.ca
[16]	Govaerts, W.J.F. (2000) Numerical Methods for Bifurcations of Dynamical Equilibria. Society for Industrial and Applied Mathematics.[CrossRef]