To determine the equilibria of system Eq. 2, we analyze the prey and predator nullclines, which are given by the following equations: x 1 − x − x y 1 + a x = 0 , b y 1 − d y x = 0 .

From these nullclines, it is evident that system Eq. 2 has a unique boundary equilibrium at E 0 0 = ( 1,0 ) . To identify potential positive equilibria, we solve the following equations: f x ≔ 1 − x − x d 1 + a x = 0 , y ≔ x d .

By substituting x = z , we transform the problem of finding positive solutions to f(x) = 0 into finding positive solutions to g(z) = adz ³ + dz ² + (1 − ad)z − d = 0. Through algebraic analysis, we can easily determine that g(z) = 0 has a unique positive solution. Consequently, system Eq. 2 has a unique positive equilibrium point, denoted as E 0 * = ( x 0 * , y 0 * ) , where y 0 * = x 0 * d , x 0 * = z 0 2 , and z ₀ is the unique positive root of g(z) = 0. Based on this analysis, we can formulate the following theorem regarding the number of equilibria in system Eq. 2.

Theorem 2.1

System Eq. 2 possesses a unique boundary equilibrium E 0 0 = ( 1,0 ) and unique positive equilibrium denoted by E 0 * = ( x 0 * , y 0 * ) .

Theorem 2.2

For all positive parameters, system Eq. 2 has a boundary equilibrium E 0 0 and a positive equilibrium E 0 * . Moreover, the boundary equilibrium E 0 0 is always a hyperbolic saddle. Furthermore, the positive equilibrium E 0 * is a sink when b > b ₀, a source when b < b ₀, or a center when b = b ₀, where b 0 = 1 − z 0 2 − z 0 2 d ( 1 + a z 0 ) 2 and z ₀ is the positive root of g(z) = 0.

Proof. The Jacobian matrix of system Eq. 2 evaluated at E 0 0 is given by J E 0 0 = − 1 − 1 1 + a 0 b .

It is obvious that the matrix J ( E 0 0 ) has a positive eigenvalue b and a negative eigenvalue −1. Therefore, E 0 0 is always a hyperbolic saddle.

The Jacobian matrix of system Eq. 2 evaluated at E 0 * is given by J E 0 * = 1 − 2 x 0 * − x 0 * 2 d 1 + a x 0 * 2 − x 0 * 1 + a x 0 * b d − b = 1 − 2 z 0 2 − z 0 2 d 1 + a z 0 2 − z 0 1 + a z 0 b d − b .

According to g(z ₀) = 0, we can obtain 1 − z 0 2 = z 0 d ( 1 + a z 0 ) . Then, D e t J E 0 * = b z 0 2 + b z 0 2 d 1 + a z 0 2 > 0 , T r J E 0 * = 1 − z 0 2 − z 0 2 d 1 + a z 0 2 − b ≜ b 0 − b .

Hence, when b > b ₀, T r [ J ( E 0 * ) ] < 0 , the matrix J ( E 0 * ) has two negative part eigenvalues, and E 0 * is a sink. When b < b ₀, then T r [ J ( E 0 * ) ] > 0 , the matrix J ( E 0 * ) has two positive part eigenvalues, and E 0 * is a source. When b = b ₀, then T r [ J ( E 0 * ) ] = 0 , the matrix J ( E 0 * ) has a pair of purely imaginary eigenvalues, and E 0 * is a center.

This section conducts a qualitative analysis of model Eq. 3. With the biological context in mind, we focus on examining the dynamics of system Eq. 3 within the closed first quadrant, denoted as R + 2 , in the (x, y) coordinate plane.

In order to obtain the equilibria of system Eq. 3, we first study the roots of the following equations: x 1 − x − x y 1 + a x − h = 0 , b y 1 − d y x = 0 . (4)

Regarding the count of boundary equilibria for system Eq. 3, we derived the subsequent theorem.

Theorem 3.1

