A Mathematical Model to Study the Combined Uses of Infected Pests and Nutrients in Crop Pest Control: Stability Changes and Optimal Control

1. Introduction

Crop pests pose a significant challenge for farmers and agricultural producers, as inadequate management can result in reduced yields, lower-quality crops, or even total crop failure [1]. Common pests that damage crops by feeding on leaves, stems, or fruits include aphids, caterpillars, beetles, and moths. The use of biological pest control agents is gaining popularity in both forests and agriculture [3], driven by environmental concerns and advances in biotechnology, which have encouraged a shift away from chemical pest control methods [2]. A key challenge in future pest management lies in the production and application of biological insect control agents. Natural organisms such as bacteria, fungi, viruses, protozoa, and nematodes are considered effective biological pest control agents [4,5]. Among various biological pest control techniques, the strategy of releasing infected pests into agricultural systems is the most effective. In this context, agricultural awareness plays a vital role.

Agricultural awareness campaigns are essential for informing and educating communities about the damage pests inflict on crops and how to protect them while minimizing pesticide use [6]. Pesticide communication initiatives inform farmers about the risks associated with chemical pesticides, encouraging a reduction in their use [8,19]. Those who understand these risks are more likely to adopt alternative pest control methods, such as Integrated Pest Management (IPM) or biopesticides [20,21]. IPM is a crop protection approach that combines various pest control techniques, selecting the most effective strategies for implementation [22,23]. Its main objective is to limit the use of pesticides and other interventions to economically justifiable levels while minimizing risks to human health and the environment [9]. Farm awareness programs are designed to promote the adoption of biocontrol methods for crop pest management [30]. Enhancing farmer awareness helps farmers prevent crop losses while reducing the negative impacts of pesticide use [7].

Mathematical modeling plays a crucial role in analyzing pest control strategies that use infected pests, typically involving viruses as biopesticides. For example, Pathak et al. [26] studied the dynamics of the pest population influenced by viral infection. The application of entomopathogenic fungi such as Beauveria bassiana and Metarhizium anisopliae within integrated pest management programs has been examined by the authors of [34,35]. Furthermore, Silva et al. [18] emphasized the importance of deploying biocontrol agents at optimal times and locations to maximize pest management efficacy while minimizing the environmental impact. Tan and Chen [36] introduced continuous and impulsive pest management techniques, with an emphasis on scenarios involving the lab release of infected insects.

The incorporation of farming awareness into pest control modeling was first introduced by Basir et al. [10], who examined the dynamics of the pest using biopesticides. Building upon this, Abraha et al. [29] proposed a mathematical model that integrates farming awareness, biopesticide application, and optimal control theory, revealing the occurrence of Hopf bifurcation. Recent studies have explored the use of predatory insects as part of farming awareness-based initiatives to better understand crop–pest dynamics [33].

Nutrients, including both macronutrients and micronutrients, play a vital role in supporting the growth and developmental processes of plants [11]. These nutrients are needed in substantial quantities and are integral to numerous biochemical functions vital to plant development and physiology [13,14]. Effective nitrogen management is key to enhancing crop yields and preserving soil fertility [12]. The appropriate application of macronutrients and micronutrients can strengthen plant disease resistance, improve crop quality, and support sustainable agricultural practices [15]. Therefore, by ensuring the proper nutrient supply, crop production can be enhanced [31], but excess amounts of nutrients can have negative effects on the crop yield [37]. Thus, control of nutrient use is required for maximum crop yields.

The studies mentioned above enhance our understanding of the intricate interactions among pests, pathogens, and various management techniques. The findings offer important guidance for designing pest control methods that are not only efficient but also sustainable and ecologically responsible. However, a mathematical model that integrates farming knowledge-based application of infected pests and nutrients remains unexplored. This study aims to develop a mathematical model that captures the dynamics of crop pest populations influenced by the release of infected pests and the application of nutrients within an agricultural setting. In addition, it examines the role of agricultural awareness in shaping the deployment of control measures, including both the release of infected pests and the input of nutrients. To study cost-effectiveness, an optimal control problem is formulated. The Pontryagin Maximum Principle is applied to identify optimal intervention strategies. Through rigorous mathematical and numerical analyses, our study yields valuable insights into effective and sustainable approaches to crop pest management.

The paper is arranged as follows. In the next section, we propose the mathematical model. The mathematical analysis is presented in Section 3. The model with optimal control approach is analyzed in Section 4. Numerical simulations and sensitivity analysis are presented in Section 5. Section 6 concluded the paper with major findings and real-world applications.

2. The Mathematical Model

We propose an eco-epidemic model for pest control, incorporating four dynamic variables to characterize the system, namely, the biomass of crop, $X\left(t\right)$, the susceptible pest, $S\left(t\right)$, the infected pest, $I\left(t\right)$, and level of agricultural awareness, $A\left(t\right)$, at any time $t>0$. The derivation of the mathematical model is based on the following facts.

H1:

The logistic growth is assumed for crop biomass with r as the intrinsic growth rate and K as the maximum crop biomass [10,27]. In addition, we consider the impact of nutrients on the quick recovery of crops. Let c be the efficacy of the nutrients to enhance the carrying capacity of crop biomass. The modified carrying capacity of crop biomass due to awareness is $K+cA$[31]. The crop population is consumed by both susceptible pests (S) and infected pests (I) at the rates $\alpha$ and $\alpha \varphi$, respectively. Here, $\varphi <1$ as the infected pests are weak and will consume the crops at a lower rate [10]. The dynamics of crop biomass is governed by the following equation:

$${\displaystyle \frac{dX}{dt}}=rX\left(1-{\displaystyle \frac{X}{K+cA}}\right)-\alpha XS-\varphi \alpha XI.$$

(1)

H2:

The susceptible pest population increases due to consumption of crop biomass, with a conversion efficiency ${\theta }_1$. However, it decreases due to natural mortality and infection caused by contact with infected pests [10]. Thus, we have the following equation for the susceptible pest population:

$${\displaystyle \frac{dS}{dt}}={\theta }_1\alpha XS-dS-\lambda SI.$$

(2)

H3:

Infected pests grow through infection of susceptible pests and interactions with the hosts. Also, the infected pests increase due to additional applications influenced by awareness campaigns. Their population is reduced by natural death and infection-induced mortality. Thus, the dynamics of infected pest population is modeled by the equation below:

$${\displaystyle \frac{dI}{dt}}=bIA+\lambda SI+{\theta }_2\varphi \alpha XI-(d+\delta )I.$$

(3)

H4:

The level of agricultural awareness is influenced by various factors. Awareness campaigns increase the level of awareness through media platforms such as social networks and television at a rate $\omega$. The impact of pest damage also increases awareness, directly proportional to the density of the total population of pests, modeled by the term $\sigma (S+I)$[10]. Over time, the level of awareness decreases due to inactivity or lack of involvement at a rate $\mu$. The local awareness rate is the rate at which local awareness is generated or influenced. These factors play a dynamic role in shaping agricultural awareness.

Based on the aforementioned assumptions, we obtain the following mathematical model:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dX}{dt}}& =& rX\left(1-{\displaystyle \frac{X}{K+cA}}\right)-\alpha XS-\varphi \alpha XI,\hfill \\ \hfill {\displaystyle \frac{dS}{dt}}& =& {\theta }_1\alpha XS-dS-\lambda SI,\hfill \\ \hfill {\displaystyle \frac{dI}{dt}}& =& bIA+\lambda SI+{\theta }_2\varphi \alpha XI-(d+\delta )I,\hfill \\ \hfill {\displaystyle \frac{dA}{dt}}& =& \omega +\sigma (S+I)-\mu A.\hfill \end{array}$$

(4)

For solving the above model, the initial values of state variables are taken as follows:

$$X\left(0\right)\ge 0,S\left(0\right)\ge 0,I\left(0\right)\ge 0,A\left(0\right)\ge 0.$$

(5)

Remark 1.

Note that the infected pests are released in the system to make the susceptible pests as they are harmful to crops. The term $bIA$ denotes “awareness-driven release of infected pests”. The term $\lambda SI$ represents the infection rate of susceptible pests when they are in contact with infected pests.

3. Analytical Results

