Stability analysis of solutions of certain May's host-parasitoid model by using KAM theory

1.   Introduction and preliminaries

In this paper, we investigate the following general May's host-parasitoid model:

where $f$ is a sufficiently smooth and strictly decreasing function such that $f:\left[0, +\infty \right) \rightarrow \left(0, +\infty \right)$, $f\left(0\right) > 0$, $f$ has only one positive root of the equation $bxf(x)-a = 0$, and $x_{-1}, x_{0}\in \left[0, +\infty \right)$ are the initial conditions. In this model, $x_n$ represents the host density and $y_n$ represents the parasitoid density at the $n$th generation. Obviously, the solutions of the system are positive for all initial conditions from $\mathbb{R}^{2}_+$.

Model (1.1) is a special case of the following general model that describes the host-parasitoid behavior in discrete time:

where $a > 0$ is the rate of increase of hosts in the absence of parasitoids, $b > 0$ is the average number of adult female parasitoids emerging from each parasitized host, $\Phi\left(x, y\right)$ represents the probability that a host escapes parasitism, and $1-\Phi\left(x, y\right)$ is the probability of parasitized hosts ($\Phi\left(x, y\right) = \frac{1}{1+f(x)}$ in the model (1.1)). Models of this type were constructed by Thompson in 1922 ^[1], Nicholson and Bailey in 1935 ^[2] and May ^[3]. May's model has the following form:

where the parameter $m > 0$ is the aggregation of parasitoid attacks and $\mu > 0$ describes the efficiency of the parasitoids' search. Note that the model (1.1) is a special case of May's model (1.2) with $m = 1$ and $\mu = f(x_n)$. Here, the function $f$ describes the efficiency of the parasitoids' search. Host-parasitoid interaction represented by the model (1.1) can be considered as a May's host-parasitoid model of biological control where the hosts represent pests.

The essential characteristics of the local and global behavior of May's model (1.2), depending on the values of parameters $m$ and $\mu$, can be found in ^[4]. The case when $m = 1$ and $\alpha > 1$ is fascinating, and when using the Kolmogorov-Arnold-Moser (KAM) theory, the conclusion is reached that the positive equilibrium is stable (but not asymptotically) ^[5]. In ^[6], the author proved that when $\alpha > 1$ and $m > 1$, solutions with the initial conditions in the complement of a bounded subset of the positive quadrant are unbounded. Also, host and parasitoid populations oscillate for these initial conditions with infinitely increasing amplitude. In ^[4], using the KAM theory, the authors investigate the stability of May's host-parasitoid model's solutions with proportional stocking of the parasitoid population. KAM theory applications can be seen in ^[7,8,9,10].

Recently, numerous authors have investigated host-parasitoid models with different characteristics.

Ladas et al., in ^[5], considered the following May's host-parasitoid model:

where $\alpha$ and $\beta$ are positive numbers, and the initial conditions $x_{0}$ and $y_{0}$ are arbitrary positive numbers. This model is a special case of model (1.1). Among other methods, the authors used KAM theory to show that the observed model exhibits very complex behavior.

Jang, in ^[11], presents two general discrete-time host-parasitoid models with Allee effects on the host:

and

where $N_{0}\geq 0$, $P_{0}\geq 0$, $a > 0$, $b > 0$, and the function $g$ satisfies certain conditions.

The author showed that both models exhibit similar asymptotic behavior. The parasitoid population will go extinct if the maximal growth rate of the host population is less than or equal to one, regardless of whether density-dependence parasitism occurs first. If this maximal growth rate exceeds one, the fate of the population is dependent on the initial condition.

Kalabušić and Pilav, and their collaborators considered several host-parasitoid models, especially those with immigration in the population. Among others, they considered the following May's host-parasitoid model with stocking (in ^[4]):

where $a$ and $b$ are positive numbers. Using the KAM theory, the authors investigated the stability of solutions of the May's host-parasitoid model with proportional supply of the parasitoid population. They showed the existence of the extinction point, the limit, and the internal equilibrium point. When the intrinsic growth rate of the host population and the release coefficient are less than one, both populations are extinct. They showed that there is an infinite number of equilibrium limit points, which are non-hyperbolic and stable. They also showed that 1:1 non-isolated resonant fixed points appear under certain conditions, and they described their nature of stability in detail. The stability of the internal equilibrium was demonstrated by using the KAM theory.