(i) When h = 1 4 , system Eq. 3 admits a unique boundary equilibrium E 0 = ( 1 2 , 0 ) .

(ii) When 0 < h < 1 4 , system Eq. 3 admits two boundary equilibria E _i = (x _i , 0) (i = 1, 2), where x 1 = 1 − 1 − 4 h 2 and x 2 = 1 + 1 − 4 h 2 .

Proof. When h exceeds 1 4 , x ̇ = x ( 1 − x ) − h − x y 1 + a x < 0 . This implies that the dynamics of system Eq. 3 within R + 2 are inherently trivial. Consequently, the prey species faces extinction, an event that subsequently triggers the extinction of the predator. It is evident that system Eq. 3 lacks any equilibrium points within the positive quadrant R + 2 .

We next explore the existence of boundary equilibria under the condition that h ≤ 1 4 . Specifically, when y = 0 in Eq. 4, we are confronted with the following equation: x 2 − x − h = 0 . (5)

It is straightforward to determine that when h < 1 4 , Eq. 5 has two distinct roots, denoted as x ₁ and x ₂. Conversely, when h = 1 4 , the equation yields a unique root, which is 1 2 .

Note that f 1 z ≜ a d z 5 + d z 4 + 1 − a d z 3 − d z 2 + a d h z + h d , g 1 z ≜ f 1 ′ z = 5 a d z 4 + 4 d z 3 + 3 1 − a d z 2 − 2 d z + a d h , h 1 z ≜ g 1 ′ z = 20 a d z 3 + 12 d z 2 + 6 1 − a d z − 2 d . (6)

Regarding the positive equilibria of system Eq. 3, we have the following theorem.

Theorem 3.2

(i) When h = h 2 * and h 2 * < h 1 * , system Eq. 3 has a unique positive equilibrium E ₃ = (x ₃, y ₃).

(ii) When h < min { h 1 * , h 2 * } , system (4) has two different positive equilibria E _i = (x _i , y _i ) (i = 4, 5),

where h 1 * = 5 a z 1 4 + 2 z 1 3 + z 1 a , h 2 * = z 3 2 ( 2 a z 3 3 + z 3 2 + 1 ) 2 a z 3 + 3 , z ₃ is the larger positive root between two positive roots of g ₁(z) = 0, and z ₁ is the only positive root of h ₁(z) = 0.

Proof. When y ≠ 0, then y = x d according to Eq. 5, and x is the positive root of the equation x 1 − x 1 + a x d − x x − d h 1 + a x = f 1 z , where z = x and f ₁(z) is given in Eq. 6.

Through algebraic analysis, we find that the equation h ₁(z) = 0 has a unique positive root, denoted as z ₁. When g ₁(z ₁) < 0, which corresponds to h < h 1 * ≜ 5 a z 1 4 + 2 z 1 3 + z 1 a , the equation g ₁(z) = 0 admits two distinct positive roots, z ₂ and z ₃ (with z ₂ < z ₃). If f ₁(z ₃) = 0, which occurs when h = h 2 * ≜ z 3 2 ( 2 a z 3 3 + z 3 2 + 1 ) 2 a z 3 + 3 , then the equation f ₁(z) = 0 has a unique positive root. Consequently, system Eq. 3 admits a unique positive equilibrium E ₃ = (x ₃, y ₃), where x 3 = z 3 2 and y 3 = x 3 d . If f ₁(z ₃) < 0, corresponding to h < h 2 * , then the equation f ₁(z) = 0 has two distinct positive roots, z ₄ and z ₅ (with z ₂ < z ₄ < z ₃ < z ₅). In this case, system Eq. 3 has two distinct positive equilibria, E ₄ = (x ₄, y ₄) and E ₅ = (x ₅, y ₅), where x 4 = z 4 2 , x 5 = z 5 2 , y 4 = x 4 d , and y 5 = x 5 d .

Theorem 3.3