It can be established that the system (4) with the initial condition (5), admits a unique solution [16]. In the following analysis, we identify the equilibrium points of the system and examine their stability properties.

3.1. Existence of Equilibria and Their Stability

We identify five equilibrium points for system (4), as outlined below:

i)

The crop- and pest-free equilibrium, $E_0=\left(0,0,0,{\displaystyle \frac{\omega }{\mu }}\right)$,

ii)

The pest-free steady state, $E_1=(X_1,0,0,A_1)$. Here, $X_1={\displaystyle \frac{K\mu +c\omega }{\mu }}$ and $A_1={\displaystyle \frac{\omega }{\mu }}$.

iii)

The susceptible pest-free equilibrium, $E_2=(\overline{X},0,\overline{I},\overline{A})$, where

$$\begin{array}{ccc}\hfill \overline{X}& =& {\displaystyle \frac{d+\delta -b\overline{A}}{{\theta }_2\varphi \alpha }},\overline{I}={\displaystyle \frac{\mu \overline{A}-\omega }{\sigma }},\hfill \end{array}$$

and $\overline{A}$ is the positive root of

$$C_0{A}^2+C_{1}A+C_2=0,$$

(6)

where, $C_0={\alpha }^{2}c\mu {\varphi }^2{\theta }_2,C_1={\alpha }^{2}k\mu {\varphi }^2{\theta }_2-br\sigma -{\alpha }^{2}c{\varphi }^2{\theta }_2\omega -\alpha c\varphi r\sigma {\theta }_2$, and $C_2=-dr\sigma -\delta r\sigma +\alpha k\varphi r\sigma {\theta }_2+{\alpha }^{2}k{\varphi }^2{\theta }_2\omega .$

Since $C_0>0$, a positive root of (6) exists when $C_2<0.$ Thus the feasibility conditions of $E_2$ are $d+\delta -b\overline{A}>0,\mu \overline{A}-\omega >0$ and $C_2<0$.

iv)

The infected-pest-free equilibrium, $E_3=(X_2,S_2,0,A_2)$, where ${\displaystyle X_2=\frac{d}{\theta \alpha },}$${\displaystyle A_2=\frac{\omega +\sigma S_2}{\mu },}$ and $S_2$ is the positive root of the quadratic equation:

$$\begin{array}{c}\hfill L_0{S}^2+L_{1}S+L_2=0,\end{array}$$

(7)

where, $L_0={\displaystyle \frac{c{\alpha }^2{\theta }_1\sigma }{\mu }},L_1={\displaystyle \frac{c{\alpha }^2{\theta }_1\omega }{\mu }}-k{\alpha }^2{\theta }_1+{\displaystyle \frac{rc\alpha {\theta }_1\sigma }{\mu }},L_2=r\left({\theta }_1\alpha k+{\displaystyle \frac{c\alpha {\theta }_1\omega }{\mu }}-d\right).$ Since $L_0>0$, a positive value of $S_2$ exists when $L_2<0$, i.e., when

$${\theta }_1\alpha \left(k+{\displaystyle \frac{c\omega }{\mu }}\right)<d\Rightarrow {\displaystyle \frac{\omega }{\mu }}<{\displaystyle \frac{d-{\theta }_1\alpha k}{{\theta }_1\alpha c}}.$$

(8)

v)

The coexisting equilibrium, $E^{*}=(X^{*},S^{*},I^{*},A^{*})$, where

$$\begin{array}{ccc}\hfill I^{*}& =& {\displaystyle \frac{A^{*}{\theta }_1\left(b\sigma +\lambda \mu \right)-\lambda {\theta }_1\omega +r{\theta }_2\varphi d-\sigma {\theta }_1\left(d+\delta \right)}{\sigma \lambda \left({\theta }_1-{\theta }_2\varphi \right)}},\hfill \\ \hfill X^{*}& =& {\displaystyle \frac{d}{{\theta }_1\alpha }}+{\displaystyle \frac{A^{*}{\theta }_1\left(b\sigma +\lambda \mu \right)-\lambda {\theta }_1\omega +r{\theta }_2\varphi d-\sigma {\theta }_1\left(d+\delta \right)}{\sigma \alpha {\theta }_1\left({\theta }_1-{\theta }_2\varphi \right)}},\hfill \\ \hfill S^{*}& =& {\displaystyle \frac{(\mu A^{*}-\omega )\lambda \left({\theta }_1-{\theta }_2\varphi \right)-\sigma \{A^{*}{\theta }_1\left(b\sigma +\lambda \mu \right)-\lambda {\theta }_1\omega +r{\theta }_2\varphi d-\sigma {\theta }_1\left(d+\delta \right)\}}{\sigma \lambda \left({\theta }_1-{\theta }_2\varphi \right)}},\hfill \end{array}$$

and $A^{*}$ satisfies the quadratic equation:

$$l_0{A}^2+l_{1}A+l_2=0,$$

(9)

with

$$\begin{array}{ccc}\hfill l_0& =& -{\alpha }^{2}c\lambda \mu \varphi {\theta }_1^2+{\alpha }^{2}bc\sigma {\theta }_1^2-{\alpha }^{2}bc\varphi \sigma {\theta }_1^2+{\alpha }^{2}c\lambda \mu \varphi {\theta }_1{\theta }_2,\hfill \\ \hfill l_1& =& -{\lambda }^2\mu r{\theta }_1-b\lambda r\sigma {\theta }_1-{\alpha }^{2}K\lambda \mu \varphi {\theta }_1^2+{\alpha }^{2}c\lambda \omega \varphi {\theta }_1^2-{\alpha }^{2}cd\sigma {\theta }_1^2\hfill \\ & & -{\alpha }^{2}c\delta \sigma {\theta }_1^2+{\alpha }^{2}bK\sigma {\theta }_1^2+{\alpha }^{2}cd\varphi \sigma {\theta }_1^2+{\alpha }^{2}c\delta \varphi \sigma {\theta }_1^2-{\alpha }^{2}bK\varphi \sigma {\theta }_1^2+\alpha c\lambda r\sigma {\theta }_1^2\hfill \\ & & +{\alpha }^{2}K\lambda \mu \varphi {\theta }_1{\theta }_2-{\alpha }^{2}c\lambda \omega \varphi {\theta }_1{\theta }_2+{\alpha }^{2}cd\varphi r{\theta }_1{\theta }_2-{\alpha }^{2}cd{\varphi }^{2}r{\theta }_1{\theta }_2-\alpha c\lambda \varphi r\sigma {\theta }_1{\theta }_2,\hfill \\ \hfill l_2& =& {\lambda }^2\omega r{\theta }_1+\delta \lambda r\sigma {\theta }_1+{\alpha }^{2}K\lambda \omega \varphi {\theta }_1^2-{\alpha }^{2}dK\sigma {\theta }_1^2-{\alpha }^2\delta K\sigma {\theta }_1^2+{\alpha }^{2}dK\varphi \sigma {\theta }_1^2\hfill \\ & & +{\alpha }^2\delta K\varphi \sigma {\theta }_1^2+\alpha K\lambda r\sigma {\theta }_1^2-d\lambda \varphi r^2{\theta }_2+d\lambda \varphi r\sigma {\theta }_2\hfill \\ & & -{\alpha }^{2}K\lambda \omega \varphi {\theta }_1{\theta }_2+{\alpha }^{2}dK\varphi r{\theta }_1{\theta }_2-{\alpha }^{2}dK{\varphi }^{2}r{\theta }_1{\theta }_2-\alpha K\lambda \varphi r\sigma {\theta }_1{\theta }_2.\hfill \end{array}$$

The existence of positive roots of equation (9) is addressed in the following proposition.

Proposition 1.

(a)

For $l_0>$, $l_1>0$ and $l_2>0$, using Descartes’ rule of sign, there is no positive root of equation (9).

(b)

For $l_0>$ and $l_2<0$, equation (9) admits a positive root.

(c)

For $l_0>$, $l_1<0,l_2<0$ or $l_0<$, $l_1>0,l_2<0$, there exist two positive roots of equation (9) provided that $l_1^2-4l_0{l}_2>0$.

Remark 2.

The susceptible pest-free equilibrium $E_2$ is the most significant for this research. We are controlling the susceptible pest using infected pest release, thus our aim will be to make the system free from susceptible pests.