By eliminating $y_{n}$ from the first equation of (1.1) and substituting in the second equation, we obtain

This is a crucial feature of this model because in ^[12], the authors noted that proper oscillations in population dynamics can only occur in density-dependent evolution in which delayed negative feedback regulates the evolution. Also, see ^[13].

In Section 2, depending on the parameters $a$ and $b$, we describe the equilibrium points and local and global stability of extinction equilibrium and boundary non-hyperbolic equilibrium points. Also, we show the existence of a 1:1 non-resonant equilibrium point and describe the dynamics of system (1.1) about this equilibrium. In Section 3, we describe the local behavior of the interior equilibrium point. For the case when the interior equilibrium is elliptic, in Section 4, we use the Birkhoff normal form and the twist KAM theorem (^{[12,14,15,16]}) to determine the stability of the interior equilibrium point. Also, we describe a structure that is close to a non-degenerate fixed point $E_{p}$. In Section 5, we apply our result to a special system of difference equations with $f(x) = \frac{1}{1+x}$. Through numerical computation, we confirm our analytic results.

To determine the stability of an elliptic fixed point, we use an appropriate coordinate transformation to simplify the nonlinear terms, that is, to obtain the use of the so-called normal form of the map.

Under certain conditions, system (1.1) has 1:1 non-isolated resonance fixed points for which we thoroughly describe the dynamics. A fixed point of a planar map is said to be 1:1 resonant if the Jacobian matrix of the map at the fixed point is similar to $\left(\begin{array}{cc} 1 & 1 \\ 0 & 1 \end{array} \right)$. A fixed point of a planar map is called isolated if there exists a neighborhood of the fixed point that does not contain any other fixed points. In all other cases, each fixed point is called non-isolated. In ^[17], is given a complete classification of all possible dynamical behavior scenarios that are valid in a neighborhood of non-isolated 1:1 resonant fixed points for planar maps that are real and analytic.

2.   The behavior of the extinction equilibrium point and the boundary equilibrium points

In this paper, we consider only non-negative equilibrium points. The equilibrium points $\left(\overline{x}, \overline{y}\right)$ of system (1.1) satisfies the following system of algebraic equations:

It is easy to see that system (1.1) always has an extinction equilibrium $E_{0} = \left(0, 0\right)$, where both populations become extinct. This equilibrium is unique if $0 < a < 1$ and $b > 0$. For $a > 1$ and $b > 0$, system (1.1) has an interior equilibrium $E_{p} = \left(\frac{a}{bf\left(\overline{x}\right) }, \frac{a-1}{f\left(\overline{x}\right) }\right)$, where $\overline{x}$ is a unique positive solution of the equation $b\overline{x}f\left(\overline{x}\right) -a = 0$ (by assumptions), and where the populations coexist. If $a = 1$ and $b > 0$, then there exist infinitely many boundary equilibriums $E_{\overline{x} } = \left(\overline{x}, 0\right)$, $\overline{x}\geq 0$ of system (1.1), where the host population survives and the parasitoid population becomes extinct.

The map associated with system (1.1) has the following form:

where $T:\left(0, \infty \right) ^{2}\rightarrow \left(0, \infty \right) ^{2}$. It is obvious that $T^{n}\left(\begin{array}{c} 0 \\ y \end{array} \right) = \left(\begin{array}{c} 0 \\ 0 \end{array} \right)$ for $n\geq 1$ and $y > 0$. Also, $T^{n}\left(\begin{array}{c} x \\ 0 \end{array} \right) = \left(\begin{array}{c} a^{n}x \\ 0 \end{array} \right)$ for $n\geq 1$ and $x > 0$, from which we obtain that

(a) $T^{n}\left(\begin{array}{c} x \\ 0 \end{array} \right) \rightarrow \left(\begin{array}{c} 0 \\ 0 \end{array} \right)$, for $0 < a < 1$, $n\rightarrow \infty$,

(b) $T^{n}\left(\begin{array}{c} x \\ 0 \end{array} \right) \rightarrow \left(\begin{array}{c} \infty \\ 0 \end{array} \right)$ for $a > 1$, $n\rightarrow \infty$, and

(c) $T^{n}\left(\begin{array}{c} x \\ 0 \end{array} \right) = \left(\begin{array}{c} x \\ 0 \end{array} \right)$ for $a = 1$.