When h = 1 4 , system Eq. 3 possesses a unique boundary equilibrium at E 0 = ( 1 2 , 0 ) , which is classified as a saddle-node. This equilibrium divides a sufficiently small neighborhood into two parts, separated by two separatrices that approach E ₀ from above and below. One part forms a parabolic sector, while the other comprises two hyperbolic sectors. Additionally, the parabolic sector is situated in the left half-plane, and the equilibrium E ₀ acts as a repelling saddle-node. The dynamics of the system in this scenario are illustrated in the phase portraits shown in Figure 2A.

Proof. The Jacobian matrix J(E ₀) of system Eq. 3 evaluated at E ₀ is J E 0 = 0 − 1 2 + a 0 b .

Given that Det[J(E ₀)] = 0, we can conclude that the equilibrium ( 1 2 , 0 ) is nonhyperbolic, indicating that it is a degenerate equilibrium. It is evident that the eigenvalues of the Jacobian matrix J(E ₀) are λ ₁ = 0 and λ ₂ = b. To ascertain the stability of E ₀, we perform a transformation by shifting the equilibrium to the origin using the transformation ( X , Y ) = ( x − 1 2 , y ) . We then expand system Eq. 3 in a power series up to the second order around the origin. Under this transformation, system Eq. 3 becomes X ̇ = − 1 2 + a Y − X 2 − 2 2 + a 2 X Y + P 1 X , Y , Y ̇ = b Y − 2 b d Y 2 + Q 1 X , Y , where P ₁(X, Y) and Q ₁(X, Y) are smooth functions of at least the third order with respect to (X, Y).

Then, making the transformation X Y = 1 − 1 2 + a 0 b x y , and introducing a new time variable τ = bt, for which we retain t to denote τ for notational simplicity, we obtain x ̇ = a 20 x 2 + a 11 x y + a 02 y 2 + P 2 x , y , y ̇ = y − 2 b d y 2 + Q 2 x , y , where P ₂(x, y) and Q ₂(x, y) are smooth functions of at least the third order with respect to (x, y) and a 20 = − a 5 + 5 2 a 4 + 20 a 3 + 20 2 a 2 + 20 a + 4 2 a + 2 5 b < 0 , a 11 = 2 a 4 − 2 b + 8 a 3 + 6 4 − b a 2 − 2 2 3 b − 8 a + 4 2 − b a + 2 5 b , a 02 = 2 a 2 b 2 d a 2 + 4 2 a + 12 + a 2 a − 2 b + 3 2 + a 16 b 2 d + 6 − 4 b + 8 b 2 d − 2 2 b + 2 2 a + 2 5 b .

According to Theorem 7.1 from Chapter 2 in [23], the equilibrium E ₀ is classified as a saddle-node. This implies that the vicinity of E ₀ is bifurcated by two separatrices that approach E ₀ from above and below. One region forms a parabolic sector, while the other comprises two hyperbolic sectors. Additionally, the parabolic sector is situated in the right half-plane, and E ₀ acts as a repelling saddle-node.

Theorem 3.4. When 0 < h < 1 4 , system Eq. 3 has two different boundary equilibria E ₁ and E ₂ . Moreover, E ₁ is always a hyperbolic unstable node, and E ₂ is always a hyperbolic saddle. The corresponding phase portraits are shown in Figure 2B and Figure 3.

Proof. The Jacobian matrix J(E ₁) of system Eq. 3 evaluated at E ₁ is J E 1 = 1 − 4 h − x 1 1 + a x 1 0 b .

We can obtain that the eigenvalues are λ 1 = 1 − 4 h > 0 and λ ₂ = b > 0. Hence, E ₁ is a hyperbolic unstable node.

The Jacobian matrix J(E ₂) of system Eq. 3 evaluated at E ₂ is J E 2 = − 1 − 4 h − x 2 1 + a x 2 0 b .

It is easy to observe that the eigenvalues are λ 1 = − 1 − 4 h < 0 and λ ₂ = b > 0. Hence, E ₂ is a hyperbolic saddle.

