Global stability of the interior equilibrium and the stability of Hopf bifurcating limit cycle in a model for crop pest control

1.   Introduction

The ruin of crops by pest invasions is a serious worldwide crisis not only in farming areas but also in woodland ecosystems. This issue connected with pests has been recognized since the cultivation of crops started. Approximately 42% of the world's food supply is destroyed because of pests ^[1]. Recently, the biological control of pests has been earning more attention among experimenters, and its practical application is also rising in the crop field. This method seeks to reduce the reliance on pesticides by emphasizing the contribution of biological control agents where living organisms are only used to control pests. Managing pests via chemical pesticides is less expensive and destroys pests rapidly but generates high environmental loss. On the other hand, natural control is a long and expensive method to use, but with very little ecological loss ^[2,3].

In this study, we focus on the application of biopesticides for crop pest control. A mathematical model on this system can help us to identify the main parameters. The stability of model system shows complex dynamics such as bifurcation and chaos. The global stability of equilibrium point and the stability of a limit cycle around it are important factors in understanding the behaviors that can be seen in a dynamical system. For example, in a biological system, the stability of a limit cycle can determine whether a population of organisms will grow or decline over time. In a biological system, the stability of a limit cycle can determine whether the system will oscillate periodically or exhibit unstable behavior ^[4,5].

In a dynamical system, a limit cycle refers to a periodic orbit that a system can exhibit ^[6]. The stability of a limit cycle determines how the system behaves near the limit cycle ^[7]. There are two types of stability for limit cycles: asymptotic stability and structural stability. Asymptotic stability means that the system tends toward the limit cycle as time goes to infinity ^[8]. This means that any initial conditions near the limit cycle will approach it over time. Asymptotic stability can be further classified into stable limit cycles and unstable limit cycles. A stable limit cycle is one where the system is attracted toward the limit cycle, while an unstable limit cycle is one where the system is repelled away from the limit cycle ^[9]. Structural stability, on the other hand, means that the qualitative behavior of the system is preserved under small perturbations to the system's parameters. In other words, if a system exhibits a limit cycle under certain conditions, then it will continue to do so under small changes to those conditions ^[10,11].

Mathematical modeling and analysis helps in understanding the dynamics of the system under consideration ^[12]. Results obtained from the modeling approach can help in proposing proper pest control strategies ^[13]. Dynamical properties, such as Hopf bifurcation, chaos, limit cycle etc., have seen studied by researchers ^[6]. Mathematical models for pest control using biopesticides have been proposed and studied by researchers ^{[14,15,16,17]}. Pest control models are developed to analyze specific dynamics, such as Hopf bifurcation via the occurrence of limit cycles ^{[18,19,20,21]}. The authors of ^[22] and ^[23] have studied the occurrences of Hopf bifurcating periodic solutions. Little research are focused on the effect of environmental fluctuations on a sterile insect release method ^[24]. Mathematical models for pest management using biological control strategies have been designed and analyzed by many researchers ^[16,25,26]. The local stability of equilibrium points and the occurrence of Hopf bifurcation has been analyzed in the articles. But there are a few articles that focus on the global stability of equilibrium points and the stability of bifurcating periodic solutions which are important properties of a dynamical system ^[27].

In ^[28], Chowdhury et al. have proposed a model for crop pest control using an integrated approach (i.e., both biopesticides and chemical pesticides are used) in an optimal manner. They have analyzed the local stability of different equilibria and provided the existence of Hopf bifurcation. But the global stability and stability of a limit cycle were not addressed. Thus, in this research, we have derived a model for pest control from the mathematical model proposed in ^[28] to study the dynamics of the crop-pest system in the presence of biopesticides. The model system is analyzed both analytically as well as numerically for the occurrence and stability of a Hopf bifurcating limit cycle using Poore's condition. Also, the global stability of the interior equilibrium point is analyzed using a suitable Lyapunov function.

The paper is organized as follows. In section 2, the formulation of the model is given. In section 4, we will analyze the global stability of the endemic equilibrium point and the occurence of Hopf bifurcation. The stability of the bifurcating limit cycle around coexisting equilibrium is analysed in section 5. Simulations are carried out in section 6 to substantiate the analytical findings. The final section contains a discussion of the main results to conclude the the paper.

2.   Formulation of the mathematical model

In this research, we have considered the mathematical model established in ^[28] and derive the mathematical model to study the global stability of the interior equilibrium and stability of the limit cycle. We provide the derivation of the model as follows.

The following hypotheses are made to derive the mathematical model:

- As system variables, the biomass of crop is denoted as $X(t)$, the susceptible pest population as $Y(t)$, the infected pest as $I(t)$ and the biopesticides (viruses) as $V(t)$.