Based on the Jacobian matrix associated with map (2.2),

we obtain the following results about extinction equilibrium point $E_{0}$.

Lemma 1. The following statements hold for the extinction equilibrium point $E_{0}$:

(i) If $0 < a < 1$, then $E_{0}$ is globally asymptotically stable.

(ii) If $a > 1$, then $E_{0}$ is unstable (a saddle point) with

as the subsets of the stable and unstable manifolds, respectively.

(iii) If $a = 1$, then $E_{0}$ is a non-hyperbolic point, which is stable but not asymptotically stable.

Proof. The Jacobian of the map $T$ at the equilibrium $E_{0} = (0, 0)$ is given by

The eigenvalues of the Jacobian at the equilibrium $E_{0} = (0, 0)$ are $\lambda _{1} = a$ and $\lambda _{2} = 0$, which implies that $E_{0} = (0, 0)$ is locally asymptotically stable for $0 < a < 1$, but is unstable (a saddle point) if $a > 1$ and a non-hyperbolic point for $a = 1$.

(ⅰ) If $0 < a < 1$, then the first equation of system (1.1) implies that $x_{n+1} < ax_{n} < a^{n+1}x_{0}$, which means that $x_{n}\rightarrow 0$ as $n\rightarrow \infty$ (since $x_n\geq 0$ for all $n = 0, 1, ...$). From the second equation of system (1.1), we have that $y_{n+1} < bx_{n}$, which implies that $y_{n}\rightarrow 0$ as $n\rightarrow +\infty$ (since $y_n\geq 0$ for all $n = 0, 1, ...$), that is, $E_{0} = (0, 0)$ is a global attractor. Since $E_{0} = (0, 0)$ is locally asymptotically stable, we conclude that it is globally asymptotically stable.

(ⅱ) The correctness of the statement follows directly from the discussion that precedes this lemma.

(ⅲ) Note that the positive $y$ axis is in the same direction as an eigenspace $E^{s}$. On the other hand, the positive $x$ axis is invariant under the map $T$, and it is in the same direction as an eigenspace $E^{c}$. It means that the positive $x$ axis is a center manifold $\mathcal{W}^{c}$, on which $x_{n+1} = x_{n}$ is valid for all $n = 0, 1, ...$, and where every point is a fixed point of the map $T$. It implies that the equilibrium point $E_{0} = (0, 0)$ is stable, but not asymptotically stable. □

If $a = 1$, then the Jacobian matrix for the equilibrium points denoted by $E_{\overline{x}} = \left(\overline{x}, 0\right)$, $\overline{x} > 0$, have the following form:

whose eigenvalues at equilibrium points are $\lambda _{1} = 1$ and $\lambda _{2} = b\overline{x}f\left(\overline{x}\right)$. It implies that each of the equilibrium points is non-hyperbolic.

Lemma 2. For the non-hyperbolic equilibrium points denoted by $E_{\overline{x}} = \left(\overline{x}, 0\right), \overline{x} > 0$, the following statements are valid:

${\bf{(i)}}$ If $\lambda _{2} = b\overline{x}f\left(\overline{x}\right) > 1$, then $E_{\overline{x}}$ is unstable.

${\bf{(ii)}}$ If $\lambda _{2} = b\overline{x}f\left(\overline{x}\right) < 1$, then $E_{\overline{x}}$ is stable.

${\bf{(iii)}}$ If $\lambda _{2} = b\overline{x}f\left(\overline{x}\right) = 1$, then $E_{\overline{x}}$ is a 1:1 resonant fixed point.

Proof. Namely, it is obvious that the statement under (ⅰ) is valid.

If $\lambda_{2} = b\overline{x}f\left(\overline{x}\right) < 1$, note that an eigenspace $E^{s}$ is in the same direction as the eigenvector $(\frac{\overline{x}f(\overline{x})}{1-b\overline{x}f(\overline{x})}, 1)$. Also, it is easy to see that the positive $x$ axis is invariant under the map $T$ and is in the same direction as an eigenspace $E^{c}$. Thus, the positive $x$ axis is a center manifold $W^{c}$. On this center manifold, it is valid that $x_{n+1} = x_{n}$ for all $n = 0, 1, ..$, and that each point of this map is a stable fixed point. Thus, the each boundary equilibrium $E_{\overline{x}} = \left(\overline{x}, 0\right)$ of the map $T$ is stable, but not asymptotically stable.