3.2. Basic Reproduction Number

We define the pest-free equilibrium (DFE) as $E_1(X_1,0,0,A_1)$, where

$$A_1={\displaystyle \frac{\omega }{\mu }},X_1=K+c{A}_1=K+{\displaystyle \frac{c\omega }{\mu }}.$$

We choose the infected subsystem as $(S,I)$ and write it in the next-generation form

$${\displaystyle \frac{d}{dt}}\left(\begin{array}{c}S\\ I\end{array}\right)=F(S,I)-V(S,I),$$

where F collects all new infection terms and V collects all other transition terms.

From the model (4), we have

$$F_1(S,I)={\theta }_1\alpha XS,F_2(S,I)=bIA+\lambda SI+{\theta }_2\varphi \alpha XI,$$

and

$$V_1(S,I)=dS+\lambda SI,V_2(S,I)=(d+\delta )I.$$

Evaluating the Jacobian matrices at the DFE $(X_1,0,0,A_1)$, we obtain

$$F={\left.{\displaystyle \frac{\partial (F_1,F_2)}{\partial (S,I)}}\right|}_{\mathrm{DFE}}=\left(\begin{array}{cc}{\theta }_1\alpha X_1& 0\\ 0& b{A}_1+{\theta }_2\varphi \alpha X_1\end{array}\right),V={\left.{\displaystyle \frac{\partial (V_1,V_2)}{\partial (S,I)}}\right|}_{\mathrm{DFE}}=\left(\begin{array}{cc}d& 0\\ 0& d+\delta \end{array}\right).$$

Thus, the next generation matrix is given by

$$K=F{V}^{-1}=\left(\begin{array}{cc}{\displaystyle \frac{{\theta }_1\alpha X_1}{d}}& 0\\ 0& {\displaystyle \frac{b{A}_1+{\theta }_2\varphi \alpha X_1}{d+\delta }}\end{array}\right).$$

Hence, the basic reproduction number is the spectral radius of K and is obtained as

$$R_0=\rho \left(K\right)=max\left\{{\displaystyle \frac{{\theta }_1\alpha X_1}{d}},{\displaystyle \frac{b{A}_1+{\theta }_2\varphi \alpha X_1}{d+\delta }}\right\}=max\left\{{\displaystyle \frac{{\theta }_1\alpha \left(K+\frac{c\omega }{\mu }\right)}{d}},{\displaystyle \frac{b\frac{\omega }{\mu }+{\theta }_2\varphi \alpha \left(K+\frac{c\omega }{\mu }\right)}{d+\delta }}\right\}.$$

(10)

3.3. Jacobian Matrix and Characteristic Equation

We denote the model variables as $x_i,i=1,2,\dots ,4$. Now, we rewrite our model system (4) as follows:

$${\displaystyle \frac{d{x}_i}{dt}}=f_i(x_1,x_2,\dots ,x_4),i=1,2,\dots ,4.$$

The Jacobian matrix is determined as

$$J\left(x\right)=\left[\begin{array}{cccc}{\displaystyle \frac{\partial f_1}{\partial x_1}}& {\displaystyle \frac{\partial f_1}{\partial x_2}}& \dots & {\displaystyle \frac{\partial f_1}{\partial x_4}}\\ {\displaystyle \frac{\partial f_2}{\partial x_1}}& {\displaystyle \frac{\partial f_2}{\partial x_2}}& \dots & {\displaystyle \frac{\partial f_2}{\partial x_4}}\\ ⋮& ⋮& \ddots & ⋮\\ {\displaystyle \frac{\partial f_4}{\partial x_1}}& {\displaystyle \frac{\partial f_4}{\partial x_2}}& \dots & {\displaystyle \frac{\partial f_4}{\partial x_4}}\end{array}\right].$$

Usually, $J\left(E\right)$ is evaluated at an equilibrium point $E(X,S,I,A)$. The characteristic equation of the Jacobian is obtained by:

$$det\left(J\right(E)-\rho I)=0,$$

(11)

where I is the identity matrix and $\rho$ are the eigenvalues. Sign of the real parts of the eigenvalues determines the stability of the equilibrium point.

3.4. Stability of Equilibria

We now derive the conditions under which an equilibrium is stable. For this, we require determining the nature of the eigenvalues of the Jacobian matrix at that equilibrium point. An equilibrium point is considered asymptotically stable if all eigenvalues have negative real parts. The following theorem addresses the local stability of the previously identified equilibria.

Theorem 1.

i)

The crop- and pest-free steady state, $E_0$ is unstable everywhere.

ii)

The pest-free equilibrium $E_1$ is stable when $R_0<1$ and unstable otherwise.

iii)

The equilibrium point $E_2$ will be stable for the following condition to hold,

$$\begin{array}{c}\hfill {\theta }_1\alpha \overline{X}-d<\lambda \overline{I},\end{array}$$

(12)

and the cubic equation is stable according to the R-H criteria that is ${\tau }_1>0$,${\tau }_2>0$,${\tau }_3>0$ and ${\tau }_1{\tau }_2>{\tau }_3$

iv)

Stability of the equilibrium $E_3$ holds under the following conditions:

$$\begin{array}{c}\hfill {{\rm Y}}_1>0,{{\rm Y}}_4>0,{{\rm Y}}_1{{\rm Y}}_2-{{\rm Y}}_3>0,{{\rm Y}}_1{{\rm Y}}_2{{\rm Y}}_3-{{\rm Y}}_3^2-{{\rm Y}}_1^2{{\rm Y}}_4>0.\end{array}$$

(13)

v)

The following conditions guarantee the stability of the coexisting equilibrium $E^{*}$:

$$\begin{array}{c}\hfill {\gamma }_1>0,{\gamma }_4>0,{\gamma }_1{\gamma }_2-{\gamma }_3>0,{\gamma }_1{\gamma }_2{\gamma }_3-{\gamma }_3^2-{\gamma }_1^2{\gamma }_4>0.\end{array}$$

(14)

Proof. 

i)

The Jacobian matrix at the crop- and pest-free steady state, $E_0(0,0,0,{\displaystyle \frac{\omega }{\mu }})$ is given by

$$J(0,0,0,{\displaystyle \frac{\omega }{\mu }})=\left[\begin{array}{cccc}r& 0& 0& 0\\ \\ 0& -d& 0& 0\\ \\ 0& 0& {\displaystyle \frac{b\omega }{\mu }}-d-\delta & 0\\ \\ 0& \sigma & \sigma & -\mu \end{array}\right].$$

The matrix $J_{E_0}$ gives the following characteristic equation:

$$\begin{array}{c}\hfill (r-\rho )(d+\rho )({\displaystyle \frac{b\omega }{\mu }}-d-\delta -\rho )(\mu +\rho )=0.\end{array}$$

(15)

One eigenvalue of the above matrix is positive, hence the axial equilibrium $E_0(0,0,0,{\displaystyle \frac{\omega }{\mu }})$ is unstable.

ii)

At the pest-free steady state, $E_1(X_1,0,0,A_1)$ the Jacobian matrix takes the form as

$$J(X_1,0,0,A_1)=\left[\begin{array}{cccc}-r& -\alpha X_1& -\varphi \alpha X_1& rc\\ \\ 0& {\theta }_1\alpha X_1-d& 0& 0\\ \\ 0& 0& b{A}_1+{\theta }_2\varphi \alpha X_1-d-\delta & 0\\ \\ 0& \sigma & \sigma & -\mu \end{array}\right].$$

This yields the characteristic equation in $\rho$ as below:

$$\begin{array}{c}\hfill (r+\rho )({\theta }_1\alpha X_1-d-\rho )(b{A}_1+{\theta }_2\varphi \alpha X_1-d-\delta -\rho )(\mu +\eta )=0.\end{array}$$

(16)

Roots of the above equation are obtained as $-r$, $-\mu$, ${\theta }_1\alpha X_1-d$ and $b{A}_1+{\theta }_2\varphi \alpha X_1-d-\delta$. Clearly, all eigenvalues will be negative whenever the following conditions met:

$$\begin{array}{c}\hfill {\theta }_1\alpha (K\mu +c\omega )<d\mu ,b\omega +{\theta }_2\varphi \alpha (K\mu +c\omega )<\mu \left(d+\delta \right).\end{array}$$

(17)

For the stability of equilibrium $E_1$, all the eigenvalues should be negative. Hence, the equilibrium $E_1$ is stable if the conditions in (17) are satisfied. Stability of the equilibrium point $E_1$ requires all eigenvalues to be negative. Accordingly, $E_1$ is stable if the criteria in (17) hold.

iii)

We construct the Jacobian matrix for the susceptible pest-free equilibrium $E_2$ as

$$J_{E_2}=\left[\begin{array}{cccc}{b}_{11}& b_{12}& b_{13}& b_{14}\\ \\ 0& b_{22}& 0& 0\\ \\ b_{31}& b_{32}& b_{33}& b_{34}\\ \\ 0& b_{42}& b_{43}& b_{44}\end{array}\right].$$

where,

$$\begin{array}{ccc}& & b_{11}=r-{\displaystyle \frac{2r\overline{X}}{K+c\overline{A}}}-\varphi \alpha \overline{I},b_{12}=-\alpha \overline{X},b_{13}=-\varphi \alpha \overline{X},b_{14}={\displaystyle \frac{r{\left(\overline{X}\right)}^{2}c}{{(K+c\overline{A})}^2}},\hfill \\ & & b_{22}={\theta }_1\alpha \overline{X}-d-\lambda \overline{I},b_{31}={\theta }_2\varphi \alpha \overline{I},b_{32}=\lambda \overline{I},b_{33}=b\overline{A}+{\theta }_2\varphi \alpha \overline{X}-d-\delta ,\hfill \\ & & b_{34}=b\overline{I},b_{42}=\sigma ,b_{43}=\sigma ,b_{44}=-\mu .\hfill \end{array}$$

Using (11), we determine the characteristic equation of the Jacobian matrix J at $E_2$ as

$$[{\theta }_1\alpha \overline{X}-d-\lambda \overline{I}-{\rho }_2][{\rho }_2^3+{\tau }_1{\rho }_2^2+{\tau }_2\rho +{\tau }_3]=0$$

(18)

$$\begin{array}{ccc}\hfill {\tau }_1& =& -b_{11}-b_{33}-b_{44}\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill {\tau }_2& =& b_{11}(b_{33}+b_{44})+b_{33}{b}_{44}-b_{34}{b}_{43}-b_{13}{b}_{31},\hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill {\tau }_3& =& -b_{11}{b}_{33}{b}_{44}+b_{11}{b}_{34}{b}_{43}+b_{13}{b}_{31}{b}_{44}-b_{14}{b}_{31}{b}_{43}\hfill \end{array}$$

(21)

The equilibrium point $E_2$ is stable if ${\theta }_1\alpha \overline{X}-d<\lambda \overline{I}$ and the cubic equation is stable according to the R-H criteria, that is ${\tau }_1>0$,${\tau }_2>0$,${\tau }_3>0$ and ${\tau }_1{\tau }_2>{\tau }_3$

iv)

We determine the Jacobian matrix at infected pest-free equilibrium $E_3$ to be:

$$J_{E_3}=\left[\begin{array}{cccc}{J}_{11}& J_{12}& J_{13}& J_{14}\\ \\ J_{21}& J_{22}& J_{23}& 0\\ \\ 0& 0& J_{33}& 0\\ \\ 0& J_{42}& J_{43}& J_{44}\end{array}\right].$$

Here,

$$\begin{array}{ccc}& & J_{11}=r-{\displaystyle \frac{2r{X}_2}{K+c{A}_2}}-\alpha S_2,J_{12}=-\alpha X_2,J_{13}=-\varphi \alpha X_2,\hfill \\ & & J_{21}={\theta }_1\alpha S_2,J_{22}={\theta }_1\alpha X_2-d,J_{23}=-\lambda S_2,\hfill \\ & & J_{14}={\displaystyle \frac{r{\left(X_2\right)}^{2}c}{{(K+c{A}_2)}^2}},J_{33}=b{A}_2+\lambda S_2+{\theta }_2\varphi \alpha X_2-d-\delta ,\hfill \\ & & J_{42}=\sigma ,J_{43}=\sigma ,J_{44}=-\mu .\hfill \end{array}$$

The matrix $J_{E_3}$ yields the characteristic equation as below:

$$\begin{array}{c}\hfill {\rho }^4+{{\rm Y}}_1{\rho }^3+{{\rm Y}}_2{\rho }^2+{{\rm Y}}_3\rho +{{\rm Y}}_4=0.\end{array}$$

(22)

Here,

$$\begin{array}{ccc}\hfill {{\rm Y}}_1& =& -(J_{11}+J_{22}+J_{33})+\mu ,\hfill \\ \hfill {{\rm Y}}_2& =& -\mu (J_{11}+J_{22}+J_{33})-J_{12}{J}_{21}+J_{11}(J_{22}+J_{33})+J_{22}{J}_{33},\hfill \\ \hfill {{\rm Y}}_3& =& \mu (J_{11}{J}_{22}-J_{12}{J}_{21}+J_{11}{J}_{33}-J_{22}{J}_{33})+J_{12}{J}_{21}{J}_{33}\hfill \\ & & -J_{11}{J}_{22}{J}_{33}-\sigma J_{21}{J}_{14},\hfill \\ \hfill {{\rm Y}}_4& =& J_{12}{J}_{21}{J}_{33}\mu +J_{11}{J}_{22}\mu \lambda S_2+J_{14}{J}_{21}{J}_{23}\sigma .\hfill \end{array}$$

(23)

Using Routh-Hurwitz criteria, we can conclude that the equilibrium $E_3$ is stable when (13) holds.

v)

At the interior equilibrium point $E^{*}(X^{*},S^{*},I^{*},A^{*})$, the Jacobian matrix corresponding to the system (4) takes the following form:

$$J_{E^{*}}=\left[\begin{array}{cccc}{a}_{11}& a_{12}& a_{13}& a_{14}\\ \\ a_{21}& a_{22}& a_{23}& 0\\ \\ a_{31}& a_{32}& a_{33}& a_{34}\\ \\ 0& a_{42}& a_{43}& a_{44}\end{array}\right].$$

(24)

Here,

$$\begin{array}{ccc}& & a_{11}=r-{\displaystyle \frac{2r{X}^{*}}{K+c{A}^{*}}}-\alpha S^{*}-\varphi \alpha I^{*},a_{12}=-\alpha X^{*},\hfill \\ & & a_{13}=-\varphi \alpha X^{*},a_{14}={\displaystyle \frac{r{\left(X^{*}\right)}^{2}c}{{(K+c{A}^{*})}^2}},a_{21}={\theta }_1\alpha S^{*},\hfill \\ & & a_{22}={\theta }_1\alpha X^{*}-d-\lambda I^{*},a_{23}=\lambda S^{*},a_{31}={\theta }_2\varphi \alpha I^{*},\hfill \\ & & a_{32}=\lambda I^{*},a_{33}=bA+\lambda S^{*}+{\theta }_2\varphi \alpha X^{*}-d-\delta ,\hfill \\ & & a_{34}=b{I}^{*},a_{42}=\sigma ,a_{43}=\sigma ,a_{44}=-\mu .\hfill \end{array}$$

We compute the characteristic equation of the matrix $J_{E^{*}}$ as below:

$$\begin{array}{c}\hfill {\rho }^4+{\gamma }_1{\rho }^3+{\gamma }_2{\rho }^2+{\gamma }_3\rho +{\gamma }_4=0.\end{array}$$

(25)

The coefficients of (22) are as follows:

$$\begin{array}{ccc}\hfill {\gamma }_1& =& -(a_{11}+a_{22}+a_{33}+a_{44}),\hfill \\ \hfill {\gamma }_2& =& a_{44}(a_{11}+a_{22})+a_{33}{a}_{44}-a_{12}{a}_{21}+a_{11}(a_{22}+a_{33})\hfill \\ & & -a_{34}{a}_{43}-a_{23}{a}_{32}+a_{22}{a}_{33}-a_{13}{a}_{31},\hfill \\ \hfill {\gamma }_3& =& a_{44}(a_{12}{a}_{21}-a_{22}{a}_{33})-a_{11}{a}_{44}(a_{22}+a_{33})\hfill \\ & & +a_{12}{a}_{21}{a}_{33}-a_{11}{a}_{22}{a}_{33}+a_{31}(a_{13}{a}_{44}-a_{12}{a}_{23}-a_{21}{a}_{32}+a_{13}{a}_{22})\hfill \\ & & -a_{14}(a_{21}{a}_{42}+a_{31}{a}_{43})+(a_{11}{a}_{43}+a_{22}{a}_{43}-a_{42}{a}_{23})a_{34}\hfill \\ & & +a_{32}(a_{11}{a}_{23}+a_{23}{a}_{44}),\hfill \\ \hfill {\gamma }_4& =& a_{21}(a_{32}{a}_{13}{a}_{44}-a_{12}{a}_{33}{a}_{44})+(a_{11}{a}_{22}{a}_{33}-a_{13}{a}_{22}{a}_{31})a_{44}\hfill \\ & & +a_{14}(a_{21}{a}_{42}{a}_{33}-a_{31}{a}_{42}{a}_{23}+a_{31}{a}_{22}{a}_{43}-a_{21}{a}_{32}{a}_{43})\hfill \\ & & +a_{34}(a_{11}{a}_{42}{a}_{23}-a_{21}{a}_{42}{a}_{13})+a_{43}(a_{12}{a}_{21}-a_{11}{a}_{22})\hfill \\ & & +(a_{31}{a}_{12}{a}_{23}-a_{11}{a}_{23}{a}_{32})a_{44}.\hfill \end{array}$$

(26)

Applying Routh–Hurwitz criterion on equation (22), we get the conditions for roots with negative real parts as follows:

$$\begin{array}{c}\hfill {\gamma }_1>0,{\gamma }_4>0,{\gamma }_1{\gamma }_2-{\gamma }_3>0,{\gamma }_1{\gamma }_2{\gamma }_3-{\gamma }_3^2-{\gamma }_1^2{\gamma }_4>0.\end{array}$$

(27)

Thus the interior equilibrium $E^{*}$ is stable if and only if the conditions in (27) are satisfied.

□

3.5. Existence of Hopf-Bifurcation

The stability condition of coexisting equilibrium is established in (27). Let $\zeta$ be any parameter of the system (4). Now, we derive the conditions for which Hopf bifurcation can occur when $\zeta$ crosses the critical value ${\zeta }^{*}$.

Let $\Psi :(0,\infty )\to R$ be a continuously differentiable function of $\zeta$ such that

$$\Psi \left(\zeta \right):={\gamma }_1\left(\zeta \right){\gamma }_2\left(\zeta \right){\gamma }_3\left(\zeta \right)-{\gamma }_3^2\left(\zeta \right)-{\gamma }_4\left(\zeta \right){\gamma }_1^2\left(\zeta \right).$$

(28)

The following theorem characterizes the occurrence of Hopf bifurcation.

Theorem 2.

The coexisting steady state $E^{*}$ of the system (4) enters into Hopf-bifurcation at $\zeta ={\zeta }^{*}$ if the following conditions hold:

(1)

$\Psi \left({\zeta }^{*}\right)=0$.

(2)

${\gamma }_1^3{\gamma }_2^{\prime }{\gamma }_3({\gamma }_1-3{\gamma }_3)\ne 2({\gamma }_2{\gamma }_1^2-2{\gamma }_3^2)({\gamma }_3^{\prime }{\gamma }_1^2-{\gamma }_1^{\prime }{\gamma }_3^2)$.

Proof. 

Occurrence of Hopf bifurcation at $\zeta ={\zeta }^{*}$ requires the existence of a pair of complex eigenvalues, say $\rho$ and $\overline{\rho }$, such that $Re\rho \left({\zeta }^{*}\right)=0,Im\rho \left({\zeta }^{*}\right)={\omega }_0>0,$ and also the following transversality condition is satisfied:

$${\displaystyle \frac{dRe\rho \left(\zeta \right)}{d\zeta }}{|}_{{\zeta }^{*}}\ne 0.$$

(29)

The remaining roots of the characteristic equation (22) should have negative real parts.

In light of condition $\Psi \left({\zeta }^{*}\right)=0$, the characteristic equation (22) becomes

$$\begin{array}{c}\hfill \left({\rho }^2+{\displaystyle \frac{{\gamma }_3}{{\gamma }_1}}\right)\left({\rho }^2+{\gamma }_1\rho +{\displaystyle \frac{{\gamma }_1{\gamma }_4}{{\gamma }_3}}\right)=0.\end{array}$$

(30)

Let ${\rho }_i$ (i=1,2,3,4) be the four roots of equation (30). Also, we assume that ${\rho }_1$ and ${\rho }_2$ are the pair of purely imaginary roots at $\zeta ={\zeta }^{*}$. Thus, we have ${\rho }_1={\overline{\rho }}_2$ and

$$\begin{array}{c}\hfill \begin{array}{ccc}{\rho }_3+{\rho }_4\hfill & =& -{\gamma }_1,\hfill \\ {\omega }_0^2+{\rho }_3{\rho }_4\hfill & =& {\gamma }_2,\hfill \\ {\omega }_0^2({\rho }_3+{\rho }_4)\hfill & =& -{\gamma }_3,\hfill \\ {\omega }_0^2{\rho }_3{\rho }_4\hfill & =& {\gamma }_4.\hfill \end{array}\end{array}$$

(31)

Here, ${\omega }_0=\mathrm{Im}{\rho }_1\left(\zeta \right)$ implies that ${\omega }_0=\sqrt{{\displaystyle \frac{{\gamma }_3}{{\gamma }_1}}}$.

When ${\rho }_3$ and ${\rho }_4$ appear as complex conjugates, equation (31) implies that $2\mathrm{Re}{\rho }_3=-{\gamma }_1$. Alternatively, if the pair consists of real values, then by referencing both (22) and (31), it can be deduced that ${\rho }_3<0$ and ${\rho }_4<0$. To complete the analysis, the next step is to verify the transversality condition.

Since $\Psi \left({\zeta }^{*}\right)$ depends continuously on its roots, an open interval $\zeta \in ({\zeta }^{*}-ϵ,{\zeta }^{*}+ϵ)$ can be identified in which ${\rho }_1$ and ${\rho }_2$ remain complex conjugates. Within this neighborhood, let their typical representations be given by

$$\begin{array}{c}\hfill \begin{array}{ccc}{\rho }_1\left(\zeta \right)\hfill & =& \varphi \left(\zeta \right)+i\nu \left(\zeta \right),\mathrm{and}{\rho }_2\left(\zeta \right)=\varphi \left(\zeta \right)-i\nu \left(\zeta \right).\hfill \end{array}\end{array}$$

Next, we proceed to examine the following transversality condition:

$${\displaystyle \frac{dRe\left({\rho }_j\left(\zeta \right)\right)}{d\zeta }}{|}_{\zeta ={\zeta }^{*}}\ne 0,j=1,2.$$

Substituting ${\rho }_j\left(\zeta \right)=\varphi \left(\zeta \right)\pm i\nu \left(\zeta \right)$ into (22) and then taking the derivative, we have

$$\begin{array}{c}\hfill \begin{array}{c}K\left(\zeta \right){\varphi }^{\prime }\left(\zeta \right)-P\left(\zeta \right){\nu }^{\prime }\left(\zeta \right)+M\left(\zeta \right)=0,\hfill \\ P\left(\zeta \right){\varphi }^{\prime }\left(\zeta \right)+K\left(\zeta \right){\nu }^{\prime }\left(\zeta \right)+Q\left(\zeta \right)=0.\hfill \end{array}\end{array}$$

Here,

$$\begin{array}{c}\hfill \begin{array}{ccc}K\left(\zeta \right)\hfill & =& 4{\varphi }^3-12\varphi {\nu }^2+3{\gamma }_1({\varphi }^2-{\nu }^2)+2{\gamma }_2\varphi +{\gamma }_3,\hfill \\ P\left(\zeta \right)\hfill & =& 12{\varphi }^2\nu +6{\gamma }_1\varphi \nu -4{\varphi }^3+2{\gamma }_2\varphi ,\hfill \\ M\left(\zeta \right)\hfill & =& {\gamma }_1{\varphi }^3-3{\gamma }_1^{\prime }\varphi {\nu }^2+{\gamma }_2^{\prime }({\varphi }^2-{\nu }^2)+{\gamma }_3^{\prime }\varphi ,\hfill \\ Q\left(\zeta \right)\hfill & =& 3{\gamma }_1^{\prime }{\varphi }^2\nu -{\gamma }_1^{\prime }{\nu }^3+2{\gamma }_2^{\prime }\varphi \nu +{\gamma }_3^{\prime }\varphi .\hfill \end{array}\end{array}$$

Solving for ${\varphi }^{\prime }\left({\zeta }^{*}\right)$, we have

$$\begin{array}{c}\hfill \begin{array}{ccc}{\displaystyle \frac{dRe\left[{\rho }_j\left(\zeta \right)\right]}{d\zeta }}{|}_{\zeta ={\zeta }^{*}}={\varphi }^{\prime }{\left(\zeta \right)}_{\zeta ={\zeta }^{*}}\hfill & =& -{\displaystyle \frac{P\left({\zeta }^{*}\right)Q\left({\zeta }^{*}\right)+K\left({\zeta }^{*}\right)M\left({\zeta }^{*}\right)}{K^2\left({\zeta }^{*}\right)+L^2\left({\zeta }^{*}\right)}}\hfill \\ & =& {\displaystyle \frac{{\gamma }_1^3{\gamma }_2^{\prime }{\gamma }_3({\gamma }_1-3{\gamma }_3)-2({\gamma }_2{\gamma }_1^2-2{\gamma }_3^2)({\gamma }_3^{\prime }{\gamma }_1^2-{\gamma }_1^{\prime }{\gamma }_3^2)}{{\gamma }_1^4{({\gamma }_1-3{\gamma }_3)}^2+4{({\gamma }_2{\gamma }_1^2-2{\gamma }_3^2)}^2}}\ne 0,\hfill \end{array}\end{array}$$

if ${\gamma }_1^3{\gamma }_2^{\prime }{\gamma }_3({\gamma }_1-3{\gamma }_3)\ne 2({\gamma }_2{\gamma }_1^2-2{\gamma }_3^2)({\gamma }_3^{\prime }{\gamma }_1^2-{\gamma }_1^{\prime }{\gamma }_3^2)$. Thus, the transversality condition for Hopf bifurcation is verified. Consequently, the occurrence of Hopf bifurcation at $\zeta ={\zeta }^{*}$ is confirmed. □

4. The Optimal Control Problem

To determine the optimal strategy for releasing infected pests, we introduce an optimal control problem. Specifically, we define two control variables: $c_1\left(t\right)$ (to regulate the release rate of infected pests) and $c_2\left(t\right)$ (to regulate the awareness campaigns), which govern the timing and intensity of pest deployment. These controls are assumed to be admissible over the interval $[t_0,t_f]$, where $t_0$ and $t_f$ denote the initial and final times of intervention, respectively. We introduce the parameter $b_0$ for the control-efficacy which will be regulated by $c_1\left(t\right)$. We also introduce one another control function $c_3\left(t\right)$ that regulates the application of nutrients. The admissibility conditions are given by $0<c_1\left(t\right)<1$ and $0<c_2\left(t\right)<1$, ensuring both practicality and biological feasibility. Our aim is to determine the optimal controls $c_1^{*}\left(t\right)$ and $c_2^{*}\left(t\right)$ that minimize a given cost functional, achieved by applying Pontryagin’s Maximum Principle [17]. Incorporating these controls leads to the following form of the associated state system:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dX}{dt}}& =& rX\left(1-{\displaystyle \frac{X}{K+c_3\left(t\right)cA}}\right)-\alpha XS-\varphi \alpha XI,\hfill \\ \hfill {\displaystyle \frac{dS}{dt}}& =& {\theta }_1\alpha XS-dS-\lambda SI,\hfill \\ \hfill {\displaystyle \frac{dI}{dt}}& =& (1-c_1\left(t\right))b_0+bIA+\lambda SI+{\theta }_2\varphi \alpha XI-(d+\delta )I,\hfill \\ \hfill {\displaystyle \frac{dA}{dt}}& =& (1-c_2\left(t\right))\omega +\sigma (S+I)-\mu A\hfill \end{array}$$

(32)

with initial conditions

$$X\left(0\right)=X_0,S\left(0\right)=S_0,I\left(0\right)=I_0,A\left(0\right)=A_0.$$

(33)

The objective function for the minimization problem is formulated as follows:

$$\begin{array}{ccc}\hfill J(c_1,c_2)& =& {\int }_{t_0}^{t_f}[D{c}_1^2\left(t\right)+E{c}_2^2\left(t\right)+F{c}_3^2\left(t\right)-PX\left(t\right)+Q{S}^2\left(t\right)]dt.\hfill \end{array}$$

(34)

Our aim in proposing the optimal control problem is to maximize crop yield while minimizing pest damage. The objective function (34) will capture the overall cost of pest management. Specifically, $E\ge 0$ represents the weighting constant associated with the susceptible pests, while $D\ge 0$ denotes the weighting constant linked to the release of infected pests. The term $D{c}_1^2\left(t\right)$ accounts for the cost incurred in introducing infected pests into the crop system. The goal is to determine the optimal controls, $c_i^{*}\left(t\right),i=1,2$, which satisfy the prescribed conditions and minimize the associated costs. Mathematically, we have

$$\begin{array}{c}\hfill J(c_1^{*},c_2^{*})=min(J(c_1,c_2):c_1,c_2\in U)\end{array}$$

(35)

subject to the state system (32), where

$$\begin{array}{ccc}\hfill U& =& \{c_1\left(t\right),c_2\left(t\right):c_j\left(t\right)ismeasurableand0\le c_j\le 1,t\in [t_0,t_f],j=1,2\}.\hfill \end{array}$$

4.1. Characteristic of the Optimal Control Triplet

According to the Pontryagin Minimum Principle [17], the optimal control problem is governed by a set of necessary conditions. This principle reformulates the optimality conditions into the task of minimizing the Hamiltonian function, $H\left(t\right)$. The Hamiltonian $H\left(t\right)$ is defined as follows:

$$\begin{array}{ccc}\hfill H& =& D{c}_1^2\left(t\right)+E{c}_2^2\left(t\right)+F{c}_3^2\left(t\right)-PX\left(t\right)+Q{S}^2\left(t\right)\hfill \\ & & +{\eta }_1\left(rX\left(1-{\displaystyle \frac{X}{K+c_3\left(t\right)cA}}\right)-\alpha XS-\varphi \alpha XI\right)\hfill \\ & & +{\eta }_2\left({\theta }_1\alpha XS-dS-\lambda SI\right)\hfill \\ & & +{\eta }_3\left(c_1\left(t\right)b_0+bIA+\lambda SI+{\theta }_2\varphi \alpha XI-(d+\delta )I\right)\hfill \\ & & +{\eta }_4\left((1-c_2)\omega +\sigma (S+I)-\mu A\right).\hfill \end{array}$$

(36)

As stated in [17], the optimal controls $c_1^{*}\left(t\right)$ and $c_2^{*}\left(t\right)$ satisfy

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial H\left(t\right)}{\partial c_1^{*}}}=0,{\displaystyle \frac{\partial H\left(t\right)}{\partial c_2^{*}}}=0.\end{array}$$

(37)

Now,

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial H}{\partial c_1}}& =& 2D{c}_1\left(t\right)+{\eta }_3\left(b_{0}IA\right),\hfill \\ \hfill {\displaystyle \frac{\partial H}{\partial c_2}}& =& 2E{c}_2^{*}\left(t\right)-{\eta }_4\omega .\hfill \end{array}$$

Setting ${\displaystyle \frac{\partial H}{\partial c_i}=0}$ ($i=1,2$) gives the optimal controls as

$$c_1^{*}\left(t\right)={\displaystyle \frac{b_0}{2D}}{\eta }_3\left(t\right),c_2^{*}\left(t\right)={\displaystyle \frac{{\eta }_4\omega }{2E}}.$$

In view of the boundedness of the controls, we have

$$\begin{array}{ccc}\hfill & & c_1^{*}\left(t\right)=min\left\{max\left\{0,\frac{b_0}{2D}{\eta }_3\left(t\right)\right\},1\right\},\hfill \\ & & c_2^{*}\left(t\right)=max\left\{0,min\left\{1,{\displaystyle \frac{{\eta }_4\omega }{2E}}\right\}\right\}.\hfill \end{array}$$

(38)

Control $c_3$ is determined as follows:

$${\displaystyle \frac{\partial \mathcal{H}}{\partial c_3}}=2F{c}_3+{\lambda }_X·{\displaystyle \frac{rcA{X}^2}{{\left(K+c_{3}cA\right)}^2}}\Rightarrow 2F{c}_3^{*}\left(t\right)+{\lambda }_X\left(t\right){\displaystyle \frac{rcA\left(t\right)X^2\left(t\right)}{{\left(K+c_3^{*}\left(t\right)cA\left(t\right)\right)}^2}}=0.$$

Thus, $c_3^{*}\left(t\right)$ is the unique solution of the above equation, projected onto $[0,1]$. We numerically find $c_3$ using the iteration method (Newton method).

4.2. The Adjoint System

In accordance with the Pontryagin Minimum Principle [17], the adjoint variables satisfy the following relations:

$$\begin{array}{c}\hfill {\displaystyle \frac{d{\eta }_i}{dt}}=-{\displaystyle \frac{\partial H}{\partial x_i}},i=1,2,3,4,\end{array}$$

(39)

where $x_i=(X,S,I,A)$ i.e., $x_1=X$, $x_2=S$, etc. Thus, we have

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\eta }_1}{dt}}& =P-{\eta }_1\left[r-{\displaystyle \frac{2rX}{K+c_{3}cA}}-\alpha S-\varphi \alpha I\right]-{\eta }_2{\theta }_1\alpha S-{\eta }_3{\theta }_2\varphi \alpha I,\hfill \\ \hfill {\displaystyle \frac{{\eta }_2}{dt}}& =-2QS-\left[{\eta }_1(-\alpha X)+{\eta }_2({\theta }_1\alpha X-d-\eta I)+{\eta }_3\left(\eta I\right)+{\eta }_4\sigma \right],\hfill \\ \hfill {\displaystyle \frac{{\eta }_3}{dt}}& =-\left[{\eta }_1(-\varphi \alpha X)+{\eta }_2(-\eta S)+{\eta }_3\left(bA+\eta S+{\theta }_2\varphi \alpha X-(d+\delta )\right)+{\eta }_4\sigma \right],\hfill \\ \hfill {\displaystyle \frac{{\eta }_4}{dt}}& =-{\eta }_1{\displaystyle \frac{r{c}_{3}c{X}^2}{{(K+c_{3}cA)}^2}}-{\eta }_{3}bI+\mu {\eta }_4.\hfill \end{array}$$

(40)

The preceding equations delineate the necessary conditions that both the optimal control functions and the state variables are required to satisfy. According to the Pontryagin Maximum Principle [17], the adjoint variables fulfill the terminal boundary conditions ${\eta }_i\left(t_f\right)=0$ for $i=1-4.$

In the following theorem, we characterize the optimal system.

Theorem 3.

Assuming the objective function $J(c_1,c_2,c_2)$ attains its minimum over the admissible set U at the optimal control pair $(c_1^{*},c_2^{*},c_3^{*})$, associated with the coexisting equilibrium point $E^{*}=(X^{*},S^{*},I^{*},A^{*})$, there exist adjoint variables ${\eta }_1,{\eta }_2,{\eta }_3$ and ${\eta }_4$ that satisfy the system of differential equations (40). The adjoint variables satisfy the boundary conditions: ${\eta }_l\left(t_f\right)=0,(l=1,2,3,4).$ The optimal controls can be expressed through the following functions:

$$\begin{array}{ccc}\hfill & & c_1^{*}\left(t\right)=max\left\{0,min\left\{1,{\displaystyle \frac{{\eta }_3{b}_0}{2D}}\right\}\right\},\hfill \\ & & c_2^{*}\left(t\right)=max\left\{0,min\left\{1,{\displaystyle \frac{{\eta }_4\omega }{2E}}\right\}\right\},\hfill \end{array}$$

(41)

and $c_3^{*}\left(t\right)$ is the unique solution of the following equation:

$$2F{c}_3^{*}\left(t\right)+{\lambda }_X\left(t\right){\displaystyle \frac{rcA\left(t\right)X^2\left(t\right)}{{\left(K+c_3^{*}\left(t\right)cA\left(t\right)\right)}^2}}=0,$$

projected onto $[0,1]$.

Remark 3.

The optimal solution is obtained using the forward–backward sweep method. In this approach, the state system (32) is first integrated forward in time with initial values (33) for the control variables. Subsequently, the adjoint system (40) is integrated backward in time. The control variables in (38) are then updated according to the optimality conditions derived from the problem. For the control $c_3\left(t\right)$, we apply the Newton iteration method, which will work along with forward–backward sweep iterative scheme. The processes of forward integration, backward integration, and control updating are repeated iteratively until convergence is achieved, thereby ensuring that the solution satisfies both the state and adjoint systems as well as the optimality conditions.

5. Numerical Simulations

In this section, we present some numerical results that support the analytical findings discussed earlier. Further, by varying key parameters, we examine the system’s behavior under different conditions. Unless otherwise mentioned in the text, the parameter values used for numerical results are same as in Table 1

5.1. Local Sensitivity for $R_0$

We derive the basic reproduction number as follows:

$$R_0=max\left\{R_{0,1},R_{0,2}\right\}$$

with

$$R_{0,1}={\displaystyle \frac{{\theta }_1\alpha \left(K+\frac{c\omega }{\mu }\right)}{d}},R_{0,2}={\displaystyle \frac{b\frac{\omega }{\mu }+{\theta }_2\varphi \alpha \left(K+\frac{c\omega }{\mu }\right)}{d+\delta }}.$$

We compute the sensitivity indices using the following set of parameter values:

$$\begin{array}{cc}\hfill & r=0.1,K=500,\alpha =0.0004,\varphi =0.5,{\theta }_1=0.6,{\theta }_2=0.59,\\ & d=0.1,\lambda =0.0025,b=0.002,\delta =0.01,\omega =0.05,\mu =0.12,c=0.2\end{array}$$

(42)

with

$$A_1={\displaystyle \frac{\omega }{\mu }}={\displaystyle \frac{0.05}{0.12}}\approx 0.4167,X_1=K+c{A}_1=500+0.2·0.4167\approx 500.0833.$$

For this set of parameter values, we have

$$R_{0,1}\approx 1.20,R_{0,2}\approx 5.37.$$

Hence, $R_0=R_{0,2}\approx 5.37$. For a parameter p, the normalized sensitivity is defined as

$$E_p={\displaystyle \frac{p}{R_{0,2}}}{\displaystyle \frac{\partial R_{0,2}}{\partial p}}.$$

Since $R_0=R_{0,2}$, we determine the sensitivity indices as follows:

$$\begin{array}{cccccc}\hfill E_b& ={\displaystyle \frac{b{A}_1}{N}},\hfill & \hfill E_{{\theta }_2}& ={\displaystyle \frac{{\theta }_2\varphi \alpha X_1}{N}},\hfill & \hfill E_{\varphi }& ={\displaystyle \frac{{\theta }_2\varphi \alpha X_1}{N}},\hfill \\ \hfill E_{\alpha }& ={\displaystyle \frac{{\theta }_2\varphi \alpha X_1}{N}},\hfill & \hfill E_K& ={\displaystyle \frac{K}{X_1}}·{\displaystyle \frac{{\theta }_2\varphi \alpha X_1}{N}},\hfill & \hfill E_c& ={\displaystyle \frac{c\omega }{\mu X_1}}·{\displaystyle \frac{{\theta }_2\varphi \alpha X_1}{N}},\hfill \\ \hfill E_{\omega }& ={\displaystyle \frac{\omega }{\mu }}·{\displaystyle \frac{b+c{\theta }_2\varphi \alpha }{N}},\hfill & \hfill E_{\mu }& =-{\displaystyle \frac{\omega }{\mu }}·{\displaystyle \frac{b+c{\theta }_2\varphi \alpha }{N}},\hfill & \hfill E_d& =-{\displaystyle \frac{d}{d+\delta }},\hfill \\ \hfill E_{\delta }& =-{\displaystyle \frac{\delta }{d+\delta }}.\hfill \end{array}$$

Substituting the values of parameters, we plot the sensitivity indices in Figure 1. From this analysis, we conclude that the most influential parameters are ${\theta }_2$, $\varphi$, $\alpha$ and K, each with elasticity close to $+1$, meaning a 10% increase in any of them raises the value of $R_0$ by about 10%. The natural death rate d has a strong negative effect ($E_d\approx -0.91$), while extra mortality rate due to infection ($\delta$) has a moderate negative effect ($E_{\delta }\approx -0.09$). Also, we note that the parameters $b,c,\omega$ and $\mu$ have less effect on pest-free system.

5.2. Global Sensitivity Analysis

To evaluate the impact of parameter uncertainty on crop production, a global sensitivity analysis based on Partial Rank Correlation Coefficients (PRCCs) is conducted. PRCC is a widely used technique for nonlinear dynamical systems, as it accounts for parameter interdependencies while quantifying their relative influence on the model outputs [38,39,40]. The parameter space is efficiently sampled using Latin Hypercube Sampling. The PRCC results for the crop biomass X are shown in Figure 2 & Figure 3. The carrying capacity (K) and intrinsic growth rate (r) are positively correlated with the crop biomass, whereas the crop consumption rate ($\alpha$) and the conversion efficiency of susceptible pests (${\theta }_1$) exhibit negative correlations. Parameters associated with infection dynamics, including $\lambda$, ${\theta }_2$ and d, show moderate positive correlations, reflecting their indirect role in reducing pest pressure. The remaining parameters display relatively small PRCC values, indicating a limited direct influence on the crop biomass. This outcome is consistent with Figure 17, which demonstrates that awareness primarily enhance crop production through indirect regulation of pest populations.

5.3. Dynamics Without Optimal Control

Figure 4 illustrates the stability of the pest-free equilibrium $E_1$, achieved under lower values of the infection rate ($\lambda$) and the rate of introduction of infected pests (b). Meanwhile, Figure 5 demonstrates that the infected pest-free equilibrium $E_2$ exists and remains stable when the consumption rate $\alpha$ is sufficiently high. Figure 6 shows the stability of susceptible pest-free equilibrium $E_2$. Solution trajectories converge to the point $E_2$ for higher values of b. Effects of the parameter b on the pest population is shown in Figure 7. Notably, the susceptible pests can be eradicated by the release of infected pests.

Figure 8 illustrates that all populations coexist within the system for different values of the infection rate $\lambda$. For each value of $\lambda$, the system exhibits damped oscillatory behavior before approaching a steady state, indicating that the system remains stable within this parameter range. The transition from stability to sustained oscillations is further confirmed in Figure 9, which shows that a Hopf bifurcation occurs at $\lambda ={\lambda }^{*}$. Again, we plot the bifurcation diagram with respect to b in Figure 10. The occurrence of Hopf bifurcation at the critical values are confirmed in Figure 11. Figure 12 highlights the influence of the rate of crop biomass consumption by pests on the system dynamics. As the consumption rate increases, the system transitions into oscillatory behavior, indicating that a higher consumption rate can induce periodic oscillations through a Hopf bifurcation. Figure 13 confirms the occurrence of a Hopf bifurcation at the critical value $\alpha ={\alpha }^{*}\approx 0.00045$.

Impact of nutrient application is shown in Figure 14. The figure shows that nutrients are useful for quick recovery from pest damage. Additionally, Figure 15 presents a Hopf bifurcation diagram with respect to c, representing the rate of introduction of infected pests. The results suggest that a higher introduction rate is unfavorable for pest control, as it destabilizes the system and increases pest management costs. The impact of local and global farming awareness campaigns $\left(\omega \right)$ on crop production is shown in Figure 16. The figure demonstrates that these control strategies can effectively reduce pest populations and enhance crop biomass, leading to improved crop yields.

5.4. Results from the Optimal Control Problem

The optimal control problem represented by (32), (36) and (38), is formulated as a two-point boundary value problem. The forward–backward sweep method treats the state system as an initial value problem while solving the adjoint system (40) as a boundary value problem. The state system (32) is solved using forward iteration whereas the system (38) is solved through the backward iteration schemes, using the parameter values from Table 1 except:

$$D=0.5,E=1.0,F=0.5,P=0.005,Q=0.2.$$

(43)

We solve the optimal control problem numerically, and the results are presented graphically in Figs. Figure 17 & Figure 18. These visual representations illustrate the significant impact of the implemented control strategies on the system’s dynamics, demonstrating their effectiveness in optimally suppressing the pest population. The results highlight that a high level of control effort is initially required to achieve substantial pest reduction. However, as time progresses, the necessity of control effort gradually decreases. This strategic approach not only minimizes the costs associated with control implementation but also enhances agricultural productivity. These findings emphasize the dual advantages of cost efficiency and improved yield achieved through optimal control strategies.

5.5. Pest Reduction Ratio, Crop Yield Increase and Control Cost Reduction Rate

The following quantitative metrics are evaluated to compare the baseline (no control) and optimal control scenarios.

(i)

Pest Reduction Ratio (PRR): It is defined as follows:

$$\mathrm{PRR}={\displaystyle \frac{{\int }_{t_0}^{t_f}{S}_{\mathrm{base}}dt-{\int }_{t_0}^{t_f}{S}_{\mathrm{opt}}dt}{{\int }_{t_0}^{t_f}{S}_{\mathrm{base}}dt}}$$

(ii)

Crop Yield Increase (CYI): It is defined as follows:

$$\mathrm{CYI}={\displaystyle \frac{{\int }_{t_0}^{t_f}{X}_{\mathrm{opt}}dt-{\int }_{t_0}^{t_f}{X}_{\mathrm{base}}dt}{{\int }_{t_0}^{t_f}{X}_{\mathrm{base}}dt}}$$

(iii)

Control Cost Reduction Rate (CCRR): It is defined as follows:

$$\mathrm{CCRR}={\displaystyle \frac{J_{\mathrm{base}}-J_{\mathrm{opt}}}{J_{\mathrm{base}}}}$$

(iv)

Awareness Efficiency (AE): Yield gain per unit awareness cost.

(v)

Infected Release Efficiency (IRE): Pest reduction per unit infected-release cost.

Using the parameter values given in Table 1 and (42), we implement the forward–backward sweep method to compute the optimal control trajectories $c_1^{*}\left(t\right)$ and $c_2^{*}\left(t\right)$ (as shown in Figure 18) and also compute the results summarized in Table = Table 2.

The results indicate that the optimal control strategy achieves approximately 47% reduction in the susceptible pest population, 22% increase in the crop yield, and 28% reduction in the overall cost compared to the baseline. The efficiency metrics further highlight that awareness campaigns and the release of infected pests contribute significantly to balancing yield improvement and pest suppression.

6. Conclusions

The escalating challenge of managing pest populations poses a significant threat to agricultural productivity, necessitating sustainable and environmentally responsible control strategies. Biological control strategies, which involve the strategic utilization of natural enemies such as predators, parasitoids, and pathogens, offer considerable advantages over conventional chemical pesticides by mitigating environmental contamination and reducing risks to human health. A key approach in biological pest control involves the introduction of infected pests carrying lethal pathogens, which effectively suppresses pest populations while minimizing chemical interventions. By promoting awareness and informed decision-making, the effectiveness of biological control measures can be significantly enhanced. In this study, we investigated a mathematical model that employs natural predators as a biologically sustainable alternative to chemical pesticides for pest regulation. Our model incorporates the interactions between pest populations, their natural predators, and nutrient dynamics within the agricultural ecosystem. Additionally, the model accounts for the influence of agricultural awareness among farmers, highlighting its critical role in optimizing pest management strategies. Furthermore, our study integrates an optimal control framework to evaluate the efficacy of biologically driven pest suppression strategies. The results provide valuable insights into the role of biological control in reducing pest densities, minimizing economic costs, and improving crop yields. The findings underscore the potential of integrating biological control with strategic management interventions to achieve sustainable agricultural practices while mitigating the detrimental effects of pests on crop productivity.