The uniqueness of limit cycles in a predator-prey system with Ivlev-type group defense

1.   Introduction

In this article, we study the limit cycle problem for the following Gause predator-prey system introduced in ^[26] by Xiao and Ruan:

where all parameters are positive: $r > 0$, $\beta > 0$, $k > 0$, $D > 0$, $\mu > 0$. For convenience, we will write $x$ and $y$ for $x(t)$ and $y(t)$, respectively. The real variables $x$ and $y$ represent the nonnegative density of the prey and predator population, respectively.

The main result of this paper is the proof of the conjecture posed in ^[26] concerning the number of limit cycles in (1.1).

Theorem 1.1. System (1.1) has at most one limit cycle. If it exists, it is stable and hyperbolic.

In a recent paper ^[15], system (1.1) was discussed for the special case of a singular perturbation, essentially studying the occurrence of so-called canard cycles for small $D$ and $\mu$.

System (1.1) is a special case of the more general Gause predator-prey system ^[5]:

In system (1.1), the natural prey growth $h(x)$ was chosen to be logistic. The function $p(x)$ is often referred to as the functional response function, describing the effectiveness with which the predator attacks the prey population. Typically, in Gause systems, the numerical response function $q(x)$ is taken to be proportional to the functional response function $p(x)$.

With the choice of $p(x)$ in system (1.1), the system exhibits group defense by the prey population. This means that if the prey population is sufficiently large, they can form a more effective defense against attacks by predators. The particular form of the functional response function $x e^{-\beta x}$ is similar to the well-known Ivlev functional response in the absence of group defense $p(x) = 1-e^{-\beta x}$; see, for example, ^[11]. The essential difference is that the traditional Ivlev functional response function is monotonically increasing in $x$, whereas the group defense functional response in (1.1) has a global maximum for some positive $x$.

From a biological point of view, one of the first studies of this kind of group defense appeared in the book by Tener ^[24]. There, a single muskox is easily eaten by wolves, but when the number of muskoxen increases, it becomes less likely to be attacked by wolves. Another fundamental paper in this area is by Holmes and Bethel ^[8], where it was discussed that parasites can influence host behavior by changing evolutionary strategies to achieve group defense. The relationship with mathematical modeling was made in the fundamental paper by Freedman and Wolkowicz ^[3]. Recent research on the topic of group defense can be found in ^[2,9,21].

From a mathematical modeling point of view, one of the first papers, to our knowledge, mentioning this effect was Wolkowicz's paper ^[25]. Her model is similar to (1.1) but with a different functional response function, $p(x) = \frac{mx}{ax^{2}+bx+1}$, where $m$, $a$, and $b$ are positive constants satisfying $p(0) = 0$, $p^{'}(0) > 0$, $p^{''}(0) < 0$, and $p(x) > 0$ for $x > 0$. For small $x$, this function approximates the Holling Type II functional response.

Mathematically, group defense implies that the predators will be less effective in their attack when there are more prey. This leads to some consequences on the behavior of the functional response function $p(x)$. The following definition was introduced by ^[3].

Definition 1.2. In Gause systems (1.2), the prey population exhibits group defense if $p(x) > 0$, $p'(x) > 0$ ($< 0$) for $0 < x < M$ ($x > M$).

Figure 1.1 shows the function $p(x)$ in the case of group defense.

The main problem in the study of these Gause systems is the determination of the number and stability of limit cycles. This is related to a more general problem for planar autonomous ordinary differential equations, which asks for the number of limit cycles for systems of the following form:

where $x\in \mathbb{R}$, $y\in \mathbb{R}$, and $t\in \mathbb{R}$. The system typically depends on real parameters, indicated by $\vec{\mu}$.

In particular, the unsolved second part of the $16^{th}$ Hilbert problem asks for an upper bound on the number of limit cycles for systems, where $P(x, y)$ and $Q(x, y)$ are polynomial functions of $x$ and $y$; see ^[1,29].

There are few results with the current technology for the study of systems that do not involve small parameters. Specifically, obtaining global results on the qualitative behavior of solutions for all parameters is a notoriously difficult problem. The most successful results have been obtained in the study of the non-existence, existence, and uniqueness of limit cycles, in particular the study of systems where at most one limit cycle exists. There are few studies involving systems with at most two limit cycles. In those few cases, there are restrictive conditions on the functions $P(x, y)$ and $Q(x, y)$, such as satisfying certain symmetry conditions; see the books ^[29,33]. For the solution of the $16^{th}$ Hilbert problem, the development of tools to determine upper bounds for the number of limit cycles for systems with at most two limit cycles would be a first natural but difficult step.

For systems with at most one limit cycle, the strongest methods have been developed for a special case of (1.3), so-called generalized Liénard systems.

More on such systems can be found in the books ^[29,33].

To apply one of the many theorems about limit cycles available for Liénard systems to (1.3), a nonlinear transformation of the variables needs to be found, changing the system into the form (1.4). This is not always possible, but for the system discussed in this paper, (1.1), it is relatively straightforward. For a more general discussion, see ^[17].

Additionally, we will compare system (1.1) with similar group defense mechanisms in other predator-prey models. One instance is when the predator growth dynamics is of Leslie-Gower type, while keeping the same dynamics for the prey. A second instance is to take a similar Gause model, but with a slightly different functional response function $p(x)$, in particular the well-known case of Wolkowicz ^[25] and Rothe-Shafer ^[22]. In both systems, we will show that typically situations occur with two limit cycles and that the uniqueness of limit cycles in system (1.1) is a rather rare occurrence in systems with group defense. Although it is difficult to formulate a general statement, we will indicate the mathematical reasons for why such differences occur.

In Section 2, we will outline the use of Liénard systems and introduce the uniqueness theorem that will be the main tool in the study of (1.1). In Section 3, we will transform system (1.1) into Liénard form and apply the uniqueness theorem to prove Theorem 1.1. In Section 4, we show that a similar system with group defense and Leslie-Gower predator dynamics can have a weak focus of order two and that two limit cycles can be created in a degenerate Hopf bifurcation. Finally, we discuss the mathematical reasons why (1.1) only has at most one limit cycle, while an apparently similar Gause model introduced by ^[25] and further studied by ^[22,27] can have two limit cycles.

2.   Prerequisites

The use of Liénard systems (1.4) in the theory of limit cycles has proved to be advantageous. They contain some of the strongest available theorems about the existence or uniqueness of limit cycles.

In practice, we can try to convert a planar system of nonlinear autonomous ordinary differential equations (1.3) into a Liénard system (1.4) and then analyze it.

To prove that a certain Liénard system has at most one limit cycle, we will adapt two theorems for our purposes. The first is a simple theorem introduced in ^[12] (Theorem 1.1) to show that for some cases, no limit cycles can occur:

Theorem 2.1. (^[12]) Suppose system (1.4) with $F(x) = \int_{x_g}^x{f(\bar{x})d\bar{x}} \in C^3$, $g(x) \in C^2$, $f(x) \in C^2$, satisfies the following conditions on the interval $r_1 < x < r_2$:

(i) $\frac{d\phi(y)}{dy} > 0$,

(ii) $x_{g}\in(r_{1}, r_{2})$ is the unique value such that $(x-x_{g})g(x) > 0$, for $x\neq x_{g}$, and $g(x_{g}) = 0$,

(iii) $x_{f}\in(r_{1}, r_{2})$ is the unique value such that $(x-x_{f})f(x) < 0$, for $x\neq x_{f}$, and $f(x_{f}) = 0$,

(iv) $\exists c\in \mathbb{R}$, such that $f(x)-cg(x)$ has no zeroes in $r_{1} < x < r_{2}$.

Then, (1.4) does not have limit cycles in the strip $r_1 < x < r_2$.

There is a graphical interpretation of condition (iv) of Theorem (2.1). It implies that the function $\frac{f(x)}{g(x)}$ has a gap. It means that there is a horizontal line $L$: $y = c$ such that part of the graph of $y = \frac{f(x)}{g(x)}$ lies above $L$ and part of it below $L$, but without any intersections. See Figure 2.1.

Condition (i) of Theorem 2.1 is a condition typically satisfied in applications. In Gause predator-prey systems, $\phi(y) = e^y$ and the condition is obviously satisfied.

Condition (ii) of Theorem 2.1 means that only limit cycles surrounding one singularity are considered. In Gause systems, this is typically true, but for exceptions (for example, in cases where an Allee effect is introduced in the natural growth of the prey), see ^[6,16].

Condition (iii) of Theorem 2.1 implies that the divergence of system (1.4) changes sign for exactly one value of $x$. This condition can be relaxed, but for this paper it is more convenient to use this form. The differentiability conditions on the three functions $F(x)$, $g(x)$, $f(x)$ and the sign of the condition $(x-x_{f})f(x) < 0$ can be relaxed as well, but were chosen to be consistent with the system studied in this paper.

The non-existence Theorem (2.1), in practice, can be combined with a uniqueness theorem for limit cycles. For our purposes, the following uniqueness theorem will be sufficient. It is a special case of a more general theorem presented in ^[34] (Theorem 3, page 485).

Theorem 2.2. (^[34]) Suppose system (1.4) with $F(x) = \int_{x_g}^x{f(\bar{x})d\bar{x}} \in C^3$, $g(x) \in C^2$, $f(x) \in C^2$, satisfies the following conditions on the interval $r_1 < x < r_2$:

(i) $\frac{d\phi(y)}{dy} > 0$,

(ii) $x_{g}\in(r_{1}, r_{2})$ is the unique value such that $(x-x_{g})g(x) > 0$, for $x\neq x_{g}$, and $g(x_{g}) = 0$,

(iii) $x_{f}\in(r_{1}, r_{2})$ is the unique value such that $(x-x_{f})f(x) < 0$, for $x\neq x_{f}$, and $f(x_{f}) = 0$, and $x_g < x_f$,

(iv) $\frac{d\frac{f(x)}{g(x)}}{dx} < 0$ in $r_{1} < x < x_g$ and $x_f < x < r_{2}$.

Then, (1.4) has at most one limit cycle in the strip $r_1 < x < r_2$, which is stable and hyperbolic if it exists.

Theorems 2.1 and 2.2 are in fact complimentary and can be combined in one theorem, which is more natural for practical applications. Conditions (i)–(iii) of Theorem 2.2 are the same as the conditions in Theorem 2.1. The difference lies in Condition (iv). As mentioned above, in Theorem 2.1, there is a geometrical gap interpretation of the condition. A similar approach is possible for Condition (iv) in Theorem 2.2.

Lemma 2.3. A sufficient condition for Condition (iv) of Theorem (2.2) to hold true is: $\forall c \in \mathbb{R}$, $f(x)-cg(x)$ has no multiple zeroes on the interval $r_{1} < x < x_g$ and $x_f < x < r_{2}$.

Proof. Suppose the function $(\frac{f(\bar{x})}{g(\bar{x})})^{'}$ has a zero $\bar{x}$, then we can write:

Denote the function value of $\frac{f(\bar{x})}{g(\bar{x})}$ at $x = \bar{x}$ by $c^{*}$, i.e.,

These equations imply that

That means $(\frac{f(\bar{x})}{g(\bar{x})})^{'}$ has a multiple zero, which is not possible. Therefore, $(\frac{f(\bar{x})}{g(\bar{x})})^{'}$ does not have a zero and, from Conditions (ii) and, (iii) of Theorem 2.2, it follows that it has to be negative. □

This simple lemma has a geometrical interpretation, as shown in Figure 2.2: If a differentiable function has a local minimum or maximum at $x = \bar{x}$, then necessarily the horizontal line $y = \frac{f(\bar{x})}{g(\bar{x})}$ corresponds to a value $c$ where $f(x)-c g(x)$ has a multiple zero, because the line is tangent to the graph of the function.

Conditions (iv) of Theorems 2.1 and 2.2 are complimentary, because of the geometrical interpretation in Lemma 2.3. We combine both theorems into one convenient theorem. This leads to a similar theorem as Lemma $1.4$ in ^[12]. The only difference with that lemma is the interval for $x$ in Condition (iv). In Lemma $1.4$ ^[12], the condition needs to hold true for the whole interval $r_1 < x < r_2$, while in Theorem (2.2), the interval is restricted. This subtle difference will simplify the proof for the uniqueness of limit cycles in (1.1).

Theorem 2.4. Suppose system (1.4) with $F(x) = \int_{x_g}^x{f(\bar{x})d\bar{x}} \in C^3$, $g(x) \in C^2$, $f(x) \in C^2$ satisfies the following conditions on the interval $r_1 < x < r_2$, $y \in \mathbb{R}$:

(i) $\frac{d\phi(y)}{dy} > 0$,

(ii) $x_{g}\in(r_{1}, r_{2})$ is the unique value such that $(x-x_{g})g(x) > 0$, for $x\neq x_{g}$, and $g(x_{g}) = 0$,

(iii) $x_{f}\in(r_{1}, r_{2})$ is the unique value such that $(x-x_{f})f(x) < 0$, for $x\neq x_{f}$, and $f(x_{f}) = 0$, and $x_g < x_f$.

(iv) Either $(iv)^{'}$: $\exists c\in \mathbb{R}$, such that $f(x)-cg(x)$ has no zeroes in $r_{1} < x < r_{2}$,

or $(iv)^{''}$: $\forall c \in \mathbb{R}$, $f(x)-cg(x)$ has no multiple zeroes on the interval $r_{1} < x < x_g$ and $x_f < x < r_{2}$.

Then, (1.4) has at most one limit cycle in the strip $r_{1} < x < r_{2}$, and if it exists, it is hyperbolic. In case $(iv)^{'}$, no limit cycles occur.

Proof. The theorem is a combination of Theorems 2.1 and 2.2, which have the same Conditions (i)–(iii). Condition $(iv)^{'}$ is Condition (iv) from Theorem 2.1 and Condition $(iv)^{''}$ is Condition (iv) from Theorem 2.2 rewritten in terms of Lemma 2.3. □

Conditions (i)–(iii) can easily be verified in our case. The critical conditions are $(iv)^{'}$ and $(iv)^{''}$. These are, in principle, two independent conditions with different implications, but in typical families (the system in this paper is an example of this), their combination is sufficient to cover all parameter values in a system.

2.1.   Singularities

Compared to the general system (1.3), Liénard systems have another advantage. It is straightforward to determine the position of the singularities, because they satisfy an equation with one variable only.

Lemma 2.5. In system (1.4) with $\frac{d\phi(y)}{dy} > 0$, a necessary and sufficient condition for the existence of a singularity at $(x_g, y_g)$ is $g(x_g) = 0$, and $\phi^{(-1)}[F(x_g)]$ exists. If $g'(x_g) > 0$ ($< 0$), then the singularity is an anti-saddle (saddle).

Proof. Singularities of the system satisfy $F(x_g)-\phi(y_g) = 0$, $g(x_g) = 0$. From this, the necessary and sufficient conditions for the coordinates of a singularity follow.

The Jacobian matrix at a singularity, where $g(x_g) = 0$, $F(x_g)-\phi(y_g) = 0$, is given by $\begin{bmatrix} F'(x_g) & -\frac{d\phi(y_g)}{dy} \\ g'(x_g) & 0 \end{bmatrix}$. The product of the eigenvalues is $\frac{d\phi(y_g)}{dy} g'(x_g)$, from which the result of the lemma follows. □

The strategy for the analysis of system (1.1) is to convert it into (1.4) first, then to study the singularities using Lemma 2.5, and finally, to use Theorem 2.4 to establish the uniqueness of limit cycles.

3.   Gause model with group defense

To apply the results of the previous section to system (1.1), the system first needs to be converted into Liénard form. This is a well-known procedure for Gause predator-prey systems ^[4]. In general, different transformations into Liénard form exist, a concept not widely understood yet. The choice of transformation is not unique, and different choices will lead to different results. For examples of how to convert predator-prey systems to a Liénard system, see ^[4,17,28,30]. In this paper, we choose the simplest transformation with the justification that it is sufficient for getting proof of the uniqueness of limit cycles.

Lemma 3.1. System (1.1) can be transformed into a Liénard system (1.4) with $F(x) = re^{\beta x}(1-\frac{x}{k})$, $g(x) = -D\frac{e^{\beta x}}{x}+\mu$, $\phi(y) = e^{y}$.

Proof. Let $t \rightarrow \frac{t}{p(x)}$, with $p(x) = x e^{-\beta x }$, and $y = e^{v}$. After the transformation, we label for convenience $v$ again by $y$ and get:

□

For future reference, we define the divergence of the vector field by:

3.1.   Singularities

To find the position of the singularities in the first quadrant of the phase plane (since $x(t) \geq0$ and $y(t) \geq 0$), we apply Lemma 2.5 to the Liénard system (3.1) in Lemma 3.1.

Lemma 3.2. In system (1.1) define $\mu^{*} = D\beta e$ and $\tilde{\mu} \equiv D\frac{e^{\beta k}}{k}$, with $0 < \mu^{*} < \tilde{\mu}$.

● If $\mu^{*} < \mu$, then the equivalent Liénard system (3.1) has an anti-saddle with $x$-coordinate $x_{AS}$ and if $\mu^{*} < \mu < \tilde{\mu}$, it has a saddle at $x = x_{S}$. For the $x$-coordinates of the singularities, the inequality $0 < x_{AS} < \frac{1}{\beta} < x_{S} < k$ holds, if the saddle exists.

● If $\mu < \mu^{*}$, then system (3.1) has no singularities.

Proof. The necessary and sufficient conditions for the existence of a singularity in Liénard system (3.1), according to Lemma 2.5, translate into the following conditions for the system in Lemma 3.1.

The number of solutions to the equation $g(x_g) = 0$ follows from looking at the derivative $g'(x) = \frac{De^{\beta x}(1-\beta x)}{x^{2}}$. A double zero can only occur when $x_g = \frac{1}{\beta}$. Substituting this into $g(x_g)$, we find the condition for a double zero to be $\mu = \mu^{*} = D\beta e$. If $\mu > \mu^{*}$, then $g(x_g) = 0$ has two solutions $x_g > 0$, and for $\mu < \mu^{*}$, there are no solutions. The Jacobian matrix of system (3.1) is:

Under the condition $\mu > \mu^{*}$, we find that for the two zeroes $x_{AS}$ and $x_{S}$ of $g(x)$, with $x_{AS} < \frac{1}{\beta} < x_{S}$ the determinant is positive and negative, respectively, i.e., $x_{AS}$ is an anti-saddle and $x_{S}$ is a saddle. Moreover, $0 < x_{AS} < \frac{1}{\beta} < x_{S} < k$.

The last inequality in (3.3) is the condition that $y_g = ln(F(x_g) = ln(re^{\beta x_g}(1-\frac{x_g}{k}))$ exists. Since $x_g$ satisfies $g(x_g) = 0$, we get that $0 < x_g < k$ implies that $\mu < \tilde{\mu} \equiv D\frac{e^{\beta k}}{k}$. For the two critical values of $\mu$, the inequality $0 < \mu^{*} < \tilde{\mu}$ holds true, because $k \beta > 1$.

The condition $\mu < \tilde{\mu}$ reflects the fact that the saddle in the original system (1.1) is located in the first quadrant. For $\mu > \tilde{\mu}$, the saddle moves into the region $y < 0$, which is irrelevant from a biological point of view. Since the Liénard system was obtained through a transformation $y = e^v$, only singularities in the first quadrant of the original system will appear in the Liénard form. □

Since limit cycles cannot appear in a system without singularities, we will impose in the following the condition $\mu > \mu^{*}$ of Lemma 3.2. It means that implicitly we will assume that $\mu > D \beta e$. The conclusions of Lemma 3.1 are in line with the results of ^[26]. In principle, there are two cases to be considered according to the lemma: $\mu < \tilde{\mu}$ (saddle exists) and $\mu > \tilde{\mu}$ (saddle does not exist). However, for the proof of the uniqueness of limit cycles, this distinction will not be relevant. In both cases, the condition $\mu > \mu^{*}$ ensures that the equation $g(x_g) = 0$ has two solutions denoted by $x_{AS}$ and $x_S$. We will prove the uniqueness of limit cycle in the extended interval $0 < x < x_{S}$ even if that zero does not represent a saddle.

The stability of the anti-saddle is determined by the sign of $f(x)$ for $x = x_{AS}$. If it is negative (positive), then the singularity is stable (unstable). In the case when $f(x_{AS}) = 0$, the sign of the first focal value will determine the stability. For a Liénard system (1.4), there exists a simple method to determine it. The unique zero of $f(x)$ is given by $x_{f} = k-\frac{1}{\beta} < k$.

Lemma 3.3. Under the condition $x_{f} = x_{AS}$, system (1.1) has a weak focus, which is stable and of order one. In terms of the parameters of the original system, this occurs when $\mu = \mu_{wk} \equiv D\beta\frac{e^{k\beta-1}}{k\beta-1}$.

Proof. According to system (3.1) and Eq (3.2), we have $f(x) = re^{\beta x}[\beta(1-\frac{x}{k})-\frac{1}{k}]$, $g(x) = \frac{\mu x-De^{\beta x}}{x}$. Then $f(x_{f}) = 0$, with $x_{f} = k-\frac{1}{\beta}$. Therefore, a weak focus occurs for $g(x_{f}) = 0$, which is equivalent to $\mu = \mu_{wk} \equiv D\beta\frac{e^{k\beta-1}}{k\beta-1}$.

Under this condition, we have $x_f = x_{AS}$ and the derivatives of the functions are given by $f^{'}(x) = r\beta e^{\beta x}[\beta(1-\frac{x}{k})-\frac{2}{k}]$, $f^{''}(x) = r\beta e^{\beta x}[\beta(1-\frac{x}{k})-\frac{3}{k}]$, which evaluated at $x_{f}$ leads to $f^{'}(x_{f}) = r\beta e^{\beta x_{f}}[\beta(1-\frac{x_{f}}{k})-\frac{2}{k}] = r\beta e^{\beta x_{f}}(-\frac{1}{k}) < 0$, $f^{''}(x_{f}) = r\beta e^{\beta x_{f}}[\beta(1-\frac{x_{f}}{k})-\frac{3}{k}] = r\beta e^{\beta x_{f}}(-\frac{2}{k}) < 0$.

Similarly, we get $g^{'}(x_{f}) = \frac{e^{\beta x_{f}}(1-\beta x_{f})}{x_{f}^{2}}$, $g^{''}(x_{f}) = \frac{e^{\beta x_{f}}}{x_{f}^{3}}(-\beta^{2}x_{f}^{2}+2\beta x_{f}-2)$.

From Lemma 3.2, we know that $0 < x_{AS} = x_f < \frac{1}{\beta}$, so $g^{'}(x_{f}) = g^{'}(x_{AS}) = \frac{e^{\beta x_{AS}}(1-\beta x_{AS})}{x_{AS}^{2}} > 0$, and $g^{''}(x_{f}) = g^{''}(x_{AS}) = \frac{e^{\beta x_{AS}}}{x_{AS}^{3}}(-\beta^{2}x_{AS}^{2}+2\beta x_{AS}-2) < g^{''}(\frac{1}{\beta}) = e \beta^3 (-1) < 0$.

Under the condition $x_{f} = x_{AS}$, we get according to ^[31] the following expression for the first focal value for a generalized Liénard system.

This implies that Liénard system (3.1) has a weak focus, which is stable and of order one.

System (1.1) is equivalent to the Liénard system (3.1). Therefore, if the original system (1.1) has a weak focus, it is stable and has order one. □

3.2.   Checking the conditions of Theorem 2.4

To apply Theorem 2.4, we need to verify that all conditions are satisfied on a relevant interval for $x$, i.e., the strip in the phase plane where limit cycles are residing. First, we determine the interval $r_{1} < x < r_{2}$. We can determine the interval by using the properties of Liénard systems.

Lemma 3.4. Limit cycles in (1.1) are restricted to the strip $0 < x < x_S$ in the phase plane. Here, $x_S$ is one of the two solutions $x_{AS} < x_{S}$ of $g(x) = 0$. If $\mu^{*} < \mu < \tilde{\mu}$, according to Lemma 3.2, $x_{S}$ represents the $x$-coordinate of the saddle.

Proof. Limit cycles surround the anti-saddle with $x$-coordinate $x_{AS}$ with $0 < x_{AS} < x_{S}$ and cannot cross the vertical line $x = x_S$ in the phase plane. This is true, because if the saddle at $x = x_{S}$ exists, then on the line $\frac{dx}{dt} = F(x_{S})-\phi(y_{S})$ changes sign exactly at the saddle. Any limit cycle crossing the line would have to contain the saddle, which is impossible, because the system only has two singularities, one anti-saddle and one saddle.

If there is no saddle at $x = x_{S}$, then it implies that $k < x_{S}$. In the original system (1.1), the line $x = k$ is a line without contact and cannot be crossed by limit cycles. Therefore, limit cycles cannot cross any other vertical line $x = c > k$. □

It means that for the application of Theorem 2.4, we will effectively use $r_{1} = 0$, $r_{2} = x_{S}$, even if $x_{S}$ does not correspond to a saddle.

Lemma 3.5. Under the Condition $\mu > D\beta e$, Liénard system (3.1) satisfies the first three conditions (i)–(iii) of Theorem 2.4 on the interval $r_{1} < x < r_{2}$, with $r_{1} = 0$, $r_{2} = x_{S}$.

Proof. Since $\phi(y) = e^{y}$, $g(x) = \frac{\mu x-De^{\beta x}}{x}$, and $F(x) = re^{\beta x}(1-\frac{x}{k})$ in system (3.1), and $r_{1} = 0$, $r_{2} = x_{S}$ according to Lemma 3.4, we get $\frac{d\phi(y)}{dy} > 0$ and $f(x) = F'(x) = re^{\beta x}[\beta(1-\frac{x}{k})-\frac{1}{k}]$. Moreover, the zeroes of $g(x)$ satisfy $0 < x_{AS} < \frac{1}{\beta} < x_{S} < k$ according to Lemma 3.2 and the unique zero of $f(x)$ is given by $x_{f} = k-\frac{1}{\beta}\in(0, x_{S})$. If $x_f > x_S$, then the divergence of the system has fixed sign on $0 < x < x_{S}$ and no limit cycles can exist. Therefore, we can safely exclude this case in the following.

If we denote $x_{g} = x_{AS}$, then $x_{g}\in(0, x_{S})$ is the unique value such that $g(x_{g}) = 0$. We know $g'(\frac{1}{\beta}) = 0$ according to Lemma 3.2. As a consequence, $(x-x_{AS})g(x) > 0$ where $x\neq x_{AS}$, i.e., $(x-x_{g})g(x) > 0$ where $x\neq x_{g}$. □

To apply Theorem 2.4, we need to investigate the zeroes of the function $f(x)-c g(x)$ and check whether one of the two Conditions $(iv)^{'}$ or $(iv)^{''}$ holds true. To simplify the analysis, we rewrite the expression in the following way:

where

Evidently, $s(x) < 0$ for $0 < x < x_S$, and to check Condition (iv) we can restrict ourselves to analyzing the zeroes of $\bar{f}(x)-\bar{c}\bar{g}(x)$. The change of variables $t = \beta x$ does not change the nature of the zeroes and simplifies the expression further. We get:

where $\bar{f} = t(t-t_{f})$, $\bar{g} = Cte^{-t}-L$, $C = \frac{\mu}{D\beta}$, $L > 0$. We only need to consider the case $L > 0$, because $\bar{f}(x)-\bar{c}\bar{g}(x)$ can only have a zero on the intervals $r_{1} < t < t_{AS}, t_{f} < t < t_{S}$, where we used the notation $t_{AS} = \beta x_{AS}$ and $t_{S} = \beta x_{S}$. In short, we will study in the following the zeroes of the function $\bar{z} = \bar{f}-L\bar{g}$ for $L > 0$.

Next, we will check Condition $(iv)$ of Theorem 2.4. As we will see, the two subcases $(iv)^{'}$ and $(iv)^{''}$ correspond to the relative positions of $t_{AS}$ and $t_f$.

3.3.   The case $x_{AS} > x_f$

In this case, we expect no limit cycles and construct a constant $c$ such that Condition $(iv)^{'}$ in Theorem 2.4 is satisfied.

Proposition 3.6. Under the condition $\mu > D\beta e$, there exists a constant $c$ such $f(x)-c g(x)$ in Liénard system (3.1) has fixed sign on the interval $0 < x < x_S$.

Proof. Following the discussion of the previous section, we will consider the equivalent expressions for the variable $t = \beta x$. Under the condition $\mu > D\beta e$, we know $t_{AS} > t_{f}$ where $t_{AS} = \beta x_{AS}$, $t_{f} = \beta x_{f}$. We want to choose $L$ in such a way that the tangent lines of $\bar{f}(t)$ at $t_{f}$ and $L\bar{g}(t)$ at $t_{AS}$ have the same slope, i.e., $\bar{f}^{'}(t_{f}) = L\bar{g}^{'}(t_{AS})$. Denote the tangent lines of $\bar{f}$ and $L\bar{g}$ by $l_{f}$ and $l_{AS}$, respectively. We get the situation of Figure 3.1.

From Figure 3.1, we see that the graph of $\bar{f}(t)$ lies above $l_{f}$, and the graph of $L\bar{g}(t)$ lies below $l_{AS}$. This would imply that $\bar{f}(t) = L\bar{g}(t)$ has no solution for this value of $L$.

To prove analytically that the function $f(t)-Lg(t)$ has fixed sign for this choice of $L$, we consider the expression for $l_{f}$ and $l_{AS}$:

where $L^{*} = \frac{\bar{f}^{'}(t_{f})}{\bar{g}^{'}(t_{AS})}$.

Next, let $H_{1}(t) = \bar{f}(t)-y_f = \bar{f}(t)-\bar{f}'(t_{f})(t-t_{f})$. Then, the condition that the graph of $\bar{f}(t)$ lies above $l_{f}$ is equivalent to $H_{1}'(t)\geqslant0$. Taking derivatives we get: $H_{1}''(t) = \bar{f}''(t) > 0$, $H_{1}(t_{f}) = 0$, and $H_{1}'(t_{f}) = 0$. From this, it follows that $H_{1}(t) = \int_{t_{f}}^{t}\int_{t_{f}}^{t_{2}} H^{''}_{1}(t_{1})dt_{1}dt_{2}\geqslant0$.

Similarly, let $H_{2}(t) = y_{AS}-L^{*}\bar{g}(t) = L^{*}(\bar{g}'(t_{AS})(t-t_{AS})-\bar{g}(t))$. The condition that the graph of $L^{*}\bar{g}(t)$ lies below $l_{AS}$ is equivalent to $H_{2}(t)\geqslant0$. For this function, we have $H_{2}(t_{AS}) = H_2'(t_{AS}) = 0$ and $H_{2}''(t) > 0$ for $0 < t < 2$. Taking derivatives leads to $H_{2}(t)\geqslant0$ in $t\in(0, 2)$, and $H_{2}(2) > 0$. It follows that $H_{2}(t) = \int_{t_{AS}}^{t}\int_{t_{AS}}^{t_{2}} H^{''}_{2}(t_{1})dt_{1}dt_{2}\geqslant0$ for $0 < t < 2$.

For $t\geqslant2$, $H_{2}'(t) = L^{*} [\bar{g}'(t_{AS})-\bar{g}'(t)] > 0$, and $H_{2}'(t)$ increases monotonically for $t\geqslant2$. Therefore, $H_{2}(t)\geqslant H_{2}(t)_{min} = H_{2}(2) > 0$. This implies that $H_{2}(t)\geqslant0$ for $t > 0$ as we intended to prove. It follows that the graph of $L^{*} \bar{g}(t)$ lies below $l_{AS}$.

We have proved that $\bar{f}(t)\geqslant\bar{f}'(t_{f})(t-t_{f}) > L^{*}\bar{g}'(t_{AS})(t-t_{AS})\geqslant L^{*}\bar{g}(t)$, which leads to the result we are looking for: There exists an $L = L^{*}$ such that $f(x)-L^{*}g(x)\geqslant 0$. □

3.4.   The case $x_{AS} < x_f$

In this case, we expect at most one limit cycle and will show that Condition $(iv)^{''}$ in Theorem 2.4 is satisfied.

Proposition 3.7. Under the condition $\mu > D\beta e$, $f(x)-c g(x)$, $\forall c \in \mathbb{R}$, in Liénard system (3.1) does not have a multiple zero on the intervals $0 < x < x_{AS}$ and $x_f < x < x_S$.

Proof. Again, we will consider the equivalent expressions for the variable $t = \beta x$. First, we prove $\nexists t^{*}\in(0, t_{AS})$ such that $\bar{z}(t^{*}) = \bar{z}^{'}(t^{*}) = 0$.

According to (3.6), we know $\bar{f}^{''}(t) = 2$, $\bar{g}^{''}(t)(t-2) < 0$ for $t\in(0, 2)$ and $t_{AS} < 2$, so $\bar{z}^{''}(t) = \bar{f}^{''}(t)-L\bar{g}^{''}(t) > 0$, for $t\in(0, 2)$.

Suppose $\exists t^{*}$, such that $\bar{z}(t^{*}) = \bar{z}^{'}(t^{*}) = 0$. Since $\bar{z}^{''}(t) > 0$ for $t\in(0, 2)$, we know that $\bar{z}^{'}(t)(t-t^{*})\geqslant0$, i.e., $\bar{z}^{'}(t) < 0(>0)$, for $0 < t < t^{*}(t^{*} < t < t_{AS} < 2)$. It follows that

However, $\bar{z}(t_{AS}) = \bar{f}(t_{AS})-L\bar{g}(t_{AS}) = \bar{f}(t_{AS}) < 0$. This contradicts the assumption that $\bar{z}(t^{*}) = \bar{z}^{'}(t^{*}) = 0$. Therefore, $\nexists t^{*}\in(0, t_{AS})$ such that $\bar{z}(t^{*}) = \bar{z}^{'}(t^{*}) = 0$. The above calculation shows that the result is true for the extended interval $(0, 2)$, i.e., $\nexists t^{*}\in(0, 2)$ such that $\bar{z}(t^{*}) = \bar{z}^{'}(t^{*}) = 0$.

Finally, we prove that under the condition $\mu > D\beta e$, $\bar{z} = \bar{f}-L\bar{g}$ does not have a multiple zero on the interval $2 < t < t_{S}$.

According to (3.6), $\bar{z}^{'}(t) = \bar{f}^{'}(t)-L\bar{g}^{'}(t)$ and $\bar{f}^{'}(t) > 0$, $\bar{g}^{'}(t) < 0$ in $(2, t_{S})$. It follows that $\bar{z}^{'}(t) > 0$ in $(2, t_{S})$ and there cannot be a $t^{*}$ such that $\bar{z}(t^{*}) = \bar{z}^{'}(t^{*}) = 0$, for $t^{*}\in(2, t_{S})$. Therefore, $y = \bar{f}(t)$ and $y = L\bar{g}(t)$ have exactly one transversal intersection in $(2, t_{S})$.

Combining the results, we find that under the condition $\mu > D\beta e$, $f(x)-c g(x)$ in Liénard system (3.1) does not have a multiple zero on the intervals $0 < x < x_{AS}$ and $x_f < x < x_S$. The images of $\bar{f}$ and $L\bar{g}$ are shown in Figure 3.2.

□

3.5.   Uniqueness of limit cycle

With the results of the previous sections, we have sufficient information to prove the main result of this paper.

Theorem 3.8. System (1.1) has at most one limit cycle. If it exists, it is stable and hyperbolic.

Proof. According to Lemma 3.5, Propositions 3.6 and 3.7, we know that Liénard system (3.1) satisfies all conditions of Theorem (2.4). Therefore, Liénard system (3.1) has at most one limit cycle.

In particular, Liénard system (3.1) has no limit cycles (at most one limit cycle) when $\mu\leq\mu_{WF}$ ($\mu > \mu_{WF})$ if the Condition $(iv)^{'}$ (Condition $(iv)^{''}$) in (2.4) is satisfied according to Proposition (3.6) (Proposition (3.7)). Here, $\mu_{WF}$ is defined in Lemma 3.3.

System (1.1) is equivalent to the Liénard system (3.1) and the uniqueness of limit cycle in the Liénard system (3.1) implies the uniqueness of the limit cycle in the original system (1.1), i.e., system (1.1) has at most one limit cycle.

In the case when $\mu > \mu_{WF}$, the divergence of the vector field $f(x)$ at $x = x_{g}$ is positive in the case where a limit cycle can exist ($x_{AS} < x_{f}$). For system (1.1), this means that the limit cycle, if it exists, must be stable, because the singularity inside the limit cycle is unstable. □

3.6.   Existence of limit cycle

Using Lemmas 3.2 and 3.3, we establish the existence of a limit cycle by a suitable perturbation of the weak focus leading to a small-amplitude limit cycle created in a Hopf bifurcation. A further change of the parameters will lead to the disappearance of the stable limit cycle in a stable saddle loop. Figure 3.3 shows a numerical scenario where the parameter $k$ is changed while the other parameters are fixed.

3.7.   Biological interpretation of the results

Even though it is difficult to show a clear pattern in the asymptotic behaviour of the solutions for all values of the parameters, the numerical example of Section 3.6 shows what can typically be expected.

In the case when $x_{AS} > x_f$, system (1.1) does not have limit cycles, as indicated in the proof of Theorem 3.8. A typical situation is depicted in Figure 3.3(a). The anti-saddle is stable and locally attracts solutions.

It means that in the positive quadrant of the phase plane, for some initial values, the solution $[x(t), y(t)]$ will tend to a positive equilibrium, while for some initial values it will tend to a stable node situated at $x = k, y = 0$, i.e., leading to extinction of the predators. The regions in the phase plane leading to extinction of the predators or a positive equilibrium are determined by the relative position of the separatrix of the saddle lying in the first quadrant. The exact positions can only be found numerically.

In the case when $x_{AS} < x_f$, system (1.1) has an unstable positive singularity. For increasing $k$, see Figure 3.3(b)–(d), a growing stable limit cycle surrounds the equilibrium. In those situations, depending on the position in the phase plane of the separatrices of the saddle and the position of the limit cycle, the solution will either tend to the stable limit cycle or to the stable node situated at $x = k, y = 0$, i.e., the predator and prey populations will approach a periodic limit set or the predator population will die out.

When $k$ is large enough, the limit cycle has disappeared into a saddle loop and the ensuing system has a positive unstable anti-saddle, no limit cycles, a saddle, and a stable node situated at $x = k, y = 0$; see Figure 3.3(e). In that case, regardless of the initial position, the predator population will die out. The biological interpretation is that if the carrying capacity is large, the prey population can grow big enough to contain the predation effect through group defense, and the predator population will die.

4.   Connections to other predator-prey systems with group defense

A natural question is how typical the behavior of system (1.1) is, compared to other systems with group defense.

The first system with group defense that was studied systematically in the literature ^[22,25,27] contained cases with two limit cycles. A first question therefore would be: Why does (1.1) only have one limit cycle? In ^[13], the group defense property of the functional response function $p(x)$ in a Gause system was discussed. It was observed that actually a weak focus of second order can only be created for those values of $x$ where $p(x)$ is increasing, implying that the group defense property is not the essential reason to why more than one limit cycle can occur.

To study this difference in the number of limit cycles between the various predator-prey models with group defense, we look at two systems that have a similar structure to (1.1).

(modified Leslie-Gower with group defense)

(Wolkowicz-Rothe-Shafer with group defense)

Both systems bear a similarity with (1.1) in the following sense:

System (4.1) has the same dynamics for the prey, i.e., the same logistic growth and same functional response function, but the predator dynamics contains a version of the carrying capacity term introduced by Leslie-Gower ^[14]. It has not occurred in the literature as far as we know. Its purpose here is purely mathematical: If the predator dynamics in system (1.1) are changed into a Leslie-Gower form while keeping the same dynamics for the prey population including group defense of Ivlev-type, what will be the consequences, especially in terms of the number of limit cycles?

System (4.2) is the system studied by Wolkowicz-Rothe-Shafer ^[22,25], while ^[27] showed that a second-order weak focus can occur in the system. The only difference with (1.1) is the denominator in the functional response function $p(x) = \frac{x}{M(x)}$: $M_{XR}(x) = e^{\beta x}$ in ^[26] versus $M_{WRS}(x) = x^2+a x +b$ in ^[22,25]. These two functions have similar behavior for $x > 0$: $M(x) > 0$ (the parameters in (4.2) satisfy conditions to ensure this; see ^[22]), and second- and higher-order derivatives cannot be negative. One difference is that the first derivative of $M_{XR}(x)$ is positive, while the first derivative of $M_{WRS}(x)$ can be negative (if $a < 0$). We note that for system (4.2), the conjecture is that for all parameters, at most two nested limit cycles can occur ^[22].

4.1.   Comparison to a modified Leslie-Gower system

Proposition 4.1. System (4.1) can have two small-amplitude limit cycles created in an Andronov-Hopf bifurcation.

Proof. First, we transform system (4.1) into a Liénard system. This is more complicated than for a standard Gause system because of the $y^2$-term in $\frac{dy}{dt}$. Following ^[30], we apply the transformation $y = \omega(x) u$, with $\omega(x) = e^{\int_{x_0}^x{\frac{\mu }{(\lambda+\bar{x}) \bar{x} e^{-\beta \bar{x}}}d\bar{x}}}$, for some constant $x_0 > 0$, and a rescaling of time $t \rightarrow \frac{t}{ xe^{-\beta x} \omega(x)}$, which leads to the Liénard system

with

with

The expressions are quite involved and their full analysis is beyond the scope of this article. To prove the existence of a second-order weak focus, we first impose the condition of a weak focus, i.e., the parameter values for which $f(x)$ and $g(x)$ have a zero in common. To construct such a situation, we introduce a new parameter $x_f$ corresponding to a zero of $f(x)$. Solving $f(x_f) = 0$ for the parameter $r$, we get:

In the same spirit, we can introduce a new parameter $x_g$ corresponding to a zero of $g(x)$ by removing the parameter $\mu$:

Next, we create a weak focus by setting $x_f = x_g$ and $\lambda = 0$ to simplify the expressions. The quotient $\frac{f(x)}{g(x)}$ becomes:

This function is useful, because it determines the focal values of the weak focus at $x = x_g = x_f$; see the proof of Lemma 3.3, Eq (3.4). We get, according to (3.4), the following expression determining the first focal value:

To create an example of a second-order weak focus, we take $\beta = \frac{1}{4}$, $x_g = 1$, $k = \frac{175 +15 \sqrt{137}}{50}$. With this choice, the first focal value $V_1$ vanishes. Using Maple, it can be shown that the second focal value is negative and therefore the original system has a weak focus of order two. Since the system can have a weak focus of order one, which is either stable or unstable ($V_1$ can be made positive or negative by choosing appropriate values of the parameters in the system), we can create a stable small-amplitude limit cycle surrounding an unstable weak focus of order one, by changing the sign of $V_1$. Then, changing the parameters in such a way that the focus becomes a stable strong focus, a second, unstable limit cycle is created inside the already existing stable limit cycle. We note that this particular system has a weak focus of at most second order and the conjecture is that for all parameters at most two nested limit cycles can occur.

In Figure 4.1, a numerical example is shown of system (4.1) with two nested limit cycles.

□

4.2.   Comparison to the Wolkowicz-Rothe-Shafer system

To compare the two systems (1.1) and (4.2), we transform to a Liénard system and compare the functions $F(x)$ and $g(x)$.

In Lemma 3.1, we obtained the Liénard system (3.1) of system (1.1). Next, we transform system (4.2) to a Liénard system in a similar way.

Let $t = \frac{t(x^{2}+ax+b)}{x}$, $y = e^{y}$. Then, system (4.1) becomes a Liénard system (1.4).

In Liénard system (3.1), we obtained $F(x) = re^{\beta x}(1-\frac{x}{k})$. Then, $f(x) = F^{'}(x) = re^{\beta x}[\beta(1-\frac{x}{k})-\frac{1}{k}]$. Let $f(x_{f}) = 0$, which leads to $x_{f} = k-\frac{1}{\beta}$. It follows that $\exists\delta_{1} > 0$, such that $F(x)$ has one local maximum $F(x_{f})$ in $(x_{f}-\delta_{1}, x_{f}+\delta_{1})$.

Similarly, in Liénard system (4.7), we have $F(x) = r(x^{2}+ax+b)(1-\frac{x}{k})$. Then, $f(x) = F^{'}(x) = \frac{r}{k}[-3x^{2}+2(k-a)x+ak-b]$. Setting $f(x) = 0$, we get two positive solutions $x_{1}$, $x_{2}$. It means that in the case of the Wolkowicz-Rothe-Shafer system, for the function $F(x)$, the situation of simultaneously one local minimum and one local maximum can occur.

5.   Conclusions

This paper showed that the conjecture in ^[26] regarding the uniqueness of limit cycles in a Gause predator-prey model with a specific functional response function for group defense is true. The method to achieve this was an adaptation of a standard uniqueness theorem in Liénard systems.

By comparing the system to similar systems, it was shown that it is difficult in general to relate the number of limit cycles to the specific form of the group defense mechanism. In relation to these observations, we make some suggestions for future research.

The functional response functions to model group defense in systems (1.1) and (4.2) had very similar properties. Since their form does not seem to have a specific justification from a biological point of view, the choice of the function in a practical situation would be based on parameter fitting to field data. However, it has been a notoriously difficult problem to fit functional response functions to field data. For example, the four Holling functional response functions, originally introduced in ^[7], were a poor fit to field data, but nevertheless have become a standard in predator-prey modeling. Similarly, in more recent biological studies ^[10,35], it was shown that obtaining good fits with popular functional response functions is not easy.

The functions $p(x) = xe^{-\beta x}$ and $p(x) = \frac{x}{x^2+a x +b}$ used in (1.1) and (4.2), respectively, have similar mathematical properties and are virtually indistinguishable when used for data fitting. This leads to the following issues, that need to be resolved in future research.

● (Mathematical) What causes the difference between systems with almost similar functional response functions, while modelling group defense? The functional response functions $p(x) = xe^{-\beta x}$ and $p(x) = \frac{x}{x^2+a x +b}$ have similar properties, but in the Liénard form, the functions $F(x)$ have one local maximum and one local maximum + one local minimum, respectively, as shown in Section 4.2. It is a worthwhile future research topic to investigate what the real mathematical factor is determining the maximum number of limit cycles.

● (Biological) It is difficult to obtain a good fit with field data for the functional response function. This could either be circumvented by finding a biological interpretation of the parameters in a suitable functional response function or a detailed statistical analysis of field data. This seems to be a general problem that needs to be resolved in many predator-prey models introduced in the literature.

● (Mathematical) It is not clear if some systems can have more than two limit cycles. Several recent papers contain examples with more than two limit cycles ^[20,23,31]. Many other systems seem to have at most two limit cycles, see ^[18,19] amongst many others. In general, no methods exist to get global upper bounds for the number of limit cycles if the number is two or more. Any progress in this direction for predator-prey systems would be of importance for the still open second part of the $16^{th}$ Hilbert problem ^[1,29,33], which asks for upper bounds for the number of limit cycles in polynomial systems.

● (Mathematical) The parameters in a predator-prey system typically do not correspond to a rotated vector field (see Chapter IV, section 3 in ^[33] or Chapter 3 in ^[29]). This means it is difficult to trace limit cycles and saddle loops as a function of the parameter. From a biological point of view, it makes it difficult to indicate in the positive quadrant $x \geq 0, y \geq 0$ which solutions stay positive and which solutions lead to extinction of either the predator or prey population. An understanding of the role of the parameters in the Gause system (1.2) or specifically in (1.1) is needed. For example, in Figure 3.3, a clear numerical pattern of the phase portraits as a function of the parameter $k$, i.e., the carrying capacity of the growth rate of the prey population in the absence of predators is observable. However, a rigorous mathematical proof of this pattern is still needed.

Author contributions

Jin Liao: Formal analysis, writing-original draft; André Zegeling: Validation, writing-review and editing; Wentao Huang: Validation, review. All authors have read and approved the final version of the manuscript for publication.

Acknowledgments

This work is partially supported by Innovation Project of Guangxi Graduate Education(YCSW2024152) and by the National Natural Science Foundation of China (No.12461099).

Conflict of interest

The authors declare that there is no conflict of interest.