MATHEMATICAL MODELLING THAT CHARACTERIZES THE QUALITATIVE BEHAVIOR OF THE INTERACTIONS OF APHIDS AND LADYBUGS AS INSECTIVOROUS ANIMALS

. The signiﬁcance of the agricultural industry is primarily driven by the growth of the human population and the imperative for sustainable food production. Within this sector, the presence of insect pests, such as aphids, demands particular attention, as they exert a detrimental inﬂuence on agricultural practices. To devise effective approaches, it is essential to possess a thorough comprehension of aphid control. Thus, this study presents a mathematical model that captures the qualitative dynamics of interactions between aphids and ladybugs. The model is composed of a system of non-linear differential equations, encompassing the interplay between the two species, the impact of insecticides on aphids, and the hibernation stages of ladybugs. The analytical approach is employed to ascertain the criteria for the existence and stability of the equilibrium derived from the constructed model. The study focused on stability and sensitivity analysis to understand the behavior of the model and develop effective


INTRODUCTION
The agricultural industry is crucial to current global civilization and growth, driven by the growing human population and their ongoing need for food.This sector is a vital component of economic development and a significant source of income for nations.Despite advancements in technology, the significance of agriculture remains unchanged [1].However, any disruption or decline in agricultural production can negatively impact economic growth.One factor contributing to a decrease in crop yields is the presence of pest insects, such as aphids.Aphids can harm plants in various ways and impact the agricultural sector.They feed on the phloem tissue of the plant and alter its nutrient flow, leading to stunted growth and decreased reproduction.Additionally, the saliva they inject while feeding can be toxic to the plant.Aphids are also vectors for many viruses, with nearly 50% of all viruses transmitted by insects.Of the 45 main insect pests that attack the top 6 food crops (corn, wheat, potatoes, sweet potatoes, wheat, and tomatoes), 26% are caused by aphids, resulting in yearly direct losses [2].This leads to the growth of a black filamentous saprophytic ascomycetes fungus, which frequently hinders photosynthetic activity, an essential aspect in nourishing the plants [3].Aphids are a significant concern as they can spread plant viruses, transmitting almost 30% of all plant virus species.
They have a sucking mouthpart structure, composed primarily of a stylet bundle, which can spread disease [4,5].
Ladybugs (Harmonia axyridis) are a natural predator of aphids on agricultural crops as they feed on aphids [6].They also feed on other hemipterans like Adelgidae, Psyllidae, and Aleyrodidae [7].The adult ladybug has the ability to consume 15-65 aphids per day, and the larvae have the potential to eat 90-370 aphids during their development stage.This makes them effective biological control agents for insect pests in agriculture [8].The hunting behavior of ladybugs is triggered by chemical signals released by damaged plants due to aphid infestation.These signals, along with the honeydew excreted by aphids and alarm pheromones, act as additional chemical clues that help the ladybugs pinpoint their prey [9].Ladybugs are also considered effective biological control agents against other invasive pests due to their emission of odors that serve as self-defense.When they are disturbed by a pheromone response, they reflexively release noxious compounds as a defense mechanism [10].The smell produced by the compound is potent and has a unique odor, reminiscent of freshly sliced green bell peppers or green beans [11,12].Ladybugs belong to the family of insects known as Coccinellidae, which includes various species that contain alkaloid toxins (such as adaline, coccinelline, exochomine, and hippodamine) that serve as defensive mechanisms.Additionally, some adult ladybugs exhibit aposematic coloration, which acts as a warning signal to predators.Consequently, the ladybugs are usually not preyed upon by predators due to their toxic and unpalatable nature [13,7].
Aside from utilizing natural predators, a variety of insecticides are commonly employed for aphid management.This approach is evident among smallholder farmers in Limpopo Province who heavily depend on synthetic insecticides to manage insect pests on their vegetables [14].
Due to their prolonged efficacy, positive impact on crop yield, and marginal rate of return, cypermethrin and acetamiprid insecticides are recommended for aphid control.The South African government's Department of Agriculture, Forestry and Fisheries (DAFF) has stated that more than 3000 pesticide products, which are synthetic in nature, have been authorized for the purpose of managing insect pests by local farmers [15].The utilization of pesticides in China has risen from 733 million kg in 1990 to 1.806 billion kg in 2012, and over 2.2 billion kg of pesticides are presently manufactured in China, a significant proportion of which is exported.China has around 2000 pesticide manufacturers [16].Smallholder farmers' frequent reliance on synthetic insecticides has resulted in several issues, including the contamination of the environment, the emergence of pest resistance, and human health problems [17].A thorough examination of 122 studies that were published after the year 2000 indicated that certain groups, particularly farmers and pesticide applicators, had a considerably higher likelihood of experiencing certain types of cancer, dementia, and respiratory symptoms [18].Athukorala et al. [19] found that in Sri Lanka, Costa Rica, and Nicaragua, 4% − 7% of agricultural workers experience health problems caused by pesticides every year.The accessibility and price of synthetic insecticides are still significant obstacles, particularly in distant rural regions where agriculture is the primary source of livelihood and food supply [20].Furthermore, numerous studies have indicated that the utilization of insecticides can lead to a resistance issue [21,22], leading to elevated expenses associated with their future usage.Nonetheless, the utilization of insecticides persists as a viable choice due to their efficacy in managing aphids.This was validated through a field experiment carried out by the Agricultural Research Institute-Tarnab Peshawar, which yielded the finding that the application of synthetic insecticides led to a significant reduction in the aphid population when compared to natural insecticides [23].
Many mathematicians have examined the interaction behavior between aphids and ladybugs by constructing mathematical models in an effort to comprehend it better.The study in [24] utilized a mathematical model to investigate the dynamics between planthopper pests as prey and Menochilus sexmaculatus and mirid ladybug as two predators, with the control of the prey being maintained through the application of pesticides.The proposed singular system model of aphid ecosystems takes into account the impact on aphid populations due to changes in the parameter related to their natural enemy population, within the context of the fold catastrophe manifold [25].The soybean aphids can influence the physiology of soybean mainly through two means, by promoting feeding and overcoming resistance, which results in further colonization by its own species.A differential equation population model with non-local features was established in [26] to examine the impact of virulent and avirulent aphids on soybean plants and the dynamics of the biological processes involved.In [27], a new mathematical model consisting of differential equations is presented to study the seasonal population dynamics of aphids, ladybugs, and ants.The model was initially presented in continuous time, then transformed into a discrete-time representation using the Nonstandard Finite Difference scheme (NSFD) method.It features four non-hyperbolic steady states.The study conducted by Anjos et al., [28] established a Bayesion framework to analyze aphid-ladybeetle interaction and dynamics by analyzing six ordinary differential equations that include three phenomenological models and three data-driven models.
In this study, we construct a mathematical model to investigate the interactions between aphid and ladybug populations, taking into account the utilization of insecticides.Furthermore, we conducted a detailed examination of ladybug population growth as insectivores, incorporating hibernation stages.While our model remains straightforward, the conducted studies encompass a comprehensive range, including equilibrium stability analysis, sensitivity analysis, and exploration of the optimal control problem.Derived from this analysis, our research will add to the current body of literature concerning the growth dynamics of aphid populations, specifically focusing on the impacts of insecticides and insectivores on these populations.