If $\lambda _{2} = b\overline{x}f\left(\overline{x}\right) = 1$, then $\lambda _{1} = \lambda _{2} = 1$ and the equilibrium points denoted by $E_{\overline{x}} = \left(\overline{x}, 0\right)$, $\overline{x} > 0$, become $E_{\overline{x} } = \left(\frac{1}{bf\left(\overline{x}\right) }, 0\right)$. The Jacobian matrix at the equilibrium points is of the form

which is similar to the $\left(\begin{array}{cc} 1 & 1 \\ 0 & 1 \end{array} \right)$ matrix because

Thus, $E_{\overline{x}} = \left(\frac{1}{bf\left(\overline{x}\right) }, 0\right)$ is a 1:1 resonant fixed point of $T$ for all $b > 0$. To study the dynamical behavior in a neighborhood of the 1:1 resonant fixed point, we will use a result from ^[17]. By performing the following change of variables: $x\longmapsto x+\overline{x}$, $y\longmapsto y$, the equilibrium point $E_{\overline{x}} = \left(\frac{1}{ bf\left(\overline{x}\right) }, 0\right) = \left(\overline{x}, 0\right)$ shifts to $\left(0, 0\right)$. Now, we have

and

Conjugating by $P$ yields

Denote

Calculating the partial derivatives of $\phi \left(x, y\right)$ and $\psi \left(x, y\right)$, we get

Thus, in order to apply Theorem 2 from ^[17], it must be that $\overline{x} f\left(\overline{x}\right) -b = 0$. Notice that $b =$ $\overline{x}f\left(\overline{x}\right)$ and $b\overline{x}f\left(\overline{x}\right) = 1$ implies that $b = \overline{x}f\left(\overline{x}\right) = 1$. Assuming that $b = \overline{x}f\left(\overline{x}\right) = 1$, we have $\psi \left(x, y\right) = -y\frac{\left(x+\overline{x}+y\right) f\left(x+\overline{x}\right) -1}{ yf\left(x+\overline{x}\right) -1}$ and $\varphi \left(x\right) = 0$. Also,

To apply Theorem 2 ^[17], the last expression should be developed into a power series by taking some specific function $f$ that satisfies the above conditions. For example, consider the following map:

which we get from (2.2) for $f\left(x\right) = \frac{1}{ \alpha x+1}$ $\left(0 < \alpha < 1\right)$.

Then, we obtain

Therefore, $l = 1$ and $Q\left(0\right) = \left(1-\alpha \right) ^{2} > 0$, and, by Theorem 3 ^[17], the dynamical behavior of (2.3) near $\left(\frac{1}{1-\alpha }, 0\right)$ corresponds to (ⅰ) of Figure 4 from ^[17]. That is, there are four sectors, in either clockwise or counterclockwise orientation, which are of elliptic, attracting parabolic, hyperbolic, and repelling parabolic type. Also, the set $\Phi \backslash \{(0, 0)\}$ has two connected components $\Phi _{u}$ and $\Phi _{s}$ such that $S_{-n}(u, v)\rightarrow (0, 0)$ for every $(u, v)\in \Phi _{u}$ and $S_{n}(u, v)\rightarrow (0, 0)$ for every $(u, v)\in \Phi _{s}$, where $\Phi$ is a real analytic curve that represents the set of fixed points of $S$. See Figure 1. □

3.   Local behavior of the interior equilibrium point

In this section, we consider the local behavior of the interior equilibrium $E_{p} = \left(\frac{a}{ bf\left(\overline{x}\right) }, \frac{a-1}{f\left(\overline{x}\right) } \right)$ that exists for $a > 1$. The Jacobian matrix evaluated at point $E_{p}$ is given by

Since $\overline{x} = \frac{a}{bf\left(\overline{x}\right) }$, then

Note that $\det J_{T}\left(E_{p}\right) = 1$, and the complex conjugate eigenvalues of $J_{T}\left(E_{p}\right)$ are given by

Lemma 3. If $a > 1$, then the following statements are valid for the equilibrium point $E_{p} = \left(\overline{x}, \frac{a-1}{f\left(\overline{x} \right) }\right)$:

(i) If $f\left(\overline{x}\right) +f^{\prime }(\overline{x}) \overline{x} < 0$, then $E_{p}$ is a saddle point.