- Logistic growth for the biomass of the crop is considered with net growth rate $r_1$ and carrying capacity $K$. $\alpha$ is the consumption rate of the crop biomass by susceptible pests.

- In biological control of pest, biopests (generally viruses) are sprayed on the plantation. In this process, susceptible pests are affected and become the infected ones. $\kappa$ is the rate of replication of the virus by lysis and $\mu_V$ is the mortality rate of the biopesticide (virus). $\pi_v$ is the constant rate of spray of virus in the environment. $\lambda_1$ is the reduction rate constant of the free virus.

- Pest consumes the crop resource at a rate $\alpha$ which is converted to the susceptible pest population with a maximum growth rate $r(X)$. Here, $r(X)$ is dependent on the density of the crop biomass and is also assumed to be a scalar multiple of the biomass of crop $X$ due to the consumption of the crop resources limited at the maximum crop cultivation capacity [$r(X)$ = $\alpha r_2 X$]. $\alpha r_2$ is the maximum growth rate of the susceptible pest. Hence the susceptible are growing at the rate governed by the consumption of crop resources, which is growing maximally at a rate $r_1$.

- The carrying capacity of the susceptible pest is assumed as $k_s$. Here, $k_s$ is dependent on the carrying capacity of the crop biomass $K$ factored with a constant term $c > 1$, i.e., $k_s$ = $cK$.

Based on the above assumptions, the following model is obtained:

with initial conditions as

The region given below is positively invariant region for the system (2.1):

where $M_1$ = max$\{K, J(0)\}$, $M_2 = \frac{\alpha r_2 M_1 c K}{4 \xi}$, and $M_3 = \frac{\pi_v+\kappa \xi M_2}{\mu_v}.$

3.   Existence of equilibria and their local stability

In this section, we analyse the existence of possible equilibrium points of the model system (2.1). Then we study the local stability properties.

3.1.   Existence of equilibria

Analyzing the system (2.1), we get three feasible steady states, namely,

(ⅰ) the plant-pest free equilibrium point $E_1(0, 0, 0, \frac{\pi_v}{\mu_v})$,

(ⅱ) the pest-free equilibrium point $E_2(K, 0, 0, \frac{\pi_v}{\mu_v})$, and

(ⅲ) the interior equilibrium point $E^*(X^*, S^*, I^*, V^*)$, where,

and $S^*$ is the positive root of the following cubic equation:

with,

Equation (3.1) is a cubic polynomial equation, thus a real root always exists. Also, if $\lambda_1 > \kappa\lambda$ holds then $m_1 < 0, \; m_2 > 0$, and $m_3 > 0$. Now $m_4 > 0$ holds when $\lambda \pi_v > \alpha r_2K\mu_v$. Hence, using Descartes' rule of signs, we have the following proposition.

Proposition 1. For $\lambda_1 > \kappa\lambda$ and $\lambda \pi_v > \alpha r_2K\mu_v$, there always exists a unique interior equilibrium point $E^*$.

3.2.   Local stability and Hopf bifurcation

Here, we check the local stability of the equilibria of the system (2.1). For this analysis, we need the Jacobian matrix of the system at any equilibrium point $E(X, S, I, V)$ is given by

where $J_{22} = \alpha r_2 X (1-\frac{2S+I}{cK})-\lambda V$.

At $E_1$, the Jacobian matrix gives the following characteristic equation in $\rho$:

whose eigenvalues are $r > 0$, $-\frac{\pi_v \lambda}{\mu_v} < 0, -\xi < 0$, and $-\mu_v < 0$. Thus, one eigenvalue is always positive and, consequently, the axial equilibrium, $E_0$, is unstable.

The Jacobian matrix of the system at pest-free equilibrium point $E_2(K, 0, 0, \frac{\pi_v}{\mu_v})$ which satisfies the following equation:

Eigenvalues of the above matrix are given as $-r$, $\alpha r_2 K-\frac{\pi_v \lambda}{\mu_v}$, $-\xi$, and $-\mu_v$. Clearly, three eigenvalues are negative and remaining eigenvalue will be negative if

The characteristic equation of $J_{E^*} = [J_{ij}]_{4\times4}$ is

where

Here,

According to the Routh-Hurwitz criterion, characteristic equation have roots with negative parts if

From Proposition 1 and from the above analysis, the following theorem is obtained.

Theorem 1. In the system (2.1),

(i) plant-pest free equilibrium $E_1$ is always unstable,

(ii) pest-free equilibrium $E_2$ is stable if $\alpha r_2\mu_v K < {\pi_v\lambda}$ and unstable otherwise,

(iii) interior equilibrium $E^*$ exists if $\alpha r_2\mu_v K < {\pi_v\lambda}$, i.e., when $E_2$ becomes unstable. $E^*$ is stable if the conditions in (3.10) are satisfied.