Now, we investigate the stability of equilibrium E ₃. The Jacobian matrix J(E ₃) of system Eq. 3 evaluated at E ₃ is J E 3 = 1 − 2 x 3 − x 3 2 d 1 + a x 3 2 − x 3 1 + a x 3 b d − b ≜ a 10 a 01 b 10 b 01 .

Straightforward computation shows that, due to x 3 = z 3 2 , D e t J E 3 = b 2 d 1 + a z 3 2 l z 3 , where l z = 4 a 2 d z 4 + 8 a d z 3 + 2 a + 4 d − 2 a 2 d z 2 + 3 − 4 a d z − 2 d .

Since f(z ₃) = 0, from Theorem 3.2, we can obtain h = − z 3 2 a d z 3 3 + d z 3 2 + 1 − a d z 3 − d d 1 + a z 3 .

Substituting g ₁(z ₃) in the above equation yields g 1 ( z 3 ) = z 3 1 + a z 3 l ( z 3 ) . Because z ₃ is a positive root of g ₁(z) = 0, g ₁(z ₃) = 0 and l(z ₃) = 0. Therefore, Det[J(E ₃)] = 0. Therefore, we can establish that equilibrium E ₃ is nonhyperbolic, which means it is a degenerate equilibrium. To elaborate further, we present the following theorem:

Theorem 3.5

(i) If a 10 + b 01 ≠ 0 ( b ≠ b 1 * ≜ 1 − 2 x 3 − x 3 2 d ( 1 + a x 3 ) 2 ) ,

(i-1) a 10 + b 01 < 0 ( b > b 1 * ) , then the parabolic sector is located in the right half-plane, and E ₃ is an attracting saddle-node. The corresponding phase portrait is shown in Figure 3B.

(i-2) a 10 + b 01 > 0 ( b < b 1 * ) , then the parabolic sector is located in the left half-plane, and E ₃ is a repelling saddle-node. The corresponding phase portrait is shown in Figure 3A.

(ii) When a 10 + b 01 = 0 ( b = b 1 * ) , E ₃ represents a degenerate critical point, specifically a cusp. Furthermore, E ₃ is a cusp of co-dimension 2 under the condition that d ≠ d 1 * and d ≠ d 2 * . However, if d = d 1 * or d = d 2 * , E ₃ is a cusp of co-dimension at least 3. The corresponding phase portraits are shown in Figure 2B, where d 1 * = − 4 a x 3 + ( 3 a − 4 ) x 3 + 1 8 x 3 ( 1 + a x 3 ) 3 and d 2 * = 6 x 3 ( 1 + a x 3 ) + 3 a x 3 + 1 8 ( 1 + a x 3 ) 3 ( 2 − 3 x 3 ) .

Proof. Since Tr[J(E ₃)] = a ₁₀ + b ₀₁, we know that the eigenvalues of J(E ₃) are λ ₁ = 0 and λ ₂ = a ₁₀ + b ₀₁.

Similar to the proof outlined in Theorem 3.3, we perform a transformation that shifts equilibrium E ₃ to the origin. We then expand system Eq. 3 in a power series up to the second order around this new origin. This transformation yields a simplified system that X ̇ = a 10 X + a 01 Y + a 11 X Y + a 20 X 2 + P 11 X , Y , Y ̇ = b 10 X + b 01 Y + b 11 X Y + b 20 X 2 + + b 02 Y 2 + Q 11 X , Y , (7) where P ₁₁(X, Y) and Q ₁₁(X, Y) are smooth functions of at least the third order with respect to (X, Y), and a 11 = − 1 2 x 3 1 + a x 3 2 , a 20 = 1 + 3 a x 3 8 d x 3 1 + a x 3 3 − 1 , b 11 = 2 b x 3 , b 20 = − b d x 3 , b 02 = − b d x 3 .