(ii) If $f\left(\overline{x}\right) +f^{\prime }(\overline{x}) \overline{x} = 0$, then $E_{p}$ is a non-hyperbolic point of parabolic type (1:1 resonant fixed point).

(iii) If $f\left(\overline{x}\right) +f^{\prime }(\overline{x}) \overline{x} > 0$, then $E_{p}$ is a non-hyperbolic point of elliptic type.

Proof. (ⅰ) Let us keep in mind the condition $a > 1$ and the fact that $f\left(x\right)$ is a positive decreasing function for all $x$, that is, $f(x) > 0$ and $f^{\prime }(x) < 0$ for all $x$. Now, under the assumption that $f\left(\overline{x }\right) +f^{\prime }(\overline{x})\overline{x} < 0$, we will prove that $\lambda _{+} > 1$ and $-1 < \lambda _{-} < 1$. Namely,

which is satisfied since $\left(a-1\right) \left(f\left(\overline{x} \right) +f^{\prime }(\overline{x})\overline{x}\right) < 0$.

Also, we obtain that

which is satisfied.

On the other hand, the condition $\lambda _{-} > -1$ is equivalent to

that is,

which is satisfied. Thus, $E_{p}$ is a saddle point.

(ⅱ) If $f\left(\overline{x}\right) +f^{\prime }(\overline{x})\overline{x} = 0$, then $\lambda _{\pm } = \frac{f\left(\overline{x}\right) -f^{\prime }(\overline{x})\overline{x}}{2f\left(\overline{x}\right) } = \frac{1}{2}\left(1-\frac{f^{\prime }(\overline{x})\overline{x}}{f\left(\overline{x}\right) } \right) = 1$. The Jacobian matrix at the equilibrium point $E_p$ has the following form:

which is similar to the $\left(\begin{array}{cc} 1 & 1 \\ 0 & 1 \end{array} \right)$ matrix because

Thus, $E_{p} = \left(\frac{a}{bf\left(\overline{x}\right) }, \frac{a-1}{ f\left(\overline{x}\right) }\right)$ is a 1:1 resonant fixed point of $T$ for all $b > a > 1$. Unfortunately, here, we cannot successfully carry out the procedure based on Theorem 2 in ^[17] for the case of the equilibrium points $E_{\overline{x}}$ since the condition $\phi \left(0, 0\right) = 0$ is not satisfied (here, $\phi \left(0, 0\right) = -\frac{a-1}{f\left(\overline{x}\right) ab^{2}} < 0$).

The proof of claim (ⅲ) is obvious and will be omitted. □

4.   Stability of the elliptic interior equilibrium point via KAM theory

The following considerations will be based on the assumption that

We can use the logarithmic change of variables to show that the map $T$ transforms into an area-preserving map with a non-degenerate elliptic fixed point $(0, 0)$. It is well known that the map $T$ is area-preserving if and only if the determinant of the Jacobian matrix of the map $T$ is equal to $1$ at every point in $(0, \infty)^{2}$. Under the logarithmic coordinate change $(x, y)\rightarrow (u, v)$, the fixed point $E_{p}$ becomes $(0, 0)$. By using the substitutions given by

system (1.1) transforms into the following system:

i.e.,

or, equivalently,

Therefore, denote

Lemma 4. The map $K$ has the following properties:

a) $K$ is globally area-preserving;

b) the map $K$ in the $\left(x, y\right)$ coordinates has an elliptic fixed point $E_{p}$ if $a > 1$ and $b > 0$. The corresponding fixed point $(x, y)$ in the $\left(u, v\right)$ coordinates is $\left(0, 0\right).$

Proof. The Jacobian matrix for the map $K$ is given by

Now, it is easy to see that $\det J_{K}\left(u, v\right) = 1$ for all $(u, v)\in (0, \infty)^2$, which proves statement a above. Also,

and, by (4.1),

The equation

is the characteristic equation of the matrix $J_{K}\left(0, 0\right)$ with the corresponding characteristic roots $\lambda$ and $\overline{\lambda }$, where

From (4.1) and given that $f^{\prime }\left(\overline{x}\right) < 0$, it follows that the expression under the square root is positive. Using the equality $\overline{y} = \frac{a-1}{f\left(\overline{x}\right) }$, the Eq (4.4) can also be written as

where