Now, we shall analyse the conditions for which $E^*$ enters into Hopf bifurcation as a model parameter varies over $\mathbb{R}$. We consider the Hopf bifurcation as a function of the generic bifurcation parameter $\eta\in \mathbb R.$

Let $\Psi$ : $(0, \infty)\rightarrow \mathbb R$ be the following continuously differentiable function of $\eta$:

The assumptions for Hopf bifurcation to occur are the usual ones and these require that the spectrum $\sigma(\eta) = \{\rho :D(\rho) = 0\}$ of the characteristic equation is such that

(A) There exists $\eta^*\in(0, \infty)$, at which a pair of complex eigenvalues $\rho(\eta^*), \bar{\rho}(\eta^*)\in\sigma(\eta)$ are such that

and the transversality condition

(B) All other elements of $\sigma(\eta)$ have negative real parts.

Thus we have the following theorem ^[28].

Theorem 2. The system (2.1) around the interior equilibrium $E^*$ enters into Hopf bifurcation at $\eta = \eta^*\in(0, \infty)$ if and only if

(i) $\Psi(\eta^*) = 0,$ and

(ii) $\sigma_1^3\sigma_2'\sigma_3(\sigma_1-3\sigma_3) > 2(\sigma_2\sigma_1^2-2\sigma_3^2)(\sigma_3'\sigma_1^2-\sigma_1'\sigma_3^2)$,

and all other eigenvalues have negative real parts, where $\rho(\eta)$ is purely imaginary at $\eta = \eta^*$.

4.   Global stability of interior equilibrium $E^*$

In the this section, we analyse the global stability of the coexisting equilibrium of system (2.1).

For the global stability, we choose the Lyapunov function as follows:

Here, $c_i > 0, \; i = 1, 2, 3, 4,$ are so chosen that $\dot\psi$ is negative definite. Now, derivative of $\psi$ along the solution of the equation $\dot{X}(t) = J_{E^*} X(t)$, where $X(t) = (X(t), S(t), I(t), V(t))^T$, is as follows

Thus symmetric matrix corresponding to $\dot{\psi}$ is given by

Here,

The positive equilibrium $E^*$ is locally asymptotically stable if $\dot{\psi}$ is negative definite which implies the matrix $M$ must be negative definite. But the matrix $M$ will be negative definite if all of the principal minors of odd rank are negative and all of the principal minors of even rank are positive. This gives rise to the following four conditions:

Clearly, first condition is automatically satisfied. Now, if we choose $c_2, c_3$ in such a way that

i.e.,

then the condition (ⅱ) and (ⅲ) are satisfied, and thus the condition (ⅳ) reduces to the following form:

The above condition is satisfied for any suitable large value of $c_1$ as first term in the above inequality is positive and other two terms are negative. The negative terms do not contain $c_1$.

Thus, we have the following theorem for the global stability of $E^*$.

Theorem 3. The equilibrium point $E^*(X^*, S^*, I^*, V^*)$ of the system (2.1) is locally asymptotically stable if the following condition holds:

where $m_{11}$ and $m_{22}$ are given in (4.2) and $c_1, c_2, c_3$, and $c_4$ are given in (4.5).

5.   Stability of bifurcating limit cycle

To investigate the orbital stability of the Hopf-bifurcating periodic solution, ref30's condition has been followed ^[29]. According to ref30's sufficient condition, the supercritical and subcritical nature of the Hopf-bifurcating periodic solution is determined respectively by the positive and negative sign of real part of the number $\digamma$, where $\digamma$ is defined by:

where the repeated indices within each term imply a sum notation and all the derivatives of $f^l$ are evaluated at the equilibrium $E^*$. $J_{E^*}$ is the variational matrix of the system (2.1) calculated at $E^*$. $u = (u_1, u_2, u_3, u_4)$ and $v = (v_1, v_2, v_3, v_4)^T$ are the left and right eigenvectors respectively of $E^*$ with respect to eigenvalues $i\omega_0$. So positivity of the real part of the above expression in parenthesis really indicates the orbital stability of the periodic solution arising out of Hopf bifurcation.

We rewrite our system (2.1) in the following form:

where $x = (X, S, I, V)$, $f = (f^1, f^2, f^3, f^4)^T$, and $f^l$, $l = 1, 2, 3, 4$ are right sides of system (2.1) i.e. $f^1 = r_1 X(1-\frac{X}{K})-\alpha XS$ etc. Now, all the second and third-order derivatives of $f^l(l = 1, 2, 3, 4)$ are as follows:

Now we calculate $M_{\omega_0} = (J_{E^*}-2i\omega_0 I)^{-1}$