When a ₁₀ + b ₀₁ ≠ 0, J(E ₃) has one zero eigenvalue and a nonzero eigenvalue. To further analyze the stability of this equilibrium, we introduce a transformation that X Y = 1 a 10 b 10 − a 10 a 01 1 x y , and introducing a new time variable τ ₁ = (a ₁₀ + b ₀₁)t, for which we retain t to denote τ ₁ for notational simplicity, we obtain x ̇ = c 20 x 2 + c 11 x y + c 02 y 2 + P 21 x , y , y ̇ = y + d 20 x 2 + d 11 x y + d 02 y 2 + Q 21 x , y , where P ₂₁(x, y) and Q ₂₁(x, y) are smooth functions of at least the third order with respect to (x, y), and c 20 = 1 a 10 + b 01 2 a 01 a 11 − b 02 + a 10 a 01 a 20 − b 11 b 10 − a 10 2 a 01 b 20 b 10 2 , c 11 = 1 a 10 + b 01 2 a 11 b 01 − a 10 a 11 − 2 b 02 − 2 b 01 b 20 a 10 2 b 10 2 + b 11 a 10 2 + 2 a 10 b 01 a 20 − 1 2 b 11 b 10 , c 02 = 1 a 10 + b 01 2 a 10 b 11 − a 01 b 20 + a 20 b 01 − a 11 b 10 − a 10 b 10 b 02 b 01 , d 20 = 1 a 10 + b 01 2 b 01 b 20 + b 10 a 20 − b 11 − b 10 2 a 11 − b 02 b 01 , d 11 = 1 a 10 + b 01 2 2 a 01 b 20 + b 01 b 11 + a 11 − 2 b 02 b 10 + 2 a 20 − b 11 a 10 − a 10 a 11 b 10 b 01 , d 02 = 1 a 10 + b 01 2 a 10 a 11 + b 01 b 02 + a 01 b 11 + a 10 2 a 20 + a 10 a 01 b 20 b 10 .

Straightforward computation shows that c 20 = b a 10 + b 01 2 ⋅ 3 a 2 d x 3 3 2 + a 8 + d x 3 + a 3 d x 3 2 + 3 a d x 3 + 3 8 d x 3 a + a x 3 3 > 0 .

Consequently, as shown by Theorem 7.1 from Chapter 2 in [23], if a ₁₀ + b ₀₁ ≠ 0, then the equilibrium E ₃ qualifies as a saddle-node. Moreover, the characteristics of this saddle-node are determined as follows: (i) if a ₁₀ + b ₀₁ > 0, the parabolic sector is situated in the upper half-plane, and E ₃ acts as a repelling saddle-node and (ii) conversely, if a ₁₀ + b ₀₁ < 0, the parabolic sector remains in the upper half-plane, but E ₃ is an attracting saddle-node.

When a ₁₀ + b ₀₁ = 0, J(E ₃) has one zero eigenvalue with two multiple. Then, making the transformation X Y = 1 1 − a 10 a 01 0 x y , and introducing a new time variable τ ₂ = −b ₀₁ t = bt. For notational simplicity, we retain t to represent τ ₂, then system Eq. 7 becomes x ̇ = y + e 02 y 2 + P 22 x , y , y ̇ = f 20 x 2 + f 11 x y + f 02 y 2 + Q 22 x , y , (8) where P ₂₂(x, y) and Q ₂₂(x, y) are smooth functions of at least the third order with respect to (x, y) and e 02 = − 1 x 3 , f 20 = a 11 b d + a 20 b , f 11 = 2 a 20 b + a 11 b d − 4 x 3 , f 02 = a 20 b + 1 x 3 .