Since $\left\vert \lambda \right\vert = 1$, the point $E_{p}$ is an elliptic fixed point. Using substitutions $u_{n} = \ln \frac{x_{n}}{\overline{x}}$ and $v_{n} = \ln \frac{y_{n}}{\overline{y}}$, it is obvious that logarithmic coordinate change transforms $E_{p}(\overline{x}, \overline{y})$ into $\left(u, v\right) = (0, 0)$. Hence, the statement b above is also valid. □

Now, we apply the KAM theory in a small neighborhood of an elliptic fixed point to determine its stability. For this purpose, we derive the Birkhoff normal form near the elliptic fixed point, and then we verify the non-resonance and twist conditions.

where

Since $a > 1$ and $f^{\prime }\left(\overline{x}\right) < 0$, it is obvious that $\operatorname{Im}\left(\lambda ^{2}\right) > 0$ and $\operatorname{Im}\left(\lambda ^{3}\right) > 0$, so $\lambda ^{2, 3}\neq 1$. Also, we have that $\lambda ^{4}\neq 1$. Namely, if we assume the opposite, i.e., $\lambda ^{4} = 1$, then it should follow that $I_{2} = 0$ and $\frac{R_{2}}{2a^{4}} = 1$. Because $a > 1$ and $f^{\prime }\left(\overline{x}\right) < 0$, $I_{2} = 0$ is only possible if

Equation (4.6) has only one positive solution:

for $a > \sqrt{2}+1$. Above, we used the fact that $\overline{y} = \frac{a-1}{f\left(\overline{x}\right) }$. If (4.7) is satisfied, then $\lambda ^{4} = 1$ implies that

which is equivalent to

Since $\left(a-\overline{x}\overline{y}f^{\prime }\left(\overline{x} \right) +1\right) ^{2} > 0$ and $\overline{x}\overline{y}f^{\prime }\left(\overline{x}\right) -3a-1 < 0$, then (4.8), using (4.7), is equivalent to

which is impossible. Therefore, $\lambda ^{4}\neq 1$.

By using $\overline{y} = \frac{a-1}{f\left(\overline{x}\right) }$, the matrix of the linearized system at the origin is given by

The eigenvalue of (4.9) is of the form

i.e.,

To obtain the Birkhoff normal form of system (4.2), we will expand the right-hand sides of the equations of the system (4.2) at the equilibrium point $\left(0, 0\right)$, as follows:

from which we obtain

By using the eigenvector

the associated matrix can be obtained as follows:

where

and $\det P = 1.$

Now, we change the coordinates as follows:

i.e.,

and bring the linear part into the Jordan normal form. The system, given the new coordinates, becomes

with

where the coefficients $\alpha _{20}$, $\alpha _{11}$, $\alpha _{02}$, $\alpha _{30}$, $\alpha _{21}$, $\alpha _{12}$, $\alpha _{03}$, $\beta _{20}$, $\beta _{11}$, $\beta _{02}$, $\beta _{30}$, $\beta _{21}$, $\beta _{12}$, $\beta _{03}$ are given in Supplementary A.

Now, the complex coordinates $z, \overline{z} = \widetilde{u}\pm i\widetilde{v}$ yield the complex form of the system:

Using the Mathematica package and

we obtain the coefficients as in the forms shown in Supplementary B.

The above normal form yields the approximation

with $c_{1} = i\lambda \alpha _{1}$, where $\alpha _{1}$ is the first coefficient. The coefficient $c_{1}$ can be evaluated by using the following formula:

derived by Wan in the context of Hopf bifurcation theory.

We apply

and

A tedious symbolic computation done with Mathematica yields

where

Finally, after painstaking calculations in Mathematica, we get that

i.e.,

where

which implies that $\tau _{1}\neq 0$ if (4.1), $a > 1$, $f\left(\overline{x}\right) > 0$, and $f^{\prime }\left(\overline{x}\right) < 0$ hold. It means that $E_{p}$ is a non-degenerative fixed point. By using Lemma 4 and the previous conclusion, we apply Theorems 2.26 and 2.27 from ^[18] for $q = 4$, $s = 1$, and $\tau _{1}\neq 0$ to get the following result.

Theorem 5. Assume that $f\in C^{1}\left(\left[0, +\infty \right) \right)$ is a strictly decreasing function such that $f:\left[0, +\infty \right) \rightarrow \left(0, +\infty \right)$ and $f\left(0\right) > 0$. If (4.1) holds, where $\overline{x}$ is an equilibrium point of system (1.1), then $E_{p}$ is a stable equilibrium point of the map corresponding to system (1.1).