Here,

and

Now, if we put $\omega_0 = 0$ in the above expressions, we get the component of $M = (J_{E^*})^{-1}.$

In the next section we find out the left eigenvector and right eigenvector of the variational matrix $J_{E^*}$ with respect to the eigenvalue $i\omega_0$, i.e., we calculate row vector $u = (u_1, u_2, u_3, u_4)$ and column vector $v = (v_1, v_2, v_3, v_4)^T$ such that

Solving the first equation of (5.5), we find the left eigenvector $u = (u_1, u_2, u_3, u_4)$ where

Solving the second equations of (5.5), we find the right eigenvector $v = \eta(v_1, v_2, v_3, v_4)^T$ where

Now for $uv = 1$ we get,

Here,

where,

Writing the expression of (5.1) in detail, we have the following:

The first term:

The second term:

The third term:

where

Putting the value of all second- and third-order derivatives of $f^{l}(l = 1, 2, 3, 4)$, $u, v$, and the components of matrix $M$ and $M_{\omega_0}$ in the first term, and placing the second term and third term in terms of the parameters of the model, we calculate the expression (5.1) and the sign of the real part of this expression. This in turn indicates the orbital stability of the limit cycle arising from Hopf bifurcation.

6.   Numerical simulations

In this section, we have analyzed the dynamics of the model system using numerical simulations in MATLAB. The analytical results obtained in the previous sections are verified using numerical calculations. Results are plotted in the figures and discussed below.

In Figure 1, the numerical solution of the model is presented. As the value of the parameter is enhanced, oscillation in the solutions is increased. Solutions become periodic when the rate is higher ($\lambda = 0.00012$).

Figure 2 shows the global stability of interior equilibrium point $E^*$ for $\lambda = 0.00008$. All phase portraits are converging to the same interior steady point $(E^*(35.98, 56.07, 28.92, 1031)$.

Figure 3 depicts the bifurcation diagram of the system. We have plotted the maximum and minimum values of the periodic solutions. When the bifurcation parameter $\lambda$ crosses the critical value $\lambda^* = 0.0000965\;$(approx.), the system bifurcates into periodic oscillations.

In Figure 4, we see the orbital stability of the bifurcating periodic solution for a fixed value of the parameter $\lambda$. We observe that when $\lambda$ passes through the value 0.0000965 (approx.), the interior equilibrium $E^*$ bifurcates toward a periodic solution (see Figure 3). From this figure, we conclude that the bifurcating limit cycle is stable (supercritical).

7.   Discussion and conclusions

Control of pest attacks is an important aspect in agriculture to obtain healthy crops as well as high yield. Mathematical modeling helps in identifying the parameters important for crop pest management. In this paper, we have considered a model for pest control using biological agents and observed the effects of biopesticide in controlling the pest. We have explored the global stability of the interior equilibrium point $E^*$. Applying Poore's criteria, we studied the stability of the limit cycle around the interior equilibrium point.

We have shown how the dynamics changes with the increase in the value of the parameter $\lambda$ (the infection rate of the pest by the virus) of the system. The model reveals that infection can be sustained only above a threshold value of the consumption rate $\lambda$. On increasing the value of $\lambda$, the endemic equilibrium bifurcates toward a periodic solution (Theorem 3.2). Numerically, we have shown that system (2.1) is stable globally asymptotically when the consumption rate by the pest is below a threshold value and, after that value, the system is unstable for some higher value of this threshold value giving a stable periodic solution (Figure 2).

In conclusion, in this research, two important dynamical behaviors, namely the global stability of the endemic equilibrium and stability of Hopf bifurcating periodic solution, have been successfully analyzed analytically and numerically using a mathematical model for biological control of crop pests. This article established that the endemic equilibrium is globally stable when the consumption rate of the crop biomass by pests is lower. Also, a Hopf bifurcating periodic solution exits for a higher consumption rate and it is stable. The results will help us in proposing a proper control strategy for pest control in crop cultivation.

Author contributions

J.C., F.A.B.: Conceptualization; J.C., F.A.B.: Methodology; A.A.R., F.A.B.: Software; A.A.R., F.A.B.: Validation; A.A.R., J.C.: Formal analysis; A.A.R, F.A.B.: Investigation; A.A.R., J.C.: Writing-original draft preparation; J.C., F.A.B.: Writing-review and editing; F.A.B.: Supervision; A.A.R.: Project administration. All authors have read and agreed to the published version of the manuscript.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The authors declare no conflict of interest.

Acknowledgments

The authors extend their appreciation to the Deanship of Research and Graduate Studies at King Khalid University, Abha, Kingdom of Saudi Arabia, for funding this work through Large Research Project under grant number RGP2/174/45.