MATHEMATICAL FORMULATION
A mathematical model that describes the interaction between predators and prey is used in studying aphid control using insectivores, in this case, the ladybug population.Suppose X(t) is the population of aphids and Y (t) is the population of ladybugs at time t.The logistic function is believed to control the growth rates of both populations as a result of internal competition within each population.Due to the use of insecticides, the aphids population will decrease by uX where u is the control parameter which represents the proportion of deaths due to insecticide exposure.The ladybug population is thought to be unharmed by the insecticide being used in this particular case.Biologically, the ladybug population is an insectivore for the aphids population.The relationship can be described as a predator-prey interaction, with ladybugs taking on the role of predator and aphids serving as the prey.However, the aphids population is very sensitive to stink carried by the ladybug population.Consequently, aphids are able to rapidly escape when a ladybug population approaches them.This condition is very reasonable if represented by a response function of Holling type II for aphids when interacting with ladybug populations.Mathematically, the rate of the aphids population will decrease by (mXY )/(η + ξ X) where m is the death rate due to interaction with the ladybug population and the ability of the aphids to escape from the ladybugs is represented by parameter ξ .Meanwhile, the ladybug population received (εmXY )/(η + ξ X) biomass due to interactions with aphids.In addition, this model also captures conditions where the ladybug population hibernates.Let Z(t) represent the hibernating ladybug population at time t.The rate of the ladybug population decreases by αY and increases by γZ due to the hibernation factor where α is the transition rate to the hibernating ladybug population and γ is the reverse transition rate.As a result, the rate of the hibernating ladybug population increases by αY and decreases by γZ.Mathematically, the model is represented by equation ( 1) as follows: Table 1 displays the parameters for the system (1).