The importance of Theorem 5, from a biological point of view, is that Theorem 5 guarantees that all orbits starting near equilibrium $E_{p}$ of the system (1.1) are on the invariant curves surrounding the interior equilibrium point $E_{p}$.

By using Theorem 4.3 from ^[13], we get the following theorem that describes the structure close to a non-degenerate fixed point $E_{p}$ in more detail.

Theorem 6. Assume that $f\in C^{1}\left(\left[0, +\infty \right) \right)$ is a strictly decreasing function such that $f:\left[0, +\infty \right) \rightarrow \left(0, +\infty \right)$ and $f\left(0\right) > 0$. If (4.1) holds, where $\overline{x}$ is an equilibrium point of system (1.1), then, in every neighborhood of $E_{p}$, periodic points of $T$ with arbitrarily large periods exist.

5.   Special case with numerical simulations

In this section, we apply Theorem 5 to system difference equations of the form (1.1). We consider the following system:

where $a, b$ are positive numbers. System (1.1) for $f\left(x\right) = \frac{1}{1+x}$ becomes system (5.1).

The equilibrium points $\left(\overline{x}, \overline{y}\right)$ of the system (1.1) satisfy the following system of algebraic equations:

It is easy to see that the system (5.1) always has an extinction equilibrium $E_{0} = \left(0, 0\right)$, where both populations become extinct. For $b > a > 1$, the system (5.1) has an interior equilibrium $E_{p} = \left(\frac{a}{b-a}, \frac{\left(a-1\right) b}{b-a}\right)$, where the populations coexist. If $a = 1$, the system (5.1) has another boundary equilibrium $E_{\overline{x}} = \left(\overline{x}, 0\right)$, $\overline{x}\in\mathbb{R}^{+}$, where the host population survives, and the parasitoid population becomes extinct.

Lemma 7. (i) If $0 < a < 1$ or $0 < b\leq a\neq 1$, then system (5.1) has a unique equilibrium point, extinction equilibrium $E_{0}$.

(ii) If $1 < a < b$, then system (5.1) has two equilibrium points: the extinction equilibrium point $E_{0}$ and the interior equilibrium point $E_{p}$.

(iii) If $a = 1$ and $b > 0$, then system (5.1) has infinitely many boundary equilibrium points denoted by $E_{\overline{x}} = \left(\overline{x}, 0\right)$, $\overline{x}\geq 0$.

For the stability of the extinction equilibrium point $E_0$, see Lemma 1, but, for non-hyperbolic equilibrium points denoted by $E_{\overline{x}} = \left(\overline{x}, 0\right)$, $\overline{x} > 0$, the following result is valid (also, see Lemma 2).

Lemma 8. (i) If $0 < b\leq1$, then $E_{\overline{x}}$ is stable.

(ii) If $b > 1$, then

$a.$ $E_{\overline{x}}$ is stable for $\overline{x} < \frac{1}{b-1}$,

$b.$ $E_{\overline{x}}$ is unstable for $\overline{x} > \frac{1}{b-1}$,

$c.$ $E_{\overline{x}} = (\frac{1}{b-1}, 0)$ is a 1:1 resonant fixed point.

Since $f\left(\overline{x}\right) +f^{\prime }(\overline{x})\overline{x} = \frac{\left(b-a\right) ^{2}}{b^{2}} > 0$, for $\overline{x} = \frac{a}{b-a}$, the condition (4.1) is satisfied; thus, $E_{p}$ is an elliptic fixed point and (4.10) has the following form:

This implies that $\tau _{1} < 0$ if $b > a > 1$. By Theorems 5 and 6, we have the next result.

Theorem 9. Assume that $1 < a < b$. Interior equilibrium point $E_{p}$ of (5.1) is an elliptic fixed point. Let $T$ be the map associated with system (5.1). Then, periodic points of the map $T$ with arbitrarily large periods exist in every neighborhood of $E_{p}$. In addition, $E_{p}$ is a stable equilibrium point of (5.1).

Figure 2 shows the phase portraits of the orbits of the map $K$ with the non-degenerate elliptic fixed point $(0, 0)$ created by transforming the map $T$ associated with system (5.1) for $a = 2$ and $b = 4$. Neither of these two plots shows any self-similar characteristic. Figure 3 shows the bifurcation diagram and corresponding Lyapunov coefficients for $a\in(0.9, 3.0)$ and $b = 3$. In both figures, it is possible to see the complexity of the behavior of the orbits for a neighbor of the positive equilibrium point.