To eliminate y ² terms in system Eq. 8, we make the following near identity transformation: x = X + f 02 2 X 2 , y = Y + f 02 X Y − e 02 Y 2 , which transforms system Eq. 8 into X ̇ = Y + P 23 X , Y , Y ̇ = f 20 X 2 + f 11 X Y + Q 23 X , Y , where P ₂₃(X, Y) and Q ₂₃(X, Y) are smooth functions of at least the third order with respect to (X, Y). Therefore, when f ₂₀ ≠ 0 ( d ≠ d 1 * = − 4 a x 3 + ( 3 a − 4 ) x 3 + 1 8 x 3 ( 1 + a x 3 ) 3 ) and f ₁₁ ≠ 0 ( d ≠ d 2 * = 6 x 3 ( 1 + a x 3 ) + 3 a x 3 + 1 8 ( 1 + a x 3 ) 3 ( 2 − 3 x 3 ) ) , according to the results presented by [24], E ₃ is a cusp of co-dimension 2. If f 20 = 0 ( d = d 1 * ) or f ₁₁ = 0 ( d = d 2 * ) , E ₃ is a cusp of co-dimension at least 3.

Theorem 3.6

When system Eq. 3 exists with two positive equilibria E ₄ and E ₅, E ₄ is always a hyperbolic saddle. Furthermore, positive equilibrium E ₅ is a sink when b > b*, a source when b < b*, or a center when b = b*, where b * = 1 − 2 z 5 2 − z 5 2 d ( 1 + a z 5 ) 2 and z ₅ is the larger positive root of f ₁(z) = 0.

Proof. The Jacobian matrix J(E ₄) of system Eq. 3 evaluated at E ₄ is J E 4 = 1 − 2 x 4 − x 4 2 d 1 + a x 4 2 − x 4 1 + a x 4 b d b .

Straightforward computation shows that, due to x 4 = z 4 2 , D e t J E 4 = b 2 d 1 + a z 4 2 l z 4 , where l z = 4 a 2 d z 4 + 8 a d z 3 + 2 a + 4 d − 2 a 2 d z 2 + 3 − 4 a d z − 2 d .

Since f(z ₄) = 0, from Theorem 3.2, we obtain h = − z 4 2 a d z 4 3 + d z 4 2 + 1 − a d z 4 − d d 1 + a z 4 .

Substituting g ₁(z ₄) by the above equation, then g 1 ( z 4 ) = z 4 1 + a z 4 l ( z 4 ) . Because z ₂ < z ₄ < z ₃, z ₂ and z ₃ are two positive roots of g ₁(z) = 0, then g ₁(z ₄) < 0 and l(z ₄) < 0. Therefore, Det[J(E ₄)] < 0. Hence, the matrix J(E ₄) has a positive eigenvalue and a negative eigenvalue. Hence, E ₄ is a hyperbolic saddle.

Similarly, we obtain D e t J E 5 = b 2 d 1 + a z 5 2 l z 5 > 0 , since l(z ₅) > 0 (z ₅ > z ₃). Moreover, we can also obtain T r J E 5 = 1 − 2 x 5 − x 5 2 d 1 + a x 5 2 − b ≜ b * − b .

Hence, when b > b*, Tr[J(E ₅)] < 0, the matrix J(E ₅) has two negative part eigenvalues, and E ₅ is a sink. When b < b*, then Tr[J(E ₅)] > 0, the matrix J(E ₅) has two positive part eigenvalues, and E ₅ is a source. When b = b*, then Tr[J(E ₅)] = 0, the matrix J(E ₅) has a pair of pure imaginary eigenvalues, and E ₅ is a center.

According to Theorem 3.1, we provided the conditions for the existence of E ₁ and E ₂ based on some restrictions. We obtain that S N 1 = a , b , d , h : h = 1 4 , a > 0 , b > 0 , d > 0 is a saddle-node bifurcation surface. Consequently, we can deduce that on the surface SN ₁, system Eq. 3 possesses a unique equilibrium E ₀, which is a saddle-node. As the parameters traverse this surface, the quantity of boundary equilibria undergoes a transition from zero to two. Biologically, this bifurcation corresponds to a critical threshold where the maximum sustainable yield (MSY) is at h MSY = 1 4 . Beyond this threshold ( h > 1 4 ) , the prey species faces extinction, leading to the collapse of the system. However, under certain initial conditions where 0 < h < 1 4 , the prey species can avoid extinction [18, 20]. This phenomenon is shown in the saddle-node bifurcation diagram, as shown in Figure 4B.

Similar to the above discussion, we obtain S N 2 = a , b , d , h : h = h 2 * , h 2 * < h 1 * , a > 0 , b > 0 , d > 0 , according to Theorem 3.2, is another saddle-node bifurcation surface. When the parameters transition from one side of the surface to the other, the number of positive equilibria in system Eq. 3 changes from zero to two. The saddle-node bifurcation is illustrated in Figure 4A.

From Theorem 3.6, we know that the positive equilibrium E ₅ is a center-type nonhyperbolic equilibrium when b = b* and h < min { h 1 * , h 2 * } . Therefore, system (4) may undergo a Hopf bifurcation. To proceed, we first confirm the transversality condition necessary for a Hopf bifurcation, d d b T r J E 5 | b = b * = − 1 ≠ 0 .

Therefore, the stability of the positive equilibrium E ₅ undergoes a change as the parameters traverse the specific critical surface H 1 = a , b , d , h : b = b * , h < min h 1 * , h 2 * , a > 0 , d > 0 .

To determine the direction of the Hopf bifurcation, we calculate the first Lyapunov number l ₁ at the equilibrium E ₅ = (x ₅, y ₅). Initially, we shift the equilibrium to the origin using the transformation ( x ̄ , y ̄ ) = ( x − x 5 , y − y 5 ) . Subsequently, we expand system Eq. 3 into a power series around the origin. Then, system Eq. 3 becomes d x ̄ d t = α 10 x ̄ + α 01 y ̄ + α 20 x ̄ 2 + α 11 x ̄ y ̄ + α 30 x ̄ 3 + α 21 x ̄ 2 y ̄ + P 21 x ̄ , y ̄ , d y ̄ d t = β 10 x ̄ + β 01 y ̄ + β 20 x ̄ 2 + β 11 x ̄ y ̄ + β 02 y ̄ 2 + β 30 x ̄ 3 + β 21 x ̄ 2 y ̄ + β 12 x ̄ y ̄ 2 + Q 21 x ̄ , y ̄ , where P 21 ( x ̄ , y ̄ ) and Q 21 ( x ̄ , y ̄ ) are smooth functions of at least the fourth order with respect to ( x ̄ , y ̄ ) and α 10 = 1 − 2 x 5 − x 5 2 d 1 + a x 5 2 , α 01 = − x 5 1 + a x 5 , α 11 = − 1 2 x 5 1 + a x 5 2 , α 20 = 1 + 3 a x 5 8 d x 5 1 + a x 5 3 − 1 , α 30 = − 5 a 2 x 5 + 4 a x 5 + 1 16 d x 5 3 2 1 + a x 5 4 , α 21 = 1 + 3 a x 5 8 x 5 3 2 1 + a x 5 3 , β 10 = b d , β 01 = − b , β 20 = − b d x 5 , β 11 = 2 b x 5 , β 02 = − b d x 5 , β 30 = b d x 5 2 , β 21 = − 2 b x 5 2 , β 12 = b d x 5 2 .

In accordance with the formula for the first Lyapunov number l ₁ presented in [24] [p.353], we proceed with the computation l 1 = − 3 π 2 α 01 Δ 3 2 M , where M = α 10 α 11 β 10 α 11 + β 02 + α 01 β 11 2 + α 20 β 11 + α 11 β 02 − 2 β 10 β 02 2 − 2 α 01 α 20 2 − β 20 β 02 − α 01 2 2 α 20 β 20 + β 11 β 20 + α 01 β 10 − 2 α 10 2 β 11 β 02 − α 11 α 20 − α 10 2 + α 01 β 10 − 3 α 01 α 30 + 2 α 10 α 21 + β 12 − β 21 β 01 , and Δ = α 10 β 01 − α 01 β 10 = D e t J E 5 > 0 , α 01 = − x 5 1 + a x 5 < 0 .