As a major pest in the global citrus industry, Toxoptera citricida pose a serious threat to citrus production. Their populations can experience explosive growth under suitable conditions, and if left uncontrolled, they can cause significant economic losses in a short period. However, the traditional chemical pesticide control model has fallen into a dilemma, while biological control methods suffer from serious response lag.

To break through the existing prevention and control bottlenecks, we have developed an intelligent and precise pesticide intervention system based on RNAi technology and established a comprehensive modeling framework with multi-scale, multi-stage, and closed-loop feedback. This framework integrates population ecology, pharmacokinetics, transmission dynamics, and multi-objective optimization theory, achieving full-chain mathematical modeling from molecular mechanisms to field applications.

The entire modeling system consists of one background module and three core functional modules:

Climatic conditions, especially temperature, have a decisive impact on aphid population dynamics. Through literature research and expert interviews, we found that there is a certain relationship between temperature and the intrinsic rate of increase of aphids. The study by Wang and Tsai^[1] provided us with the key mathematical expression:

$r_{\text{total}}(T)=e^{0.0727T}-e^{0.0727\times 37.2517-\frac{37.2517-T}{11.8862}}-0.5053$

This non-linear function accurately describes the reproductive potential of aphids at different temperatures. To verify the applicability of this model in Guangdong region, we collected the daily average temperature data for the whole year of 2024 in Guangdong Province and overlaid it on the theoretical growth rate curve for comparative analysis.

Fig 1. Relationship between temperature and aphid intrinsic growth rate, validated with 2024 daily temperature data from Guangdong Province.

Fig 2 Annual temperature dynamics and aphid growth rate in Guangdong Province during 2024.

Analysis reveals that aphids reach their optimal reproductive state (intrinsic rate of increase r ≈ 0.30) around 24.5°C, and in 2024, approximately 70% of the time in Guangdong Province falls within the suitable temperature range of 15 - 30°C, creating favorable conditions for the continuous outbreak of aphids. The annual temperature dynamics exhibit a bimodal distribution, with temperatures approaching the optimum during two periods: spring (April - June) and autumn (September - November), when aphid populations grow the fastest, while low temperatures in winter (December - February of the following year) and extreme high temperatures in summer (>30°C) naturally inhibit their reproduction. This finding clearly points out the critical prevention and control window period - preventive interventions should be implemented at the initial stage of the spring outbreak and before the autumn reproductive peak.

● Fundamentals of Ecology: S-shaped Curve

The aphid population cannot grow indefinitely, and its dynamics are subject to triple constraints:

1. Natural resource constraints: Limited supply of citrus phloem sap
2. Intraspecific competition pressure: Spatial competition for feeding check points
3. Environmental carrying capacity constraints: Constraints on plant health status

Based on the classic Logistic model of Verhulst and Pearl, combined with temperature driving factors, we established the basic growth equation:

${\displaystyle \frac{dN}{dt}}=r_{\text{total}}(T(t))\cdot N\cdot (1-{\displaystyle \frac{N}{K}})$

Where N is the aphid population size, K is the environmental carrying capacity, $r_{\text{total}}(T(t))$ is the temperature-dependent intrinsic growth rate.

This equation clearly depicts the two phases of population growth: exponential growth when N << K, and a gradual decline in the growth rate to zero as N approaches K.

● Ecological Facts: Predation and Control by Natural Enemies

To accurately describe the regulatory role of natural enemies, we introduce the Holling II functional response, which can precisely characterize the "predatory saturation phenomenon"—when aphid density increases, the predation rate of natural enemies does not increase linearly but gradually approaches saturation.

The mathematical expression of the predation term is:

${\displaystyle \frac{aN}{1+ahN}}\cdot P$

Where P is the population size of natural enemies, a is the attack rate (the ability of natural enemies to search for and discover prey per unit time), and h is the handling time (the average time required to capture, kill, and digest an aphid).

● Ecological Fact: Winged Aphid Dispersal

Aphids possess unique phenotypic plasticity: the same genotype can develop into wingless forms (suitable for rapid reproduction) or winged forms (suitable for migration and dispersal). This is a sophisticated survival strategy. The review by Braendle et al.^[2] confirms that population density is the most critical signal triggering wing form transformation. When resources become scarce and the population is overcrowded, the mother will produce more winged offspring and initiate the migration program.

Although there is currently a lack of direct quantitative research on the proportion of winged aphids and the relative density of the population, based on the ecological principle of density-dependent regulation, we propose the S-shaped threshold response hypothesis:

1. Low density period (N/K < 0.5): Winged aphids are maintained at a baseline level (~5%)
2. Threshold period (N/K ≈ 0.7): Wing form differentiation rises sharply
3. High-density period (N/K > 0.9): The proportion of winged gnats approaches the upper limit (~50%)

Figure 3 shows the proposed wing shape conversion model:

$\phi (x)=0.05+\frac{0.45}{1+e^{-10(x-0.7)}}$

Among them, $x=\frac{N}{K}$ is the relative density, and this curve has the following key characteristics:

• Baseline level stability: When the relative density < 0.5, the proportion of winged aphids remains at a low level of approximately 5%, ensuring that the population prioritizes investment in reproduction rather than dispersal when resources are abundant
• Quick Response Threshold: Near x = 0.7 (the red marked point), the slope of the curve is the largest, and at this time the proportion of winged aphids is approximately 27.5%. The ecological significance of this threshold lies in that when the population reaches 70% of the carrying capacity, the resource competition pressure begins to rise significantly, triggering the "early warning mechanism".
• Saturation Upper Limit: When the relative density is > 0.9, the proportion of winged aphids approaches the upper limit of 50%, which reflects the physiological cost constraint of wing morph differentiation - even under extreme crowding conditions, the population will not completely shift to the winged morph

This "early warning mechanism" enables the population to start preparing for migration before resources are completely depleted, which is crucial for long-term survival. This hypothesis is supported by the observational data from our laboratory.

● Local dispersal dynamics of winged aphids

The winged individuals of the Toxoptera citricida do not possess long-distance migratory ability, and their dispersal is mainly achieved through short-distance flight between adjacent plants. In a two-dimensional grid orchard, each citrus tree occupies a lattice point (i,j), and its set of neighbors is defined as:

$N_{i,j}=(i-1,j),(i+1,j),(i,j-1),(i,j+1)$

The migration flow from plant (i,j) to neighbor (k,l) is determined by 4 factors:

1. Diffusion ability: Diffusion coefficient m
2. Emigration source strength: Only winged aphids can disperse, and the number of winged aphids on plant (i,j) is $W_{i,j}(t)=\phi \left(\frac{N_{i,j}}{K_{i,j}}\right)\cdot N_{i,j}$
3. Density-driven pressure: The higher the density, the easier it is to migrate, i.e., pressure is proportional to relative density, $\sigma \left(\frac{N_{i,j}}{K_{i,j}}\right)$
4. Target Attractiveness: $1-\frac{N_{k,l}}{K_{k,l}}$

The complete expression for the migration flow is:

$J_{(i,j)\to (k,l)}=m\cdot W_{i,j}(t)\cdot \sigma \left(\frac{N_{i,j}}{K_{i,j}}\right)\cdot (1-\frac{N_{k,l}}{K_{k,l}})$

Considering that diffusion is bidirectional, the net diffusion amount of plant (i,j) is:

$D_{i,j}(t)={\sum }_{(k,l)\in {𝒩}_{i,j}}[J_{(k,l)\to (i,j)}-J_{(i,j)\to (k,l)}]$

● Comprehensive Spatiotemporal Dynamics Equation

By integrating all the above ecosystem processes, we obtain the complete dynamic equation for aphids in the two-dimensional orchard space:

${\displaystyle \frac{d{N}_{i,j}}{dt}}=r_{\text{total}}(T(t))\cdot N_{i,j}\cdot (1-{\displaystyle \frac{N_{i,j}}{K_{i,j}}})-{\displaystyle \frac{a{N}_{i,j}}{1+ah{N}_{i,j}}}\cdot P_{i,j}+D_{i,j}(t)$

This equation clearly decomposes the three drivers of population change:

• Item 1: Logistic growth of temperature regulation
• Item 2: Population loss caused by natural enemy predation
• Item 3: Winged aphid-mediated spatial redistribution

The population dynamics of natural enemies (mainly predatory insects such as ladybugs and lacewings) are driven by three key processes:

1. Predation Transformation Growth: Energy is obtained through preying on aphids, and the transformation coefficient e represents how many offspring of natural enemies can be supported by each captured aphid
2. Natural death: Mortality dp includes factors such as aging, disease, etc.
3. Active tracking and diffusion: Natural enemies do not diffuse randomly, but actively move towards areas with high aphid density (chemotaxis)

The predator flow from plant (i,j) to neighbor (k,l) is expressed as:

$J_{(i,j)\to (k,l)}^P=m_P\cdot P_{i,j}(t)\cdot \frac{N_{k,l}(t)}{N_{i,j}(t)+ϵ}$

Among them, $m_P$ is the diffusion coefficient of natural enemies, $\frac{N_{k,l}(t)}{N_{i,j}(t)+ϵ}$ which reflects the following effect (ϵ is a very small value to avoid the denominator being zero)

The final form of the predator-prey equation is:

$\frac{d{P}_{i,j}}{dt}=e\cdot \frac{a{N}_{i,j}}{1+ah{N}_{i,j}}\cdot P_{i,j}-d_p{P}_{i,j}+D_{i,j}^P(t)$

Among them, $D_{i,j}^P(t)$ is the net diffusion of natural enemies, and the calculation method is similar to that of aphids.

After communicating with experts in the field of ecosystem, we decided to make the model closer to the real orchard environment. In actual orchards, factors such as soil fertility, moisture, light, and microclimate are spatially unevenly distributed, and this spatial heterogeneity significantly affects the distribution pattern of aphids.

We define environmental carrying capacity as a function of spatial location $K_{i,j}$, whose value is jointly modulated by the baseline carrying capacity and three effects:

$K_{i,j}=K_{\text{base}}\cdot G(i,j)\cdot B(i,j)\cdot H(i,j)$

1. Gradient Effect

Simulated continuous layer changes in factors such as sunlight, water source, and soil fertility in the orchard. For example, the upper left corner has high terrain and good sunlight, while the lower right corner may be close to a water source but has poor sunlight.

$G(i,j)=1-\gamma \cdot (\frac{i}{M}+\frac{j}{N})$

Where $\gamma \in \left(\mathrm{0,1}\right)$ is the gradient decay intensity coefficient. $\frac{i}{M}$ and $\frac{j}{N}$ normalize the coordinates to the $\left[\mathrm{0,1}\right]$ interval.

2. Boundary Effect

Simulated that plants at the orchard edge may be more exposed to wind and extreme temperatures, resulting in poor health conditions.

Assume that the boundary width is b grids, and the boundary bearing capacity decays to β times that of the interior.

$B(i,j)=\{\begin{array}{cc}\beta ,& \text{if }i<b\text{ or }i>M-b-1\text{ or }j<b\text{ or }j>N-b-1\\ 1,& \text{otherwise}\end{array}$

3. Hotspot Effect

Simulated the patchy distribution of resources caused by random factors such as uneven fertilization, localized waterlogging, and particularly robust individual plants.

Randomly generate S hotspot centers, whose influence decays Gaussianly with distance.

$H(i,j)=1+\alpha \cdot {\sum }_{s=1}^S\exp (-\frac{{(i-I_s)}^2+{(j-J_s)}^2}{2{\sigma }^2})$

Where $(I_s,J_s)$ is the center coordinate of the s-th hotspot.$\alpha$ is the maximum value of hotspot intensity. $\begin{array}{l}\\ sigma\end{array}$ controls the range of hotspot influence.

When analyzing the control threshold of natural enemies against aphids, we discovered a counterintuitive fact: Although it is known that the daily average predation of ladybugs, one of the natural enemies of aphids, is 56 - 150 individuals^[3], we cannot simply infer that the control threshold is 1:56 to 1:150. The reason is that:

1. Both populations are in dynamic change
2. Factors such as temperature and density affect the growth rate in real time
3. Spatial distribution is uneven

Without considering migration (since the migration term does not change the total number of aphids in the entire system, but only changes the spatial distribution), the aphid population equation is:

${\displaystyle \frac{dN}{dt}}=r_{\text{total}}(T(t))\cdot N\cdot (1-{\displaystyle \frac{N}{K}})-{\displaystyle \frac{aN}{1+ahN}}\cdot P$

The control threshold should at least prevent the aphid population from growing, so the critical value condition is: ${\displaystyle \frac{dN}{dt}}=0$

This means:

$r_{\text{total}}(T(t))\cdot N\cdot (1-{\displaystyle \frac{N}{K}})={\displaystyle \frac{aN}{1+ahN}}\cdot P$

After transformation, we obtain dynamic control threshold:

${\left(\frac{P}{N}\right)}_{\text{crit}}=\frac{r_{\text{total}}(T)\cdot (1+a\cdot h\cdot N)\cdot (K-N)}{a\cdot K}$

Further transformation reveals that this formula well explains the reason why the threshold varies with the number of aphids:

${\left(\frac{P}{N}\right)}_{\text{crit}}=\{r_{\text{total}}(T)\}\cdot \{(1+a\cdot h\cdot N)\}\cdot (1-\frac{N}{K})\cdot \left\{\frac{1}{a}\right\}$

1. The higher the temperature, the higher the proportion of natural enemies required (aphids grow rapidly)
2. The higher the aphid density, the more obvious the predation saturation effect, and the more difficult the control
3. When approaching the carrying capacity, density dependence reduces the required proportion of natural enemies
4. The higher the natural enemy attack rate, the lower the required proportion

It can be seen that at 25°C (a common temperature in Guangdong), when the aphid density reaches 600 aphids per plant, the required natural enemy ratio is as high as 1:2, that is 1 natural enemy is needed for every 2 aphids. This high-cost, high-tech threshold operation is not realistic for ordinary fruit farmers. This finding clearly points out the inherent limitations of traditional biological control methods, while chemical control also has obvious limitations, which is also the core motivation for our development of the intelligent Toxoptera citricida management project APHiGO. Learn more in Description.

We used the RK45 algorithm to simulate the aphid-natural enemy dynamics over 20 days in a 5×5 grid orchard, based on the real temperature data of Guangdong in 2024.

Through comprehensive analysis of Figures 5 and 6, we identified the critical time window for aphid control. Figure 6A shows that the aphid population was in the growth stage from Day 0 to Day 5, during which the upper row of Figure 5 shows that aphids were only concentrated in 1-2 core areas, while Figure 6B shows that the proportion of winged aphids was only 5-10%, with relatively low dispersal pressure. However, after Day 5, a critical point, the proportion of winged aphids soared to 47%, the population skyrocketed to 27,000 within 5 days, and the spatial distribution evolved from a point-like core to orchard-wide diffusion. Comparing the natural enemy dynamics in the middle row of Figure 5, it can be seen that the response of natural enemies lagged by 5-7 days and their numbers were severely insufficient, far below the 1:2 dynamic control threshold shown in Figure 4.

This finding reveals the inherent flaws of traditional biological control and provides a scientific basis for our two-stage strategy: During the growth period (0-5 days), trunk injection of RNAi pesticides is used for targeted and precise strikes, and only 20-30% of the affected area needs to be treated to prevent over 95% of aphids from entering the outbreak stage; once this golden window is missed and the outbreak stage is entered, an emergency plan of orchard-wide spraying with RNAi-Metarhizium must be initiated, at which point the prevention and control cost will increase by 5-8 times. Therefore, early warning based on model recognition and phased intervention are the keys to achieving low-cost and high-efficiency aphid management. Learn more in Design.

Although the Background model reveals the intrinsic laws of aphid outbreaks, the severity of climate change and the inherent limitations of the model make field prediction inevitably uncertain. Traditional monitoring methods rely on manual inspections or sticky trap counts, which are not only time-consuming and labor-intensive but also suffer from severe spatio-temporal sampling biases—by the time aphids are detected, the optimal control window has often been missed.

To break through this bottleneck, we designed an innovative "aphid-microorganism-gas" three-level linkage real-time monitoring system. The core idea is to construct a complete biological signal transduction chain:

Aphids feed on citrus → excrete honeydew (rich in sucrose) → engineered Bacillus subtilis specifically recognizes and metabolizes sucrose → produces volatile signaling molecule methyl salicylate (MeSA) → gas sensor detects in real time → mathematical model back-calculates aphid population size

The unique advantage of this system lies in automated continuous monitoring, high sensitivity, and early warning. Learn more in Hardware.

We constructed a multi-level mechanistic model that decomposes complex biological processes layer by layer into quantifiable mathematical equations.

Layer 1: Aphid honeydew production

Aphids obtain large amounts of sap from the phloem through their piercing-sucking mouthparts, but only absorb the amino acids and some sugars in it, with the remaining carbohydrates excreted in the form of honeydew. Based on the normalization method of Moir ML et al.^[5], honeydew production is related to the body weight of the insects by a power function:

$H(t)=\alpha \cdot N(t)\cdot W^{\beta }$

Where α = 0.144 μg/h, β = 0.3156, and W is the average body weight of the worm

Layer 2: Total Sucrose in Honeydew

Honeydew is not a pure sucrose solution, but a complex mixture containing sugars, amino acids, and organic acids. Considering the physicochemical properties of honeydew:

$S_{\text{total}}(t)=H(t)\cdot \rho \cdot C_{\text{sucrose}}$

Where $\rho$ is the density of honeydew, $C_{\text{sucrose}}$ is the concentration of sucrose in honeydew

Layer 3: Sucrose Consumption by Bacillus subtilis

The engineered Bacillus subtilis, as the core component of the "biosensor", follows the classic Michaelis-Menten enzymatic reaction kinetics in its sucrose metabolism:

$S_{\text{consumed}}(t)={\mu }_{\text{max}}{B}_{\text{ss}}\cdot \frac{S_{\text{total}}(t)}{K_{\text{s}}+S_{\text{total}}(t)}$

Where ${\mu }_{\text{max}}$ is the maximum specific growth rate, $K_s$ is the Michaelis constant, $B_{\text{ss}}$ is the steady-state biomass

Layer 4: Generation and Volatilization of Methyl Salicylate

Product Generation: Based on metabolic flux balance, the stoichiometric relationship between sucrose consumption and MeSA production is:

$P(t)=Y_{\text{P/S}}\cdot S_{\text{consumed}}(t)$

Volatilization Kinetics: As a volatile organic compound (VOC), the mass transfer of MeSA at the gas-liquid interface follows Henry's Law and the Two-Film Theory:

$V(t)=k_{v}AP(t)\cdot \frac{H_c}{RT}$

Dynamic Equilibrium of Gas Phase Concentration: The gas phase concentration is affected by various factors such as volatilization rate, ventilation dilution, surface adsorption, etc.:

$C_{\text{gas}}(t)=\frac{V(t)\cdot V_{\text{chamber}}}{Q+k_{\text{loss}}\cdot V_{\text{chamber}}}\cdot (1-e^{-t/\tau })$

Integrated Model

By connecting the above four layers of mechanism chains in series, we obtain the complete forward prediction model:

$C_{\text{gas}}(t)=\mathrm{\Phi }\cdot N(t)\cdot W^{0.3156}\cdot \frac{S_{\text{total}}(t)}{K_s+S_{\text{total}}(t)}\cdot (1-e^{-t/\tau })$

Where:

• Φ is a constant that integrates all physical and chemical parameters
• $S_{\text{total}}(t)=0.144\cdot N(t)\cdot W^{0.3156}\cdot \rho \cdot C_{\text{sucrose}}\times {10}^{-9}$

This is the core application value of the model - Real-Time Conversion of the MeSA gas concentration detected by sensors into the field aphid population, providing data support for precise pesticide application decision-making.

Mathematical Derivation Process

Step 1: Define Simplification Parameters

$\mathrm{\Lambda }=0.144\cdot W^{0.3156}\cdot \rho \cdot C_{\text{sucrose}}\times {10}^{-9}$

Step 2: Substitute into the forward model and simplify

$C_{\text{gas}}(t)=\mathrm{\Phi }\cdot N(t)\cdot W^{0.3156}\cdot \frac{\mathrm{\Lambda }\cdot N(t)}{K_s+\mathrm{\Lambda }\cdot N(t)}\cdot (1-e^{-t/\tau })$

Step 3: Define the comprehensive coefficient

$\mathrm{\Psi }=\mathrm{\Phi }\cdot \mathrm{\Lambda }\cdot W^{0.3156}\cdot (1-e^{-t/\tau })$

Step 4: Construct a quadratic equation about N(t)

$\mathrm{\Psi }\cdot N(t{)}^2-C_{\text{gas}}(t)\cdot \mathrm{\Lambda }\cdot N(t)-C_{\text{gas}}(t)\cdot K_s=0$

Step 5: Solve for the number of aphids using the root formula

$N(t)=\frac{C_{\text{gas}}(t)\cdot \mathrm{\Lambda }+\sqrt{(C_{\text{gas}}(t)\cdot \mathrm{\Lambda }{)}^2+4\mathrm{\Psi }\cdot C_{\text{gas}}(t)\cdot K_s}}{2\mathrm{\Psi }}$

The physical meaning of this analytical expression is clear:

• The first term of the molecule represents the linear response component
• The radical term reflects the nonlinear correction of Michaelis-Menten kinetics
• The 2Ψ in the denominator normalizes the space-time scale effect

Figure 7 comprehensively presents the system validation results of the aphid number-MeSA gas concentration relationship model. The forward model (Figure 7A) exhibits typical Michaelis-Menten kinetic characteristics. The inverse model (Figure 7B) directly serves field applications—when the MeSA concentration > 2 ppb, the presence of aphids can be accurately identified. Logarithmic scale validation (Figure 7C) confirms that the model remains applicable within 4 orders of magnitude (10⁰ - 10³ individuals), conforming to the power-law scaling rule and covering the entire life cycle from initial invasion to the eve of an outbreak. Parameter sensitivity analysis (Figure 7D) reveals the key elements of system robustness, and quality control provides clear guidance: the culture conditions and inoculation concentration of engineered bacteria must be standardized, while the impact of environmental factors can be offset through regular calibration.

More importantly, the real-time pest population data provided by Model 1 will be directly input into the subsequent Model 2 and Model 3, forming a complete decision-making chain: when low-density infestation during the growth period is detected, the phloem injection plan will be initiated; when high-density spread during the outbreak period is identified, the Metarhizium anisopliae spraying strategy will be switched to, truly realizing the intelligent closed-loop management of "monitoring - diagnosis - pesticide application".

Toxoptera citricida obtain nutrients by sucking phloem sap through their piercing-sucking mouthparts during the growth period, causing serious damage to the citrus industry. Our innovative solution is to directly inject RNAi pesticides into the citrus phloem, using the aphids' feeding behavior to achieve targeted control. To guide farmers in determining the optimal injection time, injection location, and injection dosage while ensuring insecticidal efficacy and minimizing pesticide usage, we have developed a complete prediction system from "pesticide application" to "aphid death".

When building a prediction model, we face a core challenge: experiments cannot exhaustively enumerate the mortality rates under all "concentration-time" combinations. We need to predict the insecticidal effect under any condition based on limited wet experiment data. During the Human Practice communication, experts specifically reminded us that, considering the limitation of experimental data volume, we should adopt algorithms suitable for small samples, rather than relying on Machine Learning methods that depend on Big data.

Based on the cumulative characteristics of pesticide toxicity and the characteristics of experimental data, we adopted the modified Haber's Law ^[9] to establish a concentration × time (C×T) cumulative model:

$\text{Mortality Rate}=1-e^{(-k\times C^n\times t^m)}$

Where: C is the concentration (ng/μl), t is the time (h), k is the toxicity intensity parameter, n represents the degree of nonlinear influence of concentration on the toxic effect, and m represents the degree of nonlinear influence of time on the toxic effect.

The C×T model demonstrated excellent fitting performance across all four generations of our RNA pesticides, proving its applicability in describing the insecticidal kinetics of RNA pesticides. Learn more in Results.

Table 1. Fitted parameters for the C×T cumulative lethal model across four generations of RNAi pesticides

We selected the fourth-generation RNA molecule, bi-amiRNA, as the reference drug for subsequent simulations, and its predicted mortality equation is:

$\text{Mortality Rate}=1-e^{(-0.000782\times C^{0.3239}\times t^{1.1349})}$

This curve is the core foundation of the entire Model 2, enabling us to accurately predict the mortality probability of aphids based on the cumulative pesticide dose they ingest at any spatio-temporal point.

Aphids primarily feed on the phloem of fresh young leaves, which are part of the plant's "sink" organs and need to obtain nutrients from "source" organs such as mature leaves. According to the pressure flow theory (Münch Hypothesis) of plant phloem transport, phloem sap spontaneously flows from "source" to "sink". This physiological characteristic provides a natural advantage for our strategy: RNAi pesticides encapsulated by MS2 virus-like particles will automatically accumulate in aphid aggregation areas along with the pressure flow, achieving passive targeted delivery.

The concentration change of RNAi pesticides in the phloem is mainly controlled by three physical processes:

• Diffusion Process: Molecular motion driven by concentration layer
• Convection process: Pesticide transport driven by phloem sap flow
• Degradation Process: Natural Degradation of MS2 Capsid

Step 1: Select the coordinate system

Considering that the phloem is a tubular structure, we use the cylindrical coordinate system (r, θ, z) for description:

• r: radial coordinate (from tube wall to center)
• z: Axial coordinate (along the pipeline direction)
• Due to axial symmetry, the θ direction is ignored

Step 2: Establish the principle of mass conservation

According to the law of conservation of mass, the rate of change of concentration within the control volume = inflow - outflow - degradation - absorption

Step 3: Mathematical Representation of Various Physical Processes

(1) Diffusion term

Molecular diffusion follows Fick's first law, and the diffusion flux density is:

$J_x^{\text{diff}}=-D\frac{\partial C}{\partial x}$

In cylindrical coordinates, the complete diffusion term expands to:

$\nabla \cdot (D_{\text{eff}}\nabla C)=D_{\text{eff}}(\frac{{\partial }^{2}C}{\partial r^2}+\frac{1}{r}\frac{\partial C}{\partial r}+\frac{{\partial }^{2}C}{\partial z^2})$

(2) Convective term

Phloem sap flow mainly occurs axially, and the velocity distribution exhibits a parabolic shape (similar to Poiseuille flow):

$\nabla \cdot ({𝐯}_{\text{eff}}C)=v(r)\frac{\partial C}{\partial z}$

(3) Degradation term

The degradation of the MS2 virus-like particle capsid follows first-order reaction kinetics:

$-k_{\text{deg}}C,\text{ where }{k}_{\text{deg}}=\frac{\ln 2}{t_{1/2}}=\frac{\ln 2}{6.02\text{days}}$

Final Convection-Diffusion Equation:

$\frac{\partial C}{\partial t}+\nabla \cdot ({𝐯}_{\text{eff}}C)=\nabla \cdot (D_{\text{eff}}\nabla C)-k_{\text{deg}}C$

Boundary Conditions:

• Radial: ${\frac{\partial C}{\partial r}|}_{r=R}=0$
• Axial: ${C|}_{z=0}=C_0\delta (z-z_{\text{inj}})$

Aphids obtain sap from the phloem through their piercing-sucking mouthparts, a process that can be described by the Hagen-Poiseuille equation in fluid mechanics, and previous research ^[10] has confirmed that this equation is highly consistent with the phloem pressure flow theory.

Modified Hagen-Poiseuille Equation

The aphid stylet tapers gradually from the base to the tip (conical structure), and the standard circular tube formula cannot be directly applied:

• Base radius: R₀
• Tip radius: Rᵢ
• Mouthpiece needle length: L

Therefore, we have made corrections to the equation, assuming that the radius linearly tapers along the length direction of the needle. Let the coordinate along the length direction be x (0 ≤ x ≤ L), then we have:

$r(x)=R_0-\frac{R_0-R_i}{L}$

Assume that the pressure decreases linearly along the length direction of the edge needle:

$p(x)=p_0-\frac{\mathrm{\Delta }p}{L}x$

Substitute r(x) into the Hagen-Poiseuille equation:

$Q={\int }_0^L\frac{\pi }{8\mu }{[R_0-\frac{R_0-R_i}{L}x]}^4\frac{\mathrm{\Delta }p}{L}dx$

Let $u=R_0-\frac{R_0-R_i}{L}x$, then:

$du=-\frac{R_0-R_i}{L}dx⟹dx=-\frac{L}{R_0-R_i}du$

Integration boundary: when x = 0, u = R₀; when x = L, u = Rᵢ

Final flow formula:

$Q=\frac{\pi \mathrm{\Delta }P}{8\mu L}\cdot \frac{R_0^5-R_i^5}{5(R_0-R_i)}$

Calculation of Pesticide Intake

Instantaneous intake rate:

$S_{\text{sink}}(r,z,t)=Q\times C_{\text{local}}(r,z)$

Among them $C_{\text{local}}(r,z)$ is the local pesticide concentration at the location of the aphid (given by the previous spatio-temporal distribution model).

Daily cumulative intake:

$\text{Daily Intake}={\int }_{t_{\text{start}}}^{t_{\text{end}}}(\frac{\pi \mathrm{\Delta }P}{8\mu L}\cdot \frac{R_0^5-R_i^5}{5(R_0-R_i)})\times C(r,z,t)dt$

By combining the spatio-temporal distribution of aphid feeding, we can predict the actual intake dose of aphids in different regions, which provides a key input for the calculation of mortality in subsequent multi-objective optimization.

Based on the previously established model, we construct a multi-objective optimization framework to determine the optimal injection strategy. Our core idea is to minimize pesticide usage and environmental risk while ensuring the effectiveness of aphid control.

Decision Variable Definition

The injection strategy is fully described by the following decision variables:

$𝐗=\{𝐙,𝐑,𝐕,C_{\text{inj}},𝐓\}$

Where:

• $𝐙=[z_1,z_2,\dots ,z_n],\text{ where }{z}_i\in [0,L_{\text{trunk}}],\text{ represents the injection position coordinates along the phloem trunk}$
• $𝐑=[r_1,r_2,\dots ,r_n],\text{ where }{r}_i\in [0,R_{\text{phloem}}],\text{ represents the radial injection depth}$
• $𝐕=[V_1,V_2,\dots ,V_n],\text{ where }{V}_i\in [V_{\text{min}},V_{\text{max}}],\text{ represents the pesticide volume at each injection point}$
• $C_{\text{inj}}\in [C_{\text{min}},C_{\text{max}}],\text{ represents the injection liquid concentration}$
• $𝐓=[t_1,t_2,\dots ,t_n],\text{ where }{t}_i\in [0,24],\text{ represents the injection time (in hours)}$

Objective Function Construction

Objective 1: Maximize aphid mortality efficiency

$\max f_1={\int }_{{\mathrm{\Omega }}_{\text{aphid}}}{\int }_{T_{\text{effective}}}{P}_{\text{death}}(C_{\text{local}}(r,z,t),t)\cdot {\rho }_{\text{aphid}}(r,z)dtd\mathrm{\Omega }$

Where:

• ${\mathrm{\Omega }}_{\text{aphid}}$: Aphid distribution area (mainly tender leaf phloem)
• $T_{\text{effective}}$: Effective action time window
• $P_{\text{death}}(C,t)$: Death probability function based on lethal curve
• ${\rho }_{\text{aphid}}(r,z)$: Aphid spatial distribution density

Objective 2: Minimize the total cost of pesticides

$\min f_2={\sum }_{i=1}^n{V}_i\cdot C_{\text{inj}}+C_{\text{waste}}$

Among them, $C_{\text{waste}}$ is the waste cost, defined as the sum of the amount of pesticides lost to non-target areas and the amount of pesticides lost due to degradation:

$C_{\text{waste}}={\int }_{{\mathrm{\Omega }}_{\text{non-target}}}{\int }_{T}C(r,z,t)dtd\mathrm{\Omega }+{\int }_{\mathrm{\Omega }}{\int }_T{k}_{\text{deg}}C(r,z,t)dtd\mathrm{\Omega }$

Constraints

The optimization problem must satisfy the following constraints to ensure the safety, feasibility, and effectiveness of the solution:

1. Pesticide Residue Constraint: $C_{\text{fruit}}(t)\le C_{\text{MRL}},\forall t\in [0,T_{\text{harvest}}]$
2. Environmental safety constraints: ${\int }_{{\mathrm{\Omega }}_{\text{soil}}}C(r,z,t)d\mathrm{\Omega }\le C_{\text{env}}^{\text{max}},\forall t$
3. Feasibility constraints for injection location: $z_{\text{min}}\le z_i\le z_{\text{max}},r_i\le R_{\text{phloem}}(z_i)$
4. Time Window Constraint: $t_{\text{min}}\le t_i\le t_{\text{max}},|t_i-t_j|\ge \mathrm{\Delta }{t}_{\text{min}}$

Solution Strategy: NSGA-II Algorithm

This optimization problem has the following notable characteristics:

1. Multi-objective conflict: There is an inherent contradiction between f1 and f2 - improving insecticidal efficiency often requires increasing the amount of pesticide used
2. Strong Nonlinearity: The objective function involves PDE solving and a nonlinear dose-response relationship
3. High-dimensional decision space: If 8 injection points are considered, the dimension of the decision variables reaches 40 dimensions
4. Complex constraint coupling: Multiple constraint conditions are interrelated, and the feasible region is non-convex

These characteristics make traditional single-objective optimization methods (such as gradient descent, linear programming) difficult to apply. We need to find the Pareto optimal solution set, that is, the set of solutions where there is no other alternative that can improve one objective without compromising another.

Considering the complexity of the problem, we adopted NSGA-II (Non-dominated Sorting Genetic Algorithm II) to solve it.

The five Pareto optimal solutions identified by the algorithm validate the effectiveness of the model: the injection points of all solutions autonomously converge to the "source-sink" transition zone of the phloem and approach the aphid aggregation center, demonstrating that the optimization framework can identify key physiological nodes. These five solutions provide diverse decision-making options - from low-cost solutions to high-efficiency control solutions, among which Solution 3 achieves the optimal balance between efficiency and cost and can serve as a standard pesticide application strategy. Compared with the traditional experience-driven, single high-dose injection method, our framework realizes data-driven precise decision-making: dynamically adjusting the injection position according to the tree's physiological structure, optimizing pesticide configuration based on the spatial distribution of aphids, and revealing the efficiency-cost trade-off law through the Pareto front, providing a scientific decision-making basis for different demand scenarios.

It should be noted that parameters such as cost coefficients and concentration ranges in this study are all normalized relative values, mainly used to demonstrate the decision-making logic and methodological effectiveness of the optimization framework, rather than to provide specific economic accounting. This design of "universal framework, configurable parameters" makes the model highly scalable - in practical applications, simply substituting the measured physiological parameters, pesticide toxicology data, and local economic indicators into the framework can generate the optimal solution adapted to specific scenarios. This is precisely the core concept of the precision agriculture decision support system: providing scientific optimization tools to elevate trunk injection technology from "experience-based operation" to "Data drive intelligent decision-making".

When aphids enter the outbreak period, a large number of winged adults are produced and begin to spread. At this time, the phloem injection strategy has difficulty in dealing with aphid outbreaks across the entire orchard. We then adopted the RNAi-Metarhizium anisopliae spraying strategy for prevention and control. To guide farmers in determining the optimal spraying time, spraying location, and spraying dosage, and to achieve efficient, economical, and environmentally friendly aphid prevention and control, we established the following model.

Notably, Model 3 is highly consistent with Model 2 in its overall framework, also calculating aphid mortality based on the lethal curve, spatio-temporal distribution of pesticides, and pesticide intake, and finding the optimal pesticide application method through a multi-objective optimization model. The core ecological processes such as aphid growth and reproduction, natural enemy predation, winged morph transformation, and spatial diffusion are even more completely consistent with the background model in their mathematical descriptions in both models, reflecting our Modularization design concept.

Therefore, in Model 3, to avoid repetitive presentation of the same modules, we no longer conduct result simulations for each model module, but instead elaborate on the proposed model concepts and theoretical frameworks.

As an entomopathogenic fungus, the lethal mechanism of Metarhizium anisopliae is significantly different from that of traditional chemical pesticides. Metarhizium anisopliae spores first need to adhere to the aphid's body surface, then germinate and invade, and finally reproduce and cause death inside the insect body. When we integrate the RNAi mechanism into Metarhizium anisopliae, the lethal process becomes more complex: on the one hand, the hyphae of Metarhizium anisopliae damage the internal organs of aphids and consume nutrients; on the other hand, the released dsRNA interferes with the expression of key genes in aphids, forming a dual killing mechanism.

Cumulative Lethal Probability (considering spore concentration, time, and environmental factors):

$P_{\text{death}}(C_{\text{spore}},t,E)=1-\exp (-{\int }_0^t{k}_{\text{infection}}(C_{\text{spore}},E(\tau ))d\tau )$

Infection Rate (Michaelis-Menten Kinetics + Environmental Modification):

$k_{\text{infection}}(C_{\text{spore}},E)=k_0\cdot \frac{C_{\text{spore}}}{K_m+C_{\text{spore}}}\cdot f_{\text{temp}}(T)\cdot f_{\text{humidity}}(RH)$

Environmental Impact Function:

$f_{\text{temp}}(T)=\exp (-\frac{{(T-T_{\text{opt}})}^2}{2{\sigma }_T^2})$

$f_{\text{humidity}}(RH)=\frac{RH-R{H}_{\min }}{R{H}_{\text{opt}}-R{H}_{\min }},RH\ge R{H}_{\min }$

The diffusion of Metarhizium anisopliae spores as a biopesticide in the field is a complex multi-physical coupling process, involving the synergistic action of multiple mechanisms such as molecular diffusion, wind field convection, environmental degradation, and foliar adsorption. Different from chemical pesticides, Metarhizium anisopliae spores are biologically active, and their behavior in the environment is significantly affected by various environmental factors such as temperature, humidity, and UV radiation.

Based on continuum mechanics and transport phenomena theory, we develop a partial differential equation model to describe the spatiotemporal evolution of spores.

Dynamic Wind Field Simulation

In the small-scale field environment, the wind field exhibits obvious time-varying characteristics and randomness. Different from the atmospheric boundary layer flow, the wind field inside the orchard is affected by factors such as the plant canopy and microtopography, showing complex turbulent characteristics.

Based on Data Analysis simulation, we establish a dynamic wind field model that includes multiple time scales:

$v_{\text{wind}}(t)=\max (0.3,v_{\text{base}}+\mathrm{\Delta }v(t))$

Among them, the wind speed disturbance includes deterministic periodic variations and random fluctuations:

$\mathrm{\Delta }v(t)=0.4\sin \left(\frac{2\pi t}{2}\right)+0.15\xi (t)$

The wind direction in the field is affected by multiple factors, presenting a superposition of periodic drift and short-term disturbances:

${\theta }_{\text{wind}}(t)={\theta }_{\text{base}}+{\theta }_{\text{drift}}(t)+{\theta }_{\text{fluctuation}}(t)$

• ${\theta }_{\text{drift}}(t)=0.4\sin \left(\frac{2\pi t}{4}\right)$
• Short-term Disturbance: ${\theta }_{\text{fluctuation}}(t)=0.3\sin \left(\frac{2\pi t}{0.3}\right)$

Wind speed component calculation:

$v_x(t)=v_{\text{wind}}(t)\cos ({\theta }_{\text{wind}}(t))$

$v_y(t)=v_{\text{wind}}(t)\sin ({\theta }_{\text{wind}}(t))$

Control equations based on transport phenomena

Considering the evolution of spore concentration C(x,y,t) within a two-dimensional spatial domain, based on the law of conservation of mass:

$\frac{\partial C}{\partial t}+\nabla \cdot J=-R$

where $J$ is the spore flux density vector, and R is the net consumption rate

Flux density composition:

$J=J_{\text{diff}}+J_{\text{conv}}=-D_{\text{eff}}\nabla C+Cv$

Substituting into the mass conservation equation, we obtain:

$\frac{\partial C}{\partial t}+v\cdot \nabla C=\nabla \cdot (D_{\text{eff}}\nabla C)-R$

Multi-physics Process Modeling

Environmentally Modified Diffusion Process

The molecular diffusion of Metarhizium anisopliae spores is significantly affected by environmental factors, and an environmentally corrected effective diffusion coefficient is established:

$D_{\text{eff}}(t)=D_{\text{base}}\cdot f_{\text{humidity}}(t)\cdot f_{\text{surfactant}}$

Humidity influence mechanism: Spores maintain better viability and mobility in high humidity environments. A 12-hour humidity cycle model is established:

$f_{\text{humidity}}(t)=0.8+0.3\sin \left(\frac{2\pi t}{12}\right)$

Surfactant Effect: The surfactant in the formulation enhances the dispersal ability of spores:

$f_{\text{surfactant}}=1.2$

Wind-driven convective transport

Based on the theory of fluid mechanics, the convection term driven by the wind field is:

$\text{Convection term}=v_x(t)\frac{\partial C}{\partial x}+v_y(t)\frac{\partial C}{\partial y}$

This item describes the macroscopic transport of spores with the wind field, which is the main mechanism for long-distance spread.

Environmental attenuation process

UV Photolysis Attenuation: UV radiation is the primary environmental factor leading to spore inactivation. A time-varying attenuation model is established:

$R_{UV}=k_{\text{decay}}^{\text{eff}}(t)\cdot C$

The time-varying attenuation coefficient takes into account daily periodic variations and UV protection effects:

$k_{\text{decay}}^{\text{eff}}(t)=k_{\text{decay}}^{\text{base}}\cdot [1+f_{UV}(t)\cdot (1-f_{\text{protection}})]$

• UV daily cycle factor: $f_{UV}(t)=0.5+0.5\cos \left(\frac{2\pi t}{24}\right)$
• UV Protection Factor: $f_{\text{protection}}=0.7$

Leaf surface attachment process

Based on the theory of surface physical chemistry, the attachment of spores to plant leaf surfaces follows the Langmuir isothermal adsorption model:

$R_{\text{adhesion}}=k_{\text{attach}}\cdot C\cdot \frac{C}{C+C_{\text{sat}}}$

This non-linear term reflects the saturation effect of the attachment check point, where $C_{\text{sat}}=\frac{C_0}{10}$ is the saturation concentration parameter.

Complete system of governing equations

By integrating all the above physical processes, a complete mathematical model for the spatio-temporal distribution of Metarhizium anisopliae spores is obtained:

$\frac{\partial C}{\partial t}+v_x(t)\frac{\partial C}{\partial x}+v_y(t)\frac{\partial C}{\partial y}=D_{\text{eff}}(t)(\frac{{\partial }^{2}C}{\partial x^2}+\frac{{\partial }^{2}C}{\partial y^2})-k_{\text{decay}}^{\text{eff}}(t)C-k_{\text{attach}}C\frac{C}{C+C_{\text{sat}}}$

Figure 14 shows the entire spatio-temporal evolution process of Metarhizium anisopliae spores under the triangular distribution spraying strategy, revealing the rapid response capability of biological control. It only takes about 22 hours for the entire 40m×40m orchard to achieve a high coverage rate of 90.0%. This rapid diffusion process benefits from the synergistic effect of wind field-driven convective transport and molecular diffusion. This means that during the aphid outbreak period, the vast majority of winged aphids will come into contact with an effective concentration of Metarhizium anisopliae spores within 24 hours, and the spores can quickly and automatically adhere to the body surface of flying and diffusing aphids, then germinate, invade, and release dsRNA for dual killing, effectively blocking the orchard-wide spread and population explosion of aphids, demonstrating the rapid coverage, active tracking, and efficient killing unique advantages of the RNAi-Metarhizium anisopliae strategy in dealing with large-scale winged aphids during the outbreak period.

Considering that the winged individuals of the Toxoptera citricida do not possess long-distance migratory capabilities, and their spread is mainly achieved through short-distance flight (usually between adjacent plants), we established an improved transmission dynamics equation based on the classical SEIR model, combined with the local spread characteristics of aphids.

Population state classification

Divide the aphid population at each grid point (i,j) in the orchard into three core states, with the total population size satisfying $N_{i,j}(t)=S_{i,j}(t)+E_{i,j}(t)+I_{i,j}(t)+R_{i,j}(t)$

• Susceptible Aphids (Susceptible) - healthy, uninfected individuals;
• Exposed Aphids (Infected) - individuals that have been exposed to Metarhizium anisopliae but show no signs of infection.
• Infected Aphids (Infected) - individuals that have been infected by Metarhizium anisopliae but have not died;
• Removed Aphids (Removed) - individuals that died from infection or recovered.

Model Improvement Elements

1. Aphid Spread

Number of winged aphids:

$W_{i,j}(t)=\phi \left(\frac{N_{i,j}}{K_{i,j}}\right)\cdot N_{i,j}$

Among them

$\phi (x)=0.05+\frac{0.45}{1+e^{-10(x-0.7)}}$

Migration Traffic:

$J_{(i,j)\to (k,l)}=m\cdot W_{i,j}(t)\cdot \sigma \left(\frac{N_{i,j}}{K_{i,j}}\right)\cdot (1-\frac{N_{k,l}}{K_{k,l}})$

Among them $(k,l)\in N_{i,j}=\{(i-1,j),(i+1,j),(i,j-1),(i,j+1)\}$

Net Diffusion Amount:

$D_{i,j}(t)={\sum }_{(k,l)\in N_{i,j}}^{}[J_{(k,l)\to (i,j)}-J_{(i,j)\to (k,l)}]$

2. Metarhizium dissemination dynamics

The propagation coefficient takes into account the local spore concentration:

${\beta }_{i,j}(t)={\beta }_0\cdot C_{\text{spore}}(i,j,t)\cdot f_{\text{temp}}(T)\cdot f_{\text{humidity}}(RH)$

Final Spatial SEIR Dynamic Equations

For each grid point $(i,j)$:

$\frac{d{S}_{i,j}}{dt}=r_{\text{birth}}\cdot S_{i,j}\cdot (1-\frac{N_{i,j}}{K_{i,j}})-{\beta }_{i,j}(t)\frac{S_{i,j}{I}_{i,j}}{N_{i,j}}+D_{i,j}^S(t)-\mu S_{i,j}$

$\frac{d{E}_{i,j}}{dt}={\beta }_{i,j}(t)\frac{S_{i,j}{I}_{i,j}}{N_{i,j}}-\sigma E_{i,j}+D_{i,j}^E(t)-\mu E_{i,j}$

$\frac{d{I}_{i,j}}{dt}=\sigma E_{i,j}-\gamma I_{i,j}-\alpha I_{i,j}+D_{i,j}^I(t)-\mu I_{i,j}$

$\frac{d{R}_{i,j}}{dt}=\gamma I_{i,j}+\alpha I_{i,j}+D_{i,j}^R(t)-\mu R_{i,j}$

Among them, the diffusion term is allocated according to the proportion of aphids in each state:

$D_{i,j}^X(t)=D_{i,j}(t)\cdot \frac{X_{i,j}}{N_{i,j}},X\in \{S,E,I,R\}$

By defining decision variables, objective functions, and constraints, and combining them with improved algorithms, the optimal spraying strategy is solved.

Decision Variable Definition

1. Spatial Decision Variable

• $X_{\text{spray}}=\{x_1,x_2,...,x_m\}$: x-coordinate of the spraying point;
• $Y_{\text{spray}}=\{y_1,y_2,...,y_m\}$: y-coordinate of the spraying point;
• $H_{\text{spray}}=\{h_1,h_2,...,h_m\}$: Spraying height.

2. Dose decision variable

• $Q_{\text{spray}}=\{Q_1,Q_2,...,Q_m\}$: Spore release amount at each spraying point;
• $C_{\text{spore}}^{\text{spray}}$: Spore concentration in the spray liquid.

3. Time decision variable

• $T_{\text{spray}}=\{t_1,t_2,...,t_m\}$: Spraying time at each spraying point;
• $\mathrm{\Delta }{t}_{\text{interval}}$: Time interval between multiple sprays.

Objective Function Construction

Objective Function 1: Maximize Population Control Efficiency

$\max f_1=\frac{{\sum }_{i,j}^{}(N_{i,j}(0)-N_{i,j}\left(T_{\text{final}}\right))}{{\sum }_{i,j}^{}{N}_{i,j}(0)}\times 100\%$

Where:

• $T_{\text{final}}$: Evaluation time endpoint (usually 15 - 20 days);
• $w(t)$: Weight function (reflecting the importance of control effects at different times);
• $R(t)$: Number of aphids removed;
• $N(0)$: Initial aphid population size.

Simplified expression (population reduction rate):

$\max f_1=\frac{N(0)-N\left(T_{\text{final}}\right)}{N(0)}\times 100\%$

Objective Function 2: Minimization of Metarhizium anisopliae Usage

$\min f_2={\sum }_{i=1}^m{Q}_i+C_{\text{waste}}$

Among them, waste cost $C_{\text{waste}}$ includes three parts:

• Settlement Loss: $L_{\text{settling}}={\int }_{{\mathrm{\Omega }}_{\text{ground}}}^{}{\int }_0^T{k}_{\text{settle}}{C}_{\text{spore}}(x,y,t)dtdA$
• Wind Drift Loss: $L_{\text{drift}}={\int }_{{\mathrm{\Omega }}_{\text{non-target}}}^{}{\int }_0^T|v_{\text{wind}}|C_{\text{spore}}(x,y,t)dtdA$
• $L_{\text{decay}}={\int }_{\mathrm{\Omega }}^{}{\int }_0^T{k}_{\text{decay}}{C}_{\text{spore}}(x,y,t)dtdV$

Objective Function 3: Maximize Environmental Safety

$\max f_3=1-\frac{1}{|{\mathrm{\Omega }}_{\text{sensitive}}|}{\int }_{{\mathrm{\Omega }}_{\text{sensitive}}}^{}{\max }_t{C}_{\text{spore}}(x,y,t)dA$

Among them, ${\mathrm{\Omega }}_{\text{sensitive}}$ is the environmentally sensitive area (such as water bodies, natural enemy habitats), $|{\mathrm{\Omega }}_{\text{sensitive}}|$ is the area of the sensitive area.

Constraints

1. Ecosystem Safety Constraints: ${\int }_{{\mathrm{\Omega }}_{\text{beneficial}}}^{}{C}_{\text{spore}}(x,y,t)dA\le C_{\text{threshold}},\forall t$

   ${\mathrm{\Omega }}_{beneficial}$as habitats for beneficial organisms,$C_{threshold}$as the safe threshold for spores

2. Spraying Equipment Constraints: $Q_{\min }\le Q_i\le Q_{\max },h_{\min }\le h_i\le h_{\max }$
   $Q_{\min }/Q_{\max }$ represents the upper and lower limits of spore release amount, $h_{\min }/h_{\max }$ represents the upper and lower limits of spraying height
3. Spraying Equipment Constraints: $t_{\text{dawn}}\le t_i\le t_{\text{dusk}},RH\left(t_i\right)\ge R{H}_{\min }$
   Spraying time should be between day and night, and humidity should not be lower than the minimum effective humidity $R{H}_{\min }$
4. Coverage Constraint: ${\bigcup }_{i=1}^m{A}_{\text{coverage}}(x_i,y_i,h_i)\supseteq {\mathrm{\Omega }}_{\text{target}}$
   The coverage area of all spraying points $A_{\text{coverage}}$ must fully cover the target aphid area ${\mathrm{\Omega }}_{\text{target}}$

Final Spatial SEIR Dynamic Equations

Optimize the classic NSGA-II algorithm to address the randomness of Metarhizium anisopliae spread:

There are differences in the results of multiple simulations of the same strategy, so the mean value is calculated by taking multiple simulations:

${\bar{f}}_k=\frac{1}{N_{MC}}{\sum }_{j=1}^{N_{MC}}{f}_k^{(j)}$

($N_{MC}$ is the number of simulations, $f_k^{(j)}$ is the k-th objective function value of the j-th simulation)

Based on the monitoring of pest population dynamics, select different pest control methods —— the modeling system we have developed has established a full-chain mathematical modeling paradigm from mechanism understanding to precise intervention, providing a systematic solution for the prevention and control of complex agricultural pests. Modularization enables each sub-model to be independently verified and replaced, ensuring the flexibility of the system; parameters can be continuously calibrated through wet experimental data, guaranteeing the accuracy of the model; optimization algorithms such as NSGA-II are applicable to other constrained optimization problems, demonstrating the transferability of the method. This design concept makes our framework not only applicable to Toxoptera citricida control, but also extensible to other agricultural pest management scenarios.

The modeling framework still has vast room for expansion. We can introduce climate change scenario analysis to achieve IPCC model coupling, integrate remote sensing data to build a regional-scale early warning system, develop an interactive decision support system (Web/App) to serve practical applications, and explore reinforcement learning algorithms to implement self-adaptive optimization strategies.

As George Box said, "All models are wrong, but some are useful"—the value of a model lies not in being "perfect" but in being "useful". We expect this paradigm to provide methodological reference for interdisciplinary modeling research, promote mathematical modeling from describing phenomena to guiding practice, and ultimately realize the vision of smart agriculture.

All the model application code has been uploaded to Gitlab.