The eigenvalues denoted as $\lambda _{\pm }$ at the elliptic fixed point are of the form $\lambda = e^{i\theta }$ with $\theta = \arccos \frac{b-a+ab+a^{2}}{2ab}$ and $0 < \theta < \frac{\pi }{2}$. Thus, in the case that $b = 4$, the period of the motion around the fixed point must be greater than $\frac{2\pi }{ \theta } = 12,433$, so the minimal possible period for a periodic orbit in a neighborhood of the elliptic fixed point is $13$; similarly, in the case that $b = 3$, the period of the motion around the fixed point must be greater than $\frac{2\pi }{\theta } = 14,762$; see Figures 4 and 5.

Figures 6 and 7 show the times-series plots for the components $x_n$ and $y_n$ for the map $T$.

6.   Conclusions

In this paper, we investigated the stability of a general host-parasitoid model of the form (1.1) with the function $f$, the properties of which we assumed as detailed in Section 1. We confirmed the existence of extinction, interior, and boundary equilibrium points. When the rate of increase of the hosts is less than 1 ($0 < a < 1$), the extinction point of the equilibrium is globally asymptotically stable, which means that the extinction of both populations occurs. When this rate of increase is $a = 1$ and one of the boundary equilibrium points is 1:1 resonant, the other equilibrium points to the right of it are unstable, while those to the left are stable. Such behavior is illustrated in Figure 1 in the special case when $f(x) = \frac{1}{\alpha +1}$, $0 < a < 1$. For $a > 1$, we showed that the interior point of the equilibrium is elliptic, and that the corresponding map $T$ associated with the model (1.1) has the property of being area-preserving. After calculating the Birkhoff normal form, by using KAM theory, we concluded that the internal equilibrium point is stable. As an immediate consequence, we obtained a conclusion about the existence of periodic points with an arbitrary period in the vicinity of this elliptic equilibrium point. Finally, taking the special case of $f(x) = \frac{1}{1+x}$, that is, when the model (1.1) takes the form (5.1), we confirmed the previously obtained general results. In this case, our numerical simulations visually show the answer to the central question of biological significance for the observed model, which is demonstrated by the qualitative behavior of populations (hosts and parasitoids) over time, especially the stability/instability of trajectories. If the initial state of the population represents a point on a periodic orbit, on an invariant curve, or on some other invariant set, then the future evolution of the population will remain confined to that invariant set for all time. If the initial conditions correspond to a point between two invariant curves, the future evolution (the corresponding orbit) will forever remain bounded between these invariant curves. In a rough sense, the behavior of this population is stable but not asymptotically stable. On the other hand, if the initial condition lies on some invariant curve, the evolution of populations can be regular. However, it can also be chaotic if the initial condition lies in the stochastic region.

Our model is general since we also consider an arbitrary function $f$ as an integral part of the probability function that is associated with the host avoiding parasitism. The function satisfies the conditions of the natural properties that arise from their biological meaning. In this way, this model encompasses all similar models that use such specific probability functions for parasitoid avoidance and release. Therefore, the results obtained for the concrete form of the function $f$ are special cases of the results obtained in this study. This means that, if the population of parasitoids released into the existing population decreases or increases with other system parameters, it significantly determines the model's local and global dynamics.

The obtained theoretical results can be used for specific situations in biological control because they can help managers to find pest control strategies, etc. Let us emphasize that, instead of the function $f$ used in this work, some known host escape probabilities can be observed. The dynamics of the system, as shown by the results of numerical simulations, largely depend on the forms of these functions and system parameters. This gives us ideas about the possibilities of further research regarding these models (e.g., using the general Beverton-Holt function).

Also, we performed several significant visual simulations by using the Mathematica software package.

Use of AI tools declaration

We did not use artificial intelligence tools in the creation of this article.

Acknowledgments

The authors are grateful to the anonymous referee for number of helpful and constructive suggestions which improve the presentation of results.

This research was funded by FMON of Bosnia and Herzegovina, project number 01-3840-Ⅵ-1/23.

Conflict of interest

The authors declare no conflicts of interest.

Supplementary A

Supplementary B