MATHEMATICAL ANALYSIS
3.1.Equilibrium and their stability.We have four equilibria, three equilibria in boundedness and one interior equilibrium.
E 0 = (0, 0, 0), represents the extinction of all population represents the extinction of aphids

the extinction of the ladybug population
To investigate the stability of the three equilibria, the Jacobian matrix of system (1) is evaluated one by one at each equilibrium.From the Jacobian matrix evaluated at equilibrium E 0 yield the polynomial characteristics, The first term of the polynomial equation (3) reveals that an eigenvalue will be negative when r < u.However, achieving stability at the equilibrium point E 0 is not possible due to the second term of the polynomial (3), which results in two eigenvalues with different signs.This discrepancy is evident from the negative sign of the polynomial constant.Consequently, based on this analysis, it can be concluded that the equilibrium point E 0 is unstable regardless of the conditions.
Furthermore, the stability of the equilibrium E 1 is determined by the following characteristic polynomial equation obtained from the Jacobian matrix which has been evaluated at equilibrium Part A of the polynomial p 1 produces two eigenvalues which are negative because all the coefficients of the polynomial are positive.Therefore, the stability of the equilibrium E 1 is determined by the equation of part B from the polynomial p 1 .Part B of the polynomial p 1 produces an eigenvalue λ = −µ + (r − u) which is negative if (r − u) < µ which correspond to the stability of the equilibrium E 1 and positive if (r − u) > µ which correspond to the instability of the equilibrium E 1 .The stability from the equilibrium E 1 is preserved when u = 0 if the condition r < µ is satisfied.Meanwhile, from the Jacobian matrix evaluated at equilibrium E 2 , three eigenvalues are obtained.One explicit eigenvalue, namely λ = −(r − u) which in itself is always negative because of the conditions that must be met for the existence of an equilibrium E 2 , and the other two eigenvalues are the roots of the following characteristic polynomial equation, It can be obviously verified that the roots of the polynomial p 2 have different signs because the coefficient of the second degree is always negative while the constant of the polynomial is positive.Therefore, under its conditions of existence, the equilibrium E 2 is always unstable.
Furthermore, the existence of an interior equilibrium that represents the coexistence between aphids and ladybug populations will be explored.The existence of an interior equilibrium is determined by the third-degree polynomial equation in (5) which is the remainder of the elimination results from the system (1) for the variable x, where and x * is determined by the polynomial (5) β c 2 r From y * we can know that x * is in the interval 0 < x * < 1 − u/r provided that r > u.If r < u then the existence of x * is on the negative axis and that is impossible from a biological perspective.
Under r > u + µ we have that A 2 and A 1 are positive, but A 0 is negative.This situation indicates that the change in the sign of the coefficients of the polynomial (5) occurs once.According to Descartes' rule of signs [29], polynomials have one positive root x * .The requirements for existence typically correspond to the conditions ensuring equilibrium stability.In this instance, we will verify the stability through numerical methods simulation.The outcomes of this analysis enable us to visualize the parameter range that characterizes both the existence and stability of the equilibrium, as depicted in Figure 1.
Parameter space indicating the region of existence and equilibrium stability Region I represents a stable domain for the E 1 equilibrium, indicating that when the aphids growth rate is lower than the death rate resulting from interactions with ladybugs, the aphid populations will eventually perish.Aphid populations can sustain only if their growth rate surpasses the death rate, whether through insecticide exposure, interactions with ladybugs, or a combination of both.This scenario is referred to as having the "rate in" more significant than the "rate out".Region II represents a shift resulting from insecticide usage, revealing that such application can reduce the coexistence area (Region III).The extent of insecticide use correlates with the contraction of the interior equilibrium stability region.Additionally, insecticide application expands the size of the E 1 equilibrium stability region, signifying a decline in aphid populations.These findings align with the earlier interpretation.
Region III represents an area where the interior equilibrium point E 3 is stable while the E 1 equilibrium is unstable.This implies that coexistence becomes possible with reduced insecticide usage.However, if the objective is to eradicate aphids, the analysis results in the diagram indicate that both the death rate due to insecticide application and the predation by ladybugs must be intensified.Specifically, this involves either using higher doses or increasing the frequency of insecticide application for the former.As for the latter, enhancing mortality through interactions with ladybugs can be achieved by augmenting the population of ladybugs in agricultural regions.

NUMERICAL SIMULATION AND BIOLOGICAL INTERPRETATION
4.1.Orbit solution.The numerical solution of system (2) will be accomplished using the Runge-Kutta method to validate the analysis outcomes.In this section, we will verify the stability of two equilibrium points: firstly, the extinction point of the aphid population, E 1 , and secondly, the stability of the interior equilibrium E 3 , which was implicitly demonstrated in the analysis.The simulation results for both stability conditions are depicted in Figure 2.
In the preceding analysis, it was established that the condition where the "rate in" exceeds the "rate out" in the aphid population (r > u + µ) ensures the presence of an interior equilibrium.The numerical findings further validate the existence of a stable interior equilibrium, as illustrated in Figures 2(a  it is asserted that the aphid population will increase if insecticides are not used.Notably, the ladybug population exhibits a coexistence scenario with the aphids due to the "rate in" being greater than the "rate out."However, when the "rate out" surpasses the "rate in" in the aphid population, the interior equilibrium E 3 vanishes, leading to a stable E 1 equilibrium and signifying the extinction of the aphid population.This understanding provides valuable insights for controlling aphids as plant pests.Before delving into the ladybug population control strategy, it is essential to examine the effects of altering each parameter on the compartments.Hence, the following section will address the sensitivity of parameter changes.

Sensitivity analysis.
In this section, the sensitivity of the parameters to changes in the behavior of the solution of system (2) will be investigated.The approach employed for this purpose resembles the methodology utilized in various articles, as demonstrated in [30,22,31].
Interesting parameters to observe for sensitivity include parameters r, u, µ, c, β , and α.Suppose a vector of variables is V = (x, y, z), a vector of parameters is Φ = {α, β , r, u, c, µ}, and a vector of the right-hand side of the system (2) is F. Define sensitivity is S = ∂V ∂ Φ .By total differential, we have: where J(V ) is the 3 × 3 Jacobian matrix of the system (2), S and The numerical simulation results using the Runge-Kutta method for the system of differential equations ( 6) which are calculated based on the parameter data when the interior equilibrium is stable are shown in Figure 3. Figure 3 is organized into two rows and three columns, where each column visually represents the sensitivity level of the parameters investigated for the variable.In the situation of equilibrium, the parameters u and µ have the most significant impact on the alterations in compartment x for which the effect of their modification is inversely proportional to the changes in the compartment (Figure 3 (a)).This implies that higher dosages result in a decreased aphid population.Similarly, an increase in the number of ladybug populations or a rise in their aggressiveness will result in a decline in the aphid population.Meanwhile, the parameters that exhibit the highest sensitivity to changes in the solution variable y are the parameters u and r (Figure 3 (e)).While insecticide exposure does not directly impact the aphid population, it still leads to a decrease in the ladybug population.Additionally, the growth rate of aphids, represented by parameter r, also influences the growth rate of the ladybug population (Figure 3 (e)).This comprehensive description outlines the complete model.Although the birth rate of ladybugs indeed positively impacts their population, its visibility is not significant under equilibrium conditions (Figure 3 (b)).Similarly, the hibernation rate negatively affects variable y, but its influence is not prominently evident in equilibrium conditions (Figure 3 (b)).Instead, the hibernating ladybug population is primarily influenced by the hibernation transition rate itself, which crucially determines the population size during hibernation (Figure 3 (c)).A higher rate results in a larger hibernating ladybug population.Moreover, the impact of insecticide exposure on hibernating ladybug populations is still observed, albeit at half the magnitude compared to non-hibernating ladybug populations (Figure 3 (f)).

4.3.
Characterization of the optimal control problem.Here, we aim to reduce the aphid population by presenting multiple scenarios in this model.The main idea of the interventions is to reduce the population of aphid as much as possible, but with as low as possible of intervention of insecticides.Hence, we threat u as a time-dependent parameter u(t).To achieve this, we establish the objective function of the optimal control problem as described in equation ( 9). ( 9) where ζ is the weight parameter of u(t).This ζ parameter should balance each term of the function that we want to minimize, or in the other hand ζ should satisfy 2x(t) For the sake of written simplification, we write u(t) as u, and other variables also.To characterize the optimal control problem, lets define the Hamiltonian function of our problem as follows: Using the Pontryagin Minimum Principle [32,33,34], we taking the negative partial derivative of H with respect to each state variables yield: 11) where the transversality condition λ i (T f ) = 0 for i = 1..3.The optimal control for using insecticide effect u is obtained by solving dH du = 0 respect to u, which gives us: With a specific upper and lower bound of u, we have the optimal value of u is given by ( 12) u opt = max u min , min u max , λ 1 x ζ 4.4.Optimal control simulation.In this section, we conduct numerical simulations for the optimal control problem.To solve our optimal control problem, we have a state system described by the model in system (2), with non-negative initial conditions given, an adjoint system in (11) with given transversality conditions, and the optimal control dynamics given by u opt in (12).To solve this problem, we use a forward-backward sweep method [35,36].This method has been widely used by many authors [37,38,39,40,41].
We start by providing an initial guess of u(t), let's call it u 0 , which remains constant for t ∈ [0, T f ].Hence, we can solve the model in (2) forward in time to determine x(t), y(t), and z(t) for all t ∈ [0, T f ].Next, we solve the adjoint system in (11) backward in time, with a transversality condition λ i (T f ) = 0 for i = 1, 2, 3. Consequently, we can update u(t) using the formula in (12).
We repeat this procedure until the convergence criterion is achieved (||J i+1 − J i || < 10 −5 , with i being the i-th iteration), or it reaches the maximum iteration limit (n = 100).It is clear to see that the intervention of insecticides successfully reduces the proportion of aphids, ladybugs, and hibernating ladybugs.The dynamics of controls start increasing at the beginning of the simulation time to reduce the outbreak of x.After t = 15, the intensity of control starts to decrease and reaches zero at the final time of the simulation.Although insecticides succeed in significantly reducing x(t), the dynamics of y(t) and z(t) do not show a huge impact, which indicates that the population of ladybugs is relatively stable and persistent due to the disturbance in the population of aphids.This phenomenon occurs because the recruitment of the ladybug population is not solely dependent on predation to aphids but also relies on other resources described by the logistic equation for the recruitment of the ladybug population.To describe our claim, we run an optimal control simulation where ladybug only depend on predation to aphids to survive.The results is given in Figure 5 where β = 0, while other parameters remain the same as in Figure 4.   Conversely, when we use five different initial conditions for the ladybug population, as shown in Figure 7, fewer ladybugs at t = 0 require less intensity of insecticide to be effective.Hence, the cost function when y(0) = 0.1 is the cheapest compared to when y(0) is larger than 0.1.The cost functions for the cases with y(0) equal to 0.1, 0.2, 0.3, 0.4, and 0.5 are 0.0045, 0.0041, 0.0036, 0.0030, and 0.0022, respectively.

CONCLUSION
A mathematical model was developed to depict the interplay between aphids and ladybugs, incorporating the usage of insecticides.Explicitly, the model generates three equilibria, signifying the extinction of all populations, the extinction of aphids, and the extinction of the ladybugs population.Among these explicit equilibria, only one is stable, specifically the equilibrium corresponding to the extinction of aphids.In the meantime, the presence of the interior equilibrium is deduced indirectly by analytically determining its conditions for existence.Subsequently, its stability is verified through numerical assessment, utilizing parameter values that fulfill the conditions for existence.Optimal control simulations are conducted to find the best solutions in various possible scenarios.From the simulation of the energy conversion factor's impact, it is found that the presence of active ladybugs in predating aphids can assist the success of insecticides in controlling the aphids population.Furthermore, it is also discovered that high ) and 2(b).

Figure 2 (
a) represents the scenario when u = 0, while Figure 2 (b) portrays the situation when u = 0. Observing the change in the equilibrium position from Figure 2(a) to Figure 2(b), a shift to the right on the x-axis becomes evident.Biologically,

TABLE 1 .
(1)cription of the parameters and their values.Before analyzing model(1), normalization of the model is carried out first to simplify the observation interval.The model is normalized to the carrying capacity of each population.
ξ the proportion of aphids that can escape ladybugs non-dimensional β growth rate of ladybug populations time −1 η half-saturation constant number ε coefficient of energy conversion for ladybug from predation activity non-dimensional