EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 4, Article Number 6691 ISSN 1307-5543 – ejpam.com Published by New York Business Global Impact of Awareness Programs on Whitefly Dynamics in Coconut Farming: A Unified DDE and SDE Modeling Approach B. Dhivyadharshini1, R. Senthamarai1,∗ 1 Department of Mathematics, College of Engineering and Technology, SRM Institute of Science and Technology, Kattankulathur-603203, Tamilnadu, India Abstract. In this article, we proposed a comprehensive pest control model that integrated both delay differential equations (DDEs) and stochastic processes to mitigate the spread of whiteflies in coconut plantations. The delay model introduced a time lag in the implementation of awareness programs and investigated its impact on the system’s equilibrium stability. Numerical simulations validated the model’s effectiveness in enhancing pest management. A stochastic component, formulated using a Wiener process within a system of non-linear ordinary differential equations (ODEs), was used to estimate the probability of disease elimination under environmental fluctuations. The study was novel in combining both deterministic delays and stochastic effects in a unified framework, offering deeper insights into timing and control efficiency. Although focused on coconut farming, the modeling approach and techniques had broader applicability in agricultural pest control scenarios. The findings enhanced our understanding of pest dynamics under uncertainty and delay, providing a foundation for more informed and effective control interventions in agricultural ecosystems. 2020 Mathematics Subject Classifications: 34A34, 92D45, 34K07, 34K20, 37H10, 65C30 Key Words and Phrases: Pest control model, farming awareness, delay differential equation, stochastic differential equation, numerical simulation 1. Introduction Cocos nucifera, the coconut tree, is a vital crop with diverse benefits, from its nutrient-rich fruit to its versatile coir [1]. Coconut diseases impacted livelihoods, despite its antiviral properties; key components included endosperm (kernel), endocarp (shell), and mesocarp (coir) [2]. India, the third-largest producer, relied heavily on coconut-based industries, with products ranging from food and cosmetics to artistic crafts and eco-friendly materials [3]. This resilient palm provided antiviral properties, essential minerals, and economic ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.6691 Email addresses: db0558@srmist.edu.in (B. Dhivyadharshini), senthamr@srmist.edu.in (R. Senthamarai) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 2 of 28 sustainability [4]−[6]. Celebrated on World Coconut Day (September 2) [7], intercropping with floriculture enhanced farming efficiency [8], [9]. With applications in ropes, mats, insulation, and even rehydration therapy, the coconut truly embodied a functional, life-enriching resource [10]−[14]. The rugose spiraling whitefly (RSW), Aleurodicus rugioperculatus, first appeared in scientific literature in 2004 when Martin identified it and initially named it the gumbo limbo spiraling whitefly [15]. Highly polyphagous, it fed on over 118 hosts across 43 plant families. It was first detected in Tamil Nadu and Kerala in 2016 [16], RSW damaged coconut plantations by extracting nutrients, excreting honeydew, and promoting sooty mold growth, thereby reducing photosynthesis [17]. Infestation signs included egg spirals, white wax deposits, sticky honeydew, and black sooty mold [18]−[20]. Escalating infestations threatened coconut trees, emphasizing the need for farmer awareness to sustain production. Mathematical models helped analyze interacting populations, particularly in epidemiology, by revealing interactions and parameter impacts [21]. These models aided in understanding epidemic dynamics and control strategies [22], [23]. Studies on RSW infestations in coconut trees suggested that reducing the contact rate effectively controlled infection [24]. Awareness programs played a crucial role in promoting knowledge and understanding among farmers, leading to the adoption of improved agricultural practices, increasing productivity, and providing significant economic and social benefits to farming communities [25]. Research has explored media-driven vaccination awareness, farming awareness in pest control [26], [27]. The reproduction number, analyzed via the next-generation matrix, was found to be crucial in epidemiology [28]. A mathematical model represents real-world systems using variables, equations, and parameters to describe relationships and predict behavior under different conditions. Equilibrium points help analyze long-term dynamics [29]−[31]. DDEs incorporate time lags, making them suitable for systems where current changes depend on past states [32], [33]. Stochastic differential equations (SDEs) introduce randomness, accounting for uncertainty in system behavior [34]−[37]. Our study integrated DDEs and SDEs to model RSW infestations and awareness programs. Inspired by [38], [39], we extended the model in [40] with DDEs and SDEs to control RSW infestations and safeguard coconut populations through awareness programs. Our research is focused specifically on the Pollachi locality of Tamil Nadu. Thus, the parameter values were selected based on the same zone. Section 2 begins with the development of the mathematical model, incorporating DDEs. It explores key properties such as positive invariance, boundedness, the basic reproduction number, and equilibrium analysis. The equilibrium analysis includes subsections on the tree-pest free equilibrium, pest-free equilibrium, and coexistence equilibrium, along with stability assessments under delay conditions. Section 3 introduces the SDE model, while Section 4 investigates parameter sensitivity. Section 5 presents results and discussion to validate the theoretical findings. Finally, Section 6 summarizes the key results and conclusions of the study. B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 3 of 28 2. Mathematical formulation of the problem In [40], G. Suganya and R. Senthamarai studied the impact of awareness on the dynamics of pest control in coconut trees. Here, we have extended the ODE model to both DDE and SDE model. The tree population is classified into healthy trees H and infected trees E to evaluate the whitefly’s effect on coconut trees. L refers to the population of whiteflies, and A indicates the number of awareness programs carried out over the time period t. 2.1. DDE model We have added a time delay in the awareness programs within the system of equations, which impacts the infected tree density and whitefly density in efforts to control infectious diseases. The delayed model is given by: The population of healthy coconut trees H(t) grows logistically, reduced by infection from whiteflies with a Holling type II functional response. The logistic term accounts for total tree density (healthy plus infected) competing for site resources; specifically, parameter s denotes the plantation’s tree density (see Table 1) and represents the total available resources (space, nutrients and light). Since infected trees, although reduced in health and productivity, remain physically present and continue to occupy space and consume resources, their presence contributes to the density-dependent limitation on healthy tree growth. dH(t) dt = rH(t) ( 1− H(t) + E(t) s ) − ϕH(t)L(t) 1 + γL(t) . (1) The population of infected trees E(t) increases by infections of healthy trees by whiteflies, modeled using a Holling type II functional response, and decreases due to natural mortality at rate k and awareness-driven control measures (with delay τ) at rate bA(t−τ)E(t). The awareness term bA(t− τ)E(t) represents the removal or recovery of infected trees through interventions such as roguing and replacement, which specifically target infected trees and not the vector population. In contrast, the whitefly population in Eq.(3) is affected by awareness through a separate term vA(t − τ)L(t), representing vector-specific control measures (e.g., pesticide application, biological control). Therefore, bA(t−τ)E(t) appears only in Eq.(2), while Eq.(3) includes its own awareness-related removal mechanism. Hence, dE(t) dt = ϕH(t)L(t) 1 + γL(t) − kE(t)− bA(t− τ)E(t). (2) The whitefly population L(t) is generated from infected trees at a per capita rate α, and declines through natural mortality at rate µ and awareness-mediated control actions that increase whitefly mortality after a delay. The term vA(t− τ)L(t) represents the reduction of the whitefly population due to awareness-driven interventions. Increased awareness (from A(t − τ)) can lead to timely implementation of control strategies such as targeted pesticide application, release of natural predators, or improved sanitation of infested areas. B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 4 of 28 These measures directly reduce whitefly numbers rather than only indirectly affecting them through host removal; hence, the term is modeled as a direct mortality rate proportional to both A(t− τ) and the current whitefly density L(t). Therefore, dL(t) dt = αE(t)− µL(t)− vA(t− τ)L(t). (3) The awareness variable A(t) increases due to a baseline implementation of awareness programs and local reporting of infected trees, and decays over time in the absence of reinforcement. This dynamics is described by dA(t) dt = δ + βE(t)− ξA(t). (4) With the initial values H(θ) = l, E(θ) = m, L(θ) = n, A(θ) = q, θ ∈ [−τ, 0] . (5) We assumed that l > 0, m > 0, n > 0, and q > 0 as these represent the initial positive population densities of healthy trees, infected trees, whiteflies, and the awareness level, respectively. This assumption ensured both the biological realism and the mathematical feasibility of the model. The same constant values l,m, n, q were assigned for θ ∈ [−τ, 0] to represent a steady, uniform initial history before t = 0, implying no changes in the state variables during the delay period. This common assumption in delay differential equation models isolates the effect of the delay from additional variability in the initial conditions. Figure 1: Schematic diagram B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 5 of 28 2.2. Positive invariance The Eqs.(1) − (4) can be represented compactly as follows: dU dt = Ψ(U). (6) Here, let U(θ) = (H(θ), E(θ), L(θ), A(θ))T ∈ C, where C = C ( [−τ, 0],R4 + ) represents the Banach space of continuous functions. The function Ψ = (Ψ1,Ψ2,Ψ3,Ψ4) T denotes the the right-hand terms in the system of Eqs.(1)−(4). Based on the fundamental theory of functional differential equations [41], a unique solution (H(t), E(t), L(t), A(t)) exists for the system described by Eqs.(1)−(4) given the initial conditions specified in Eq.(5). Theorem 1. All the solutions of Eqs.(1)−(4) with initial conditions Eq.(5) are positive. Proof. This theorem was validated by the method described in [42] and [43]. Indeed, it is straightforward to verify in Eqs.(1) − (4) that anytime determining H(θ) ∈ R+ with H = 0, E = 0, L = 0, A = 0. Then Ψi(U) |Ui=0,U∈R4 + ≥ 0. Using Lemma 2 from [43] and Theorem 1.1 from [42], any solution x(t) = x(t, x(θ)) of Eqs.(1) − (4) with x(θ) ∈ C is such that U(t) ∈ R4 + for all t ≥ 0. Since all solutions to Eqs. (1) − (4) are non - negative for any t > 0, we can deduce that the solution lies in the region R4 +. Consequently, the positive cone R4 + serves as an invariant region. 2.3. Boundedness Boundedness in a system reflects its proper behavior, as it guarantees that none of the interacting sub populations can rise to infinity in a finite amount of time. If we define N = H + E to be the total plant biomass at any time t, then combining the Eqs.(1) and (2) results in: dN dt = rE [ 1− H + E s ] − kE − bA(t− τ)E ≤ rN ( 1− N s ) , thus, it follows lim t→∞ supN ≤M := max{N(0), s}. Similarly, taking into account the aware population, we get dA dt = δ + βE − ξA ≤ δ + βE − ηA, giving again an upper bound lim t→∞ supA ≤ δ + βs η . B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 6 of 28 Again, dL dt = αE − µL− vA(t− τ)L ≤ αs− µL, and it follows lim t→∞ supL ≤ αs µ . Thus, the set determines the region of attraction D = { (H,E,L,A) ∈ R4 + : 0 ≤ H + E ≤M, 0 ≤ L ≤ αs µ , 0 ≤ A ≤ δ + βs η } , which attracts all solutions from the positive cone’s interior and is positively invariant. 2.4. Reproduction number At pest-free equilibrium, when τ = 0, the disease class Eqs.(2) and (3) characterize the population inputs in terms of the population size and the initial number of infections. The reproduction number is derived using the next-generation matrix approach introduced by Diekmann et al. [44]. R0 = ϕsα ξ2 (ξk + bδ)(µξ + vδ) . The product ϕsα in the numerator captures the combined effects of transmission from whiteflies to healthy trees (ϕs) and the production of new whiteflies by infected trees (α). The factor ξ2 arises from substituting the pest-free awareness level A∗ = δ/ξ into the next-generation matrix formulation, reflecting the role of awareness decay in both host and vector components. The denominator (ξk + bδ)(µξ + vδ) represents the effective removal rates of infected trees and whiteflies, combining natural mortality with awareness-driven control actions. 2.5. Equilibria and stability assessment of delayed system The following matrix form represents the linearization of the system defined by Eqs.(1)−(4): d dt  H(t) E(t) L(t) A(t)  = F  H(t) E(t) L(t) A(t) +G  H(t− τ) E(t− τ) L(t− τ) A(t− τ)  , where F and G are, F =  r − r(2H+E) s − ϕL 1+γL − rH s − ϕH (1+γL)2 0 ϕH 1+γL −k ϕH (1+γL)2 0 0 α −µ 0 0 β 0 −ξ  , G =  0 0 0 0 0 0 0 −bE 0 0 0 −vL 0 0 0 0  . (7) B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 7 of 28 Table 1: Model parameters and their corresponding values. Symbol Meaning Unit Value used for analysis [40] Range [40] s Tree density acre−1 70 60 − 70 r Replanting rate day−1 0.0005 0 − 0.003 k Mortality rate of tree day−1 0.0002 0 − 0.002 ϕ Contact rate pest−1day−1 0.0002 0 − 0.002 α Whitefly birth rate day−1 0.2 0.1 − 0.3 µ Death rate of whitefly day−1 0.06 0.06 − 0.1 γ Saturation constant – 0.2 0.2 b Maximum activity rate day−1 0.0001 0 − 0.0001 v Death rate of whitefly due to awareness day−1 0.005 0.005 δ Rate of awareness programs day−1 0.03 0.03 β Rate of local awareness day−1 0.025 0.025 ξ Fading rate of awareness day−1 0.015 0.015 The characteristic equation of the system can be written as Υ(σ) =| σI − F − e−στG |= 0. (8) The system’s equilibrium points can be represented as: (i) Tree−pest free equilibrium: E0 = (0, 0, 0, δξ ), (ii) Pest−free equilibrium: E1 = (s, 0, 0, δξ ), (iii) Coexistence equilibrium: Ē = (H∗, E∗, L∗, A∗), where, H∗ = (k + bA∗) β(µ+ vA∗) + γα(ξA∗ − δ) ϕβα , E∗ = ξA∗ − δ β , L∗ = α(ξA∗ − δ) β(µ+ vA∗) . The coexistence equilibrium is determined by the positive root A∗ of the cubic equation p3A 3 + p2A 2 + p1A+ p0 = 0, (9) where the coefficients pi are algebraic combinations of model parameters: p3 = rbµγvβ + rbγ2ξ2, p2 = rkβγξv + rkγ2αξ2 + rbβ2vµ+ rbµξγβ − rbµγδβv + rγϕαξ2 + rbγξαβ, p1 = rkµγξβ + rkβ2v + rkγαξβ + rbµβ2 + rbγ2αδ2 + ξrϕβα− rsϕαξγ − rkβ2vγδ − 2rkγ2αξδ − rbµδγβ − rbγαδβ − ξrϕαγδ − rδαϕγξ, p0 = rsϕβαγδ + rkµβ2 − rkµγδβ − rsϕβ2α− rkγαδβ + rkγ2αδ2 − rkϕβα+ rδ2γϕα+ sδϕ2βα. B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 8 of 28 The biologically admissible root must satisfy A∗ > δ ξ , so that E∗ > 0. Among the real roots, we selected the unique value fulfilling this condition as the biologically meaningful equilibrium for subsequent stability and simulation analyses. 2.5.1. Stability analysis at the tree–pest free equilibrium: E0 = ( 0, 0, 0, δξ ) Lemma 1. The tree–pest free equilibrium E0 is always unstable [24]. Proof. At E0, the Jacobian matrix of the system has at least one eigenvalue with a positive real part for all biologically admissible parameter values. This implies that any small perturbation from E0, such as the introduction of a small number of whiteflies or infected trees, will grow over time. Biologically, in the absence of healthy and infected trees as well as pests, the introduction of pests leads to growth due to the lack of competition and absence of control measures. Hence, E0 is unstable for all parameter values. 2.5.2. Stability analysis at the pest−free equilibrium: E1 = ( s, 0, 0, δξ ) By substituting the equilibrium point into the matrices presented in Eq.(7), we can derive F and G: F =  −r −r −ϕs 0 ϕs −k ϕs 0 0 α −µ 0 0 β 0 −ξ  , G =  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  . (10) Substituting F and G in Eq.(8), the transcendental equation derived from the matrix is ψ(σ, τ) = σ4 +A1σ 3 +A2σ 2 +A3σ +A4 = 0, (11) with coefficients: A1 = k + µ+ r + ξ, A2 = kµ− αϕs+ kr+ µr+ ϕrs+ kξ + ξµ+ rξ, A3 = αϕ2s2 − αϕrs+ kµr + µϕrs− αϕsξ + kµξ + krξ + µrξ + ϕrsξ, A4 = αϕ2−αϕrsξ+kµrξ+µϕrsξ. Theorem 2. The pest free equilibrium E1 is locally asymptotically stable (LAS) for R0 < 1 and unstable otherwise. Proof. At E1, the system is locally asymptotically stable if it satisfies the Routh–Hurwitz conditions: A4 > 0, A1A2 −A3 > 0, and (A1A2 −A3)A3 −A2 1A4 > 0. These inequalities hold precisely when R0 < 1, ensuring E1 is locally asymptotically stable; otherwise, it is unstable. B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 9 of 28 2.5.3. Stability analysis at the Coexistence equilibrium: Ē = (H∗, E∗, L∗, A∗) When we substitute the equilibrium point into the matrices from Eq.(7), we obtain F and G: F =  r − r(2H∗+E∗) s − ϕL∗ 1+γL∗ − rH∗ s − ϕH∗ (1+γL∗)2 0 ϕH∗ 1+γL∗ −k ϕH∗ (1+γL∗)2 0 0 α −µ 0 0 β 0 −ξ  , G =  0 0 0 0 0 0 0 −bE∗ 0 0 0 −vL∗ 0 0 0 0  . (12) Substituting F and G in Eq.(8), the transcendental equation derived from the matrix is, ψ(σ, τ) = σ4 +A1σ 3 +A2σ 2 +A3σ +A4 + e−στ [B1σ 2 +B2σ +B3] = 0, (13) with coefficients: A1 = ξ + µ+ k − r + r(2H∗ + E∗) s , A2 = µξ+kξ+kµ− ϕH∗α (1 + γL∗)2 − rξ− rµ− rk+ r(2H∗ + E∗)ξ s + r(2H∗ + E∗)µ s + r(2H∗ + E∗)k s + rH2∗ϕ s(1 + γL∗) + ϕ2H2∗ (1 + γL∗)3 , A3 = µkξ− ϕH∗αξ (1 + γL∗)2 − ϕH∗βvL∗ (1 + γL∗)2 −rµξ−rkξ−rkµ+ rϕH∗α (1 + γL∗)2 + r(2H∗ + E∗)µξ s + r(2H∗ + E∗)kξ s + r(2H∗ + E∗)kµ s − r(2H∗ + E∗)(ϕH∗α) s(1 + γL∗)2 + rH2∗ϕξ s(1 + γL∗) + rH2∗ϕµ s(1 + γL∗) + ϕH∗ξ 1 + γW ∗ , A4 = ϕrH∗αξ (1 + γL∗)2 −rµkξ− r(2H∗ + E∗)ϕH∗αξ s(1 + γL∗)2 − r(2H∗ + E∗)ϕH∗αξ s(1 + γL∗)2 + rH2∗ϕµξ s(1 + γL∗) , B1 = −bE∗β, B2 = rbE∗β−bE∗βµ−r(2H ∗ + E∗)bE∗β s , B3 = ϕrH∗βvL∗ (1 + γL∗) + rbE∗βµ − r(2H ∗ + E∗)bE∗βµ s + ϕH∗βvL∗ (1 + γL∗) . (14) In the presence of τ , determining root signs was challenging. We knew that Ē was locally asymptotically stable if all its characteristic equation roots had negative real parts. When a root had a positive real part, the system became unstable. If all roots were imaginary, stability switches occurred. The transcendental Eq. (13) possessed an infinite number of B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 10 of 28 complex roots, in which the coefficients B1, B2, and B3 arose from the delay-dependent part of the characteristic equation and were multiplied by the exponential term e−στ , explicitly capturing how τ influenced stability through delayed reproduction, mortality, and maturation effects. Biologically, B1 reflected the influence of delayed reproduction on system inertia, B2 described how delayed feedback altered mortality and transition rates, and B3 captured the direct impact of delayed maturation and reproduction on growth dynamics. Large absolute values of these coefficients indicated greater sensitivity of the coexistence equilibrium to delays, potentially leading to destabilization and sustained oscillations when τ exceeded critical thresholds. To handle the infinite number of roots arising from the exponential terms, we followed the standard approach for delay differential equations: we first set the delay to zero to establish baseline stability, then considered purely imaginary roots (σ = iω) by separating the real and imaginary parts of Eq.(13) to solve for ω and τ . This allowed us to determine the critical delay values at which a pair of complex conjugate roots crossed the imaginary axis, indicating a Hopf bifurcation, and to deduce stability for τ below and above these thresholds without computing the full infinite spectrum. Case I: when τ = 0 Theorem 3. The coexistence equilibrium point Ē is stable if the following conditions are satisfied, where the relevant quantities are defined in the proof below: A4 > 0, A1A2 −A3 > 0, (A1A2 −A3)A3 −A2 1A4 > 0. (15) Proof. Without the delay, the characteristic Eq.(13) simplifies to: σ4 +A1σ 3 + (A2 +B1)σ 2 + (A3 +B2)σ + (A4 +B3) = 0. (16) Eq.(16) with coefficients in Eq.(14) has negative or imaginary roots with negative real parts, provided the Routh–Hurwitz conditions in Eq.(15) are satisfied, making the system stable. Case II: when τ > 0 In this case, Eq.(13) possesses infinitely many roots. Stability requires roots of the characteristic Eq.(13) with negative real parts. Stability changes of Ē happens when the characteristic Eq.(13) possesses purely imaginary solutions. Consider a root of Eq.(13) to be iω. So we get ω4 −A2ω 2 +A4 = [ −ω2B1 −B3 ] cosωτ − [ωB2] sinωτ, (17) A1ω 3 −A3ω = [ −ω2B1 −B3 ] sinωτ − [ωB2] cosωτ. (18) Squaring and adding the two above mentioned equations ω8 + ζ1ω 6 + ζ2ω 4 + ζ3ω 2 + ζ4 = 0. (19) B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 11 of 28 Substituting ω2 = ℓ in Eq.(19), introduced to simplify the stability analysis with respect to the delay parameter τ . This transformation reduces the degree of the characteristic equation from being in terms of ω8 (Eq. (19)) to a quartic equation in ℓ (Eq. (20)). A quartic form is easier to handle analytically and allows direct application of the Routh–Hurwitz criterion to determine stability. Additionally, since ω2 is always non-negative for real ω, working with ℓ eliminates the need to separately track positive and negative frequency roots, making it straightforward to assess stability changes and detect the onset of purely imaginary roots. ℓ4 + λ1ℓ 3 + λ2ℓ 2 + λ3ℓ+ λ4 = 0, (20) where λ1 = A2 1 − 2A2, λ2 = A2 2 + 2A4 − 2A1A3 − 2B2 1 , λ3 = A2 3 − 2A2A4 − 4B1B3 − 2B2 2 , λ4 = A2 4 + 2B2 3 . Eq.(20) exhibits roots with negative real parts if and only if the Routh-Hurwitz condition in Eq.(21) holds true. In this case, for the roots we find ω2 k = ℓk < 0, implying that the roots of Eq.(13) are real, iωk = ∓ √ ℓk, rather than purely imaginary. In summary, we have the theorem and it is stated as follows. Theorem 4. In case of the delayed model Eq.(13), the endemic steady-state Ē is locally asymptotically stable for all τ > 0, if the following conditions are satisfied: λ1 > 0, λ4 > 0, λ1λ2 − λ3 > 0, (λ1λ2 − λ3)λ3 − λ21λ4 > 0. (21) Proof. Instead, if λ4 < 0, Eq.(20) has at least one positive root ω0 > 0. We find that ±i√ω0 is a root of Eq.(13) for the delay τ∗. By Butler’s lemma given in [45], the endemic equilibrium Ē is stable when τ < τ∗, which is a critical delay value for the delayed system of Eqs.(1)-(4) remains stable. Using Eq.(17), we can determine cosωτ = 1 △ ∣∣∣∣D1ω 4 −A2ω 2 +A4 −ωB2 A1ω 3 −A3ω −ω2B1 −B3 ∣∣∣∣ , = −ω6B1 − ω4B3 +A2ω 4B1 +A2B3ω 2 − ω2A4B1 −A4B3 +ϕ4B2A1 − ω2A3B2, (22) where △ = ∣∣∣∣−B1 B2ω B2ω B1 ∣∣∣∣ , = B2 1 +B2 2ω 2 > 0. sinωτ = 1 △ ∣∣∣∣−ω2B1 −B3 ω4 −A2ω 2 +A4 −ωB2 A1ω 3 −A3ω ∣∣∣∣ , = −ω5B1A1 + ω3B1A3 −B3A1ω 3 +A3B3ω + ω5B2 B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 12 of 28 −A2B2ω 3 + ωB2A4, (23) where △ = ∣∣∣∣−ω2B1 −B3 −ωB2 −ωB2 −ω2B1 −B3 ∣∣∣∣ , = [ω2B1 +B3] 2 − ω2B2 2 . (24) τ∗ = 1 ω0 cos−1 ( NC(ω0) ∆(ω0) ) + 2πn ω0 , n = 0, 1, 2, . . . (25) where, ∆(ω) = [ ω2B1 +B3 ]2 − ω2B2 2 , NC(ω) = −ω6B1 − ω4B3 +A2ω 4B1 +A2B3ω 2 − ω2A4B1 −A4B3 + ϕ4B2A1 − ω2A3B2. The set of ordered pair is (ω0, τ ∗). Additionally, we can confirm that the following transversality requirement: d dt Reλ (τ)||τ=τ∗ = d dt η (τ) ∣∣∣∣ τ=τ∗ > 0. (26) This condition is satisfied. Owing to continuity, the real part of λ (τ) becomes positive once τ > τ∗ resulting in an unstable steady state. This indicates a bifurcation occurring at τ = τ∗, as demonstrated by Eqs.(1)−(4). 0 200 400 600 800 1000 Time(t) days 0 100 200 300 400 500 600 700 H ea lt h y T re es H (t ) Decreases = 0.0001, 0.0002, 0.0003, 0.0005, 0.001 (a) 0 200 400 600 800 1000 Time(t) days 20 25 30 35 40 45 50 H ea lt h y T re e P o p u la ti o n H (t ) r = 0.0001, 0.0002, 0.0003, 0.0004, 0.001 Increases (b) Figure 2: Healthy tree profile H(t) versus time (1000 days) in the presence of a time delay τ = 0.5 (a) for various values of contact rate ϕ (b) for various values of replanting rate r with other parameters held at their fixed values. B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 13 of 28 0 200 400 600 800 1000 Time(t) days 0 20 40 60 80 In fe ct ed T re e P o p u la ti o n E (t ) = 0.0001, 0.00015, 0.0002, 0.0003, 0.0005 Increases (a) 0 200 400 600 800 1000 Time (t) days 0 10 20 30 40 50 60 In fe ct ed T re e P o p u la ti o n E (t ) Decreases = 0.06, 0.07, 0.08, 0.09, 0.1 (b) Figure 3: Infected tree profile E(t) versus time (1000 days) in the presence of a time delay τ = 0.5 (a) for various values of contact rate ϕ, (b) for various values of death rate of whitefly µ with other parameters held at their fixed values. 0 200 400 600 800 1000 Time (t) days 5 6 7 8 9 10 11 12 In fe ct ed T re e P o p u la ti o n E (t ) Increases = 0.1, 5, 15, 25, 35 (a) Figure 4: Infected tree profile E(t) versus time (1000 days) in the presence of a time delay τ = 0.5 (a) for various values of time delay parameter τ with other parameters held at their fixed values. B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 14 of 28 0 200 400 600 800 1000 Time(t) days 0 10 20 30 40 50 W h it ef ly P o p u la ti o n L (t ) Increases = 0.1, 0.15, 0.2, 0.25, 0.3 (a) 0 200 400 600 800 1000 Time(t) days 5 10 15 20 25 30 W h it ef ly P o p u la ti o n L (t ) Decreases = 0.06, 0.07, 0.08, 0.09, 1.00 (b) 0 200 400 600 800 1000 Time(t) days 5 10 15 20 25 30 35 W h it ef ly P o p u la ti o n L (t ) Increases = 5, 7, 10, 13, 15 (c) Figure 5: Whitefly profile L(t) versus time (1000 days) in the presence of a time delay τ = 0.5 (a) for various values of its birth rate α, (b) for various values of death rate µ, (c) for various values of time delay parameter τ with other parameters held at their fixed values. B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 15 of 28 0 200 400 600 800 1000 Time (t) (days) 0 20 40 60 80 100 120 A w a rn es s p ro g ra m s A (t ) = 0.02, 0.021, 0.023, 0.025, 0.027 Increases (a) 0 200 400 600 800 1000 Time (t) (days) 0 20 40 60 80 100 A w a rn es s p ro g ra m s A (t ) = 0.01, 0.011, 0.013, 0.015, 0.017 Decreases (b) 0 200 400 600 800 1000 Time (t) (days) 0 20 40 60 80 100 A w a rn es s p ro g ra m s A (t ) Increases = 1, 5, 10, 15, 20 (c) Figure 6: Awareness programs profile A(t) versus time (1000 days) in the presence of a time delay τ = 0.5 (a) for various values of its local awareness rate β, (b) for various values of fading away rate ξ, (c) for various values of time parameter τ with other parameters held at their fixed values. 3. SDE model ODE models have limitations, as real biological systems are often influenced by factors that are difficult to predict or quantify. Therefore, it is necessary to enhance these models to include more complex models that can capture the variability in biological dynamics. Consequently, understanding randomness as a key factor in the evolution of biological systems under stochastic influences is crucial. Here, we built an SDE model by B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 16 of 28 expanding the ODE model with the addition of stochastic processes to the derived system equations.The text below outlines our approach to deriving an SDE model from the ODE model, using the method introduced by Yuvan et al. [46]. Let X(t) = (X1(t), X2(t), X3(t), X4(t)) T be a continuous random variable corresponding to (H(t), E(t), L(t), A(t))T , where T denotes the transpose of a matrix. Additionally, let represent the random vector capturing the changes in the random variables over the time interval △t. Every possible change between states in the SDE model is represented by the transition map. Based on ODE model in [40], there are 13 possible state changes that can occur within a small time interval △t. Table 2 provides an overview of the state changes along with their associated probabilities. As an example, if one tree gets infected by whiteflies, the corresponding state change △X is expressed as △X = (−1, 1, 0, 0) and the probability of this change can be calculated using Table 2. △X = X(t+△t)−X(t) = Table 2: Possible state changes and their corresponding probabilities. Possible State Change Probability of State Change (∆X)1 = (−1, 1, 0, 0)T : Healthy trees interact with whiteflies and turn into infected trees P1 = ϕX1X3△t+ o(△t) (∆X)2 = (0,−1, 0, 0)T : Natural death of infected trees P2 = kX2△t+ o(△t) (∆X)3 = (0, 0, 1, 0)T : Spread of whiteflies from infected trees P3 = αX2△t+ o(△t) (∆X)4 = (0,−1, 0, 0)T : Eradication of infected trees by awareness programs P4 = bX4X2△t+ o(△t) (∆X)5 = (−1, 0, 0, 0)T : Eradication of healthy trees by awareness programs P5 = γX1X2 3△t+ o(△t) (∆X)6 = (0, 0, 0, 1)T : Awareness programs initiated by uninfected trees P6 = δ△t+ o(△t) (∆X)7 = (0, 0, 0, 1)T : Awareness programs initiated by infected trees P7 = βX2△t+ o(△t) (∆X)8 = (0, 0,−1, 0)T : Natural death of whiteflies P8 = µX3△t+ o(△t) (∆X)9 = (0, 0, 0,−1)T : Awareness programs fade over time P9 = ξX4△t+ o(△t) (∆X)10 = (1, 0, 0, 0)T : Growth of healthy trees (logistic model) P10 = rX1△t+ o(△t) (∆X)11 = (−1, 0, 0, 0)T : Death of healthy trees (logistic model) P11 = ( rX1(X1+X2) s ) △t+ o(△t) (∆X)12 = (1, 0, 0,−1)T : Death of whiteflies due to awareness programs P12 = µX3△t+ o(△t) (∆X)13 = (0, 0, 0, 0)T : No change in state P13 = 1− ∑12 i=1 Pi + o(△t) (△X1,△X2,△X3,△X4) Prob{ (△X1, △X2, △X3, △X4) =(-1,1,0,0)| (X1, X2, X3, X4) }=P1=ϕ HW△ t+o(△ t). Based on the stochastic modeling technique introduced by Allen et al. [47], the system of stochastic model equations is derived and given by the following expressions. dX⃗ = f⃗(t, X⃗(t))dt+B(t, X⃗(t))dW⃗ (t) X⃗(0) = [X1(0), X2(0), X3(0), X4(0)] T } . B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 17 of 28 Here we define the drift vector as f⃗ = 13∑ j=1 Pj λ⃗j . (27) Here, λ⃗j represents the random changes, and Pj denotes the corresponding transition probabilities (refer to Table 2). The drift vector f⃗ is given by f⃗ = P1λ⃗1 + P2λ⃗2 + P3λ⃗3 + P4λ⃗4 + P5λ⃗5 + P6λ⃗6 + P7λ⃗7 + P8λ⃗8 + P9λ⃗9 + P10λ⃗10 + P11λ⃗11 + P12λ⃗12, (28) f⃗ =  rX1 ( 1− X1+X2 s ) − ϕX1X2 1+γX3 ϕX1X3 1+γX3 − kX1 − bX4 αX1 − µX3 − vX4X3 δ + βX1 − ξX4  . (29) We shall derive the covariance matrix which is defined as B⃗ =  V11 V12 0 V14 V21 V22 0 0 0 0 V33 0 V41 0 0 V44  , (30) where, V11 = P1 + P5 + P10 + P11 + P8 = ϕX1X3 + γX1X 2 3 + rX1 + rX1(X1 +X2) s + µX3, V22 = P1 + P2 + P4 = ϕX1X3 + kX2 + bX4X2, V33 = P3 + P8 = αX2 + µX3, V44 = P6 + P7 + P9 = δ + βX2 + ξX4, V12 = V21 = −P1 = −ϕX1X3, V14 = V41 = −P8 = −µX3. The square root of this covariance matrix must be computed. However, for an n × n positive semi - definite matrix with order greater than two, there is no explicit formula for calculating its square root. Therefore, in order to proceed, we will implement the second modeling procedure, which generates a diffusion matrix for which a square root may not be required. Following the second modeling approach developed by Allen et al. [47], the corresponding stochastic model is presented as follows. dX⃗ = f⃗(t, X⃗(t)dt+G(t, X⃗(t))dW⃗ (t), (31) B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 18 of 28 X⃗(0) = [X1(0), X2(0), X3(0), X4(0)] T . (32) The drift vector f⃗ is the same as that obtained in the first modeling procedure. The diffusion matrix G is defined by G = λi,jP 1/2 j , j = 1, 2, ........, 12, i = 1, 2, ......, 4 (33) G =  − √ ϕX1X3 0 0 0 − √ γX1X2 3 0 0 0 0 √ rX1 − √ rX1(X1 +X2) s √ µX3 √ ϕX1X3 − √ kX2 0 − √ bX4X2 0 0 0 0 0 0 0 0 0 0 √ αX2 0 0 0 0 − √ µX3 0 0 0 0 0 0 0 0 0 √ δ √ βX2 0 − √ ξX4 0 0 − √ µX3  In the SDE formulation, the diffusion matrixG encodes the stochastic fluctuations associated with each biological transition in the system. Each column of G corresponds to a specific reaction or event (e.g., infection, recovery, natural death, recruitment), while the non-zero entries indicate how that event changes the state variables. The square-root terms arise from the standard deviation of the Poisson process governing each transition rate. For example, the term − √ ϕX1X3 in the first row represents the random decrease in susceptibles due to infection, while the corresponding + √ ϕX1X3 in the second row represents the simultaneous stochastic increase in exposed individuals from the same event. In this way, G directly links the random noise terms to the underlying biological processes. Therefore, the Itô stochastic differential model is expressed in the following way: dX(t) = f(X1, X2, X3, X4) dt+GdW (t), (34) which can also be written as dX(t) = f(t,X(t)) dt+B(t,X(t)) dW (t), with initial conditions X(0) = ( X1(0), X2(0), X3(0), X4(0) )T , and W (t) = ( W1(t),W2(t),W3(t),W4(t) )T is denoting a vector of independent Wiener processes (standard Brownian motions) representing stochastic fluctuations. B(t,X(t)) is the diffusion matrix describing how these stochastic effects influence the state variables. Therefore, the SDE model is subsequently constructed as follows: dH = ( rH ( 1− H + E s ) − ϕHL 1 + γL ) dt− √ ϕX1X3dW1 − √ γX1X2 3dW5 + √ rX1dW10 − √ (rX1(X1 +X2)/s)dW11 − √ µX3dW12, (35) B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 19 of 28 dE = ( ϕHL 1 + γL − kE − bAE ) dt+ √ ϕX1X3dW1− √ kX2dW2− √ bX4X2dW4, (36) dL = (αE − µL− vAL) dt+ √ αdW3− √ µX3dW8, (37) dA = (δ + βE − ξA) dt+ √ δdW6+ √ δdW7− √ ξX4dW9− √ µX3dW12. (38) 0 200 400 600 800 1000 Time (t) (days) 20 25 30 35 40 45 50 55 H ea lt h y T re es H (t ) SDE DDE (a) 0 200 400 600 800 1000 Time (t) (days) 0 5 10 15 20 In fe ct ed T re e P o p u la ti o n E (t ) SDE DDE (b) Figure 7: The plot illustrates the variation that investigates the difference between DDE and SDE of (a) healthy tree population, (b) infected tree population. 0 200 400 600 800 1000 Time (t) (days) 5 10 15 20 25 30 35 40 W h it ef ly p o p u la ti o n L (t ) SDE DDE (a) 0 200 400 600 800 1000 Time (t) (days) 0 5 10 15 20 25 30 35 A w a re n es s P ro g ra m s A (t ) SDE DDE (b) Figure 8: The plot illustrates the variation that investigates the difference between DDE and SDE of (a) whitefly population, (b) awareness programs. B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 20 of 28 4. Sensitivity analysis This Section discusses the impact of varying parameter values on the functional value of the reproduction number R0. Since the critical parameter may act as a critical threshold for illness therapy, it must be determined [48]. The sensitivity indices of R0 with respect to the parameters ϕ, s, α, ξ, k, b, δ, µ, v are represented algebraically as follows: ∂R0 ∂ϕ = sαξ2 (ξk + bδ) (µξ + vδ) , ∂R0 ∂s = ϕαξ2 (ξk + bδ) (µξ + vδ) , ∂R0 ∂α = ϕsξ2 (ξk + bδ) (µξ + vδ) , ∂R0 ∂ξ = 2ϕsαξ (ξk + bδ) (µξ + vδ) − ϕsαξ2k (ξk + bδ)2 (µξ + vδ) − ϕsαξ2µ (ξk + bδ) (µξ + vδ)2 , ∂R0 ∂k = −ϕsαξ3 (ξk + bδ)2 (µξ + vδ) , ∂R0 ∂b = −ϕsαξ2δ (ξk + bδ)2 (µξ + vδ) , ∂R0 ∂δ = −ϕsαξ2b (ξk + bδ)2 (µξ + vδ) − ϕsαξ2v (ξk + bδ) (µξ + vδ)2 , ∂R0 ∂µ = −ϕsαξ3 (ξk + bδ) (µξ + vδ)2 , ∂R0 ∂v = −ϕsαξ2δ (ξk + bδ) (µξ + vδ)2 . The analysis concludes that certain partial derivatives are positive, and that increasing any of the positive parameters ϕ, s, α, ξ causes an increase in the basic reproductive number R0. Elasticity is estimated by assessing the proportional response to proportional perturbations. We’ve Eϕ = ϕ R0 ∂R0 ∂ϕ = ( ϕ(ξk + bδ)(µξ + vδ) ϕsαξ2 )( sαξ2 (ξk + bδ) (µξ + vδ) ) = 1, Es = s R0 ∂R0 ∂s = ( s(ξk + bδ)(µξ + vδ) ϕsαξ2 )( ϕαξ2 (ξk + bδ) (µξ + vδ) ) = 1, Eα = α R0 ∂R0 ∂α = ( α(ξk + bδ)(µξ + vδ) ϕsαξ2 )( ϕsξ2 (ξk + bδ) (µξ + vδ) ) = 1, Eξ = ξ R0 ∂R0 ∂ξ = ξ(ξk + bδ)(µξ + vδ) ϕsαξ2 · ϕsαξ (ξk + bδ) (µξ + vδ) ( 2− ξk ξk + bδ − ξµ µξ + vδ ) = 0.6428. Eϕ, Es, Eα, and Eξ are all positive, as evidenced by the above expressions. This indicates that an increase in the values of the parameters ϕ, s, α, ξ leads to an increase in the basic reproduction number R0. The fundamental reproduction number can experience B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 21 of 28 significant variation from even the smallest alterations in these parameters. It is essential to precisely calculate sensitive parameters, as minor changes can cause major quantitative changes in the system. Figure 9: Coconut trees Figure 10: Life cycle of Rugose Spiralling Whitefly B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 22 of 28 1 2 3 4 5 6 7 8 9 10 (Contact Rate) 10 -4 50 100 150 200 250 300 350 400 450 500 550 R 0 ( R e p ro d u c ti o n N u m b e r) (a) 50 55 60 65 70 75 80 85 90 95 100 s (Tree Density) 350 400 450 500 550 600 650 700 750 R 0 ( R e p ro d u c ti o n N u m b e r) (b) Figure 11: The plot illustrates reproduction number corresponding to the sensitive parameter (a) ϕ, (b) s. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 (Whitefly Birth Rate) 0 500 1000 1500 2000 2500 3000 3500 4000 R 0 ( R e p ro d u c ti o n N u m b e r) (a) 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 (Saturation Constant) 2500 3000 3500 4000 4500 5000 5500 6000 6500 R 0 ( R e p ro d u c ti o n N u m b e r) (b) Figure 12: The plot illustrates reproduction number corresponding to the sensitive parameter (a) α, (b) ξ. B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 23 of 28 (a) (b) Figure 13: Surface plot illustrating the reproduction number R0 in relation to (a) ξ and α, (b) s and ϕ. 5. Results and discussion We have selected the values in a manner that allows numerical analysis to be performed using a reference point for each parameter. We consider an initial population of 50 healthy trees, 5 infected trees, 10 whiteflies per tree, and 10 awareness programs. So, we start with the following initial values for H(0) = 50, E(0) = 5, L(0) = 10, and A(0) = 10. Figure 1 displays the schematic diagram of our model. Figure 2a shows the variation of H(t) with different contact rates ϕ while keeping τ = 0.5 fixed. As ϕ increased, the rate at which healthy trees became infected, leading to a more rapid decline in the healthy tree population. Figure 2b depicts the variation of H(t) with different replanting rates r. Higher values of r counteracted the loss of healthy trees by replacing infected ones, resulting in a slower decline in H(t). Figure 3a shows the effect of varying the contact rate ϕ on E(t). An increase in ϕ accelerated the infection process, causing E(t) to rise to higher peak values before stabilizing. Figure 3b shows the influence of the whitefly death rate µ on E(t). Larger µ values reduced the number of vectors, thereby lowering the infected tree population over time. Figure 4a demonstrates the effect of changing the time delay τ on E(t). Increasing τ slightly delayed the peak of infection, with larger delays associated with a marginally higher infected tree population during the peak period. Figure 5a shows the variation in L(t) for different whitefly birth rates α. Higher α values led to more rapid growth of the whitefly population and larger peak values. Figure 5b illustrates the effect of increasing the whitefly death rate µ. Larger µ values reduced both the peak and steady-state levels of L(t). Figure 5c examines the influence of time delay τ , showing that longer delays caused higher peaks in L(t) before reaching equilibrium. Figure 6a shows the impact of increasing the local B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 24 of 28 awareness rate β. Larger β values led to faster and higher growth of awareness programs. Figure 6b depicts the effect of the fading rate ξ. Higher ξ values decreased the maximum awareness level and accelerated its decline over time. Figure 6c presents the variation with different time delays τ , showing that increasing τ resulted in higher peaks in A(t) before settling. Figure 7a compares H(t) for DDE and SDE models. Both showed similar overall trends, but the SDE model captured small stochastic fluctuations around the deterministic trajectory. Figure 7b compares E(t) under DDE and SDE, where the SDE results exhibited higher variability due to random effects. Figure 8a shows L(t) for both models, with stochastic variation being more pronounced compared to H(t) and E(t). Figure 8b compares A(t) between DDE and SDE, revealing noticeable fluctuations in the stochastic case. The benefit of comparing the DDE and SDE formulations was that it highlighted how stochastic fluctuations influenced system dynamics, providing deeper insights into variability, uncertainty, and the robustness of the model predictions. Figure 9 displays the picture of coconut trees taken in a coconut plantation near Pollachi, Tamil Nadu. Figure 10 displays the life cycle of rugose spiralling whitefly. In Figure 11a, illustrates when ϕ increases, the reproduction number grows proportionally, indicating that higher contact rates between trees and whiteflies result in greater potential for whitefly infestation spread, 11b illustrates the density of trees s increases, the reproduction number also rises, indicating that higher tree densities facilitate the spread of whitefly infestations. Figure 12a it is inferred that the whitefly birth rate α rises, the basic reproduction number grows significantly, highlighting the strong influence of whitefly reproduction on the potential for infestation spread, 12b displays the higher fading rate of awareness ξ significantly influences the potential for disease or infestation spread, as a greater R0 indicates a stronger possibility of an outbreak or epidemic. The graph shows the critical role of awareness decay in controlling the spread within a population. Figure 13a shows that R0 increases with both higher contact rates ϕ and greater tree density s, indicating that a higher density of trees, combined with more frequent interactions, can lead to an increase in the spread of an infestation. The color bar on the right represents the corresponding values of R0, with warmer colors (yellow and green) denoting higher reproduction numbers. This emphasizes the role of environmental and contact factors in influencing potential outbreak dynamics. Figure 13b shows that as α increases, R0 also rises, particularly when ξ is small. This indicates that a higher whitefly birth rate α can significantly contribute to the spread of an infestation. However, when ξ increases, representing a higher rate of awareness or intervention, the value of R0 is mitigated, underscoring the importance of awareness in controlling the spread. The color bar on the right highlights the value of R0, with warmer colors (yellow and green) signifying higher reproduction numbers. Table 1 shows the parameter values used for the analysis and Table 2 shows the various state changes along with their respective probabilities. 6. Conclusion This study developed a comprehensive mathematical model using both DDEs and SDEs to analyze the complex dynamics of whitefly infestations in coconut farming. The B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 25 of 28 DDE model captured the delayed effects of awareness programs and, through equilibrium and stability analysis, demonstrated that the time delay parameter τ played a critical role in determining whether the infestation persisted or was eliminated. Numerical simulations and sensitivity analyses further showed that even moderate delays in implementing awareness programs significantly affected outbreak control. The inclusion of the SDE model addressed environmental and demographic variability, offering a more realistic representation of the system’s behavior under uncertainty. This integrated approach strengthened the predictive capabilities of the model and supported the development of more resilient pest management strategies. By bridging mathematical analysis with ecological dynamics and farming interventions, this work offered detailed insights into the timing and effectiveness of awareness-based controls for pest outbreaks. By capturing both delayed responses and environmental variability, the study paved the way for more resilient and adaptive pest management. Future research can extend this model by incorporating spatial dynamics to capture the geographical spread of infestations across plantations, as well as economic factors to optimize resource allocation for control strategies. Coupling the model with economic optimization could support cost-effective intervention planning. The inclusion of multiple pest species and natural predator–prey interactions would improve the biological realism of the framework. Furthermore, integrating real-time field data through data assimilation techniques could allow the development of adaptive, evidence-based policies that respond dynamically to emerging infestations. Such advancements would enhance the model’s applicability for decision-making in sustainable agricultural pest management. Acknowledgements We sincerely thank the reviewers for their valuable comments, which were of great help in revising the manuscript. The authors are pleased to acknowledge the financial support of the Selective Excellence Research Initiative (SRMIST/R/AR(A)/SERI2023/174/31). It is our pleasure to thank the College of Engineering and Technology, SRM IST for its valuable support and constant encouragement. Credit Authorship Contribution Statement B. Dhivyadharshini: Investigation, Methodology, Software, Writing – original draft, Visualization. R. Senthamarai: Conceptualization, Methodology, Validation, Resources, Writing – review & editing, Supervision. Conflict of Interest The authors declare that there is no conflict of interest. B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 26 of 28 References [1] C.C. Sreejith, C. Muraleedharan, P. Arun, Life cycle assessment of producer gas derived from coconut shell and its comparison with coal gas: An Indian perspective. Int. J. Energy Environ. Eng., 4: 1–22, 2013. [2] A. Snehalatharani, H. P. Maheswarappa, V. Devappa, S.k. Malhotra. Status of coconut basal stem rot disease in India - A review. Indian J. Agric. Sci., 86: 1519–1529, 2016. [3] A. Varghese, J. Jacob. A study of physical and mechanical properties of the Indian coconut for efficient dehusking. Journal of Natural Fibers, 14(3), 390-399, 2017. [4] F. M. Dayrit, T. N. Mary, The Potential of coconut oil and its derivatives as effective and safe antiviral agents against the novel coronavirus. Indian Coconut Journal, 62: 21-23, 2020. [5] G. Suganya and R. Senthamarai. Analytical Approximation of a Nonlinear Model for Pest Control in Coconut Trees by the Homotopy Analysis Method. Computer Research and Modeling, 14(5): 1093–1106, 2022. [6] B. Dhivyadharshini and R. Senthamarai. Modeling the indirect impact of rhinoceros beetle control on red palm weevils in coconut plantations. Computer Research and Modeling, 17(4): 737 − 752, 2025. [7] D. Narayana, K. N. Nair, Trends in area, production and productivity of coconuts in Kerala. Indian Journal of Agricultural Economics, 44(902-2018-2692): 159-167, 1989. [8] k. Nihad, A. Haris, S. Kalavathi, Scope of floriculture in coconut garden. Indian Coconut Journal, 5-8, 2020. [9] J.U. Chikaire, J.O. Ajaero, C.N. Atoma, Socio-economic Effects of Covid-19 pandemic on rural farm families’ well-being and food systems in Imo State, Nigeria. Journal of Sustainability and Environmental Management, 1(1), 18-21. 2022. [10] M. Abad, P. Noguera, R. Puchades, A. Maquieira, V. Noguera Physico-chemical and chemical properties of some coconut coir dusts for use as a peat substitute for containerised ornamental plants. Biores. Technol., 82:241-245, 2002. [11] O. A. Carrijo, R.S. Liz, N. Makishima. Green coconut husk fiber as an agricultural substrate. Hort. Bras., 20:533-535, 2002. [12] M.U.C. Nunes. Coconut husk fiber and dust: products of great importance for industry and agriculture. In: Aragão WM (Ed.) Coconut: Post-Harvest, Embrapa Information Technology, Brasília (Brazilian Fruit Series), 29:66-71, 2002. [13] J.E.G. van Dam, M.J.A. van den Oever, E.R.P, Keijsers. Production process for high density high performance binderless boards from whole coconut husk. Ind. Crops Prod., 20:97-101, 2004. [14] S. Nadanasabapathy, R. Kumar. Physico-chemical constituents of tender coconut (Cocos nucifera) water. The Indian J. Agric. Sci., 69: 750-51, 2013. [15] J.H. Martin Zootaxa, 681: 1-119, 2004. [16] K.S. Karthick, C. Chinniah, P. Parthiban, A. Ravikumar. Newer report of Rugose Spiraling Whitefly, Aleurodicus rugioperculatus Martin (Hemiptera: Aleyrodidae) in India. International Journal of Research Studies in Zoology 4(2), 2018, 12-16. B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 27 of 28 [17] K. Elango, S. Jeyarajan Nelson, S. Sridharan, V. Paranidharan and S. Balakrishnan. Biology, Distribution and host range of new invasive pest of India coconut rugose spiralling whitefly aleurodicus rugioperculatus martin in Tamil Nadu and the status of its natural enemies. International Journal of Agriculture Sciences, 11(9): 8423-8426, 2019. [18] G. Suganya, E. Jenitta and R. Senthamarai. A study on the dynamics of pest population with biocontrol using predator, parasite in presence of awareness. Computer Research and Modeling, 16(3): 713-729, 2024. [19] A. Josephrajkumar, C. Mohan, V. Krishnakumar. Parasitism induced bio-suppression of coconut whitefly in Kerala. Kerala Karshakan e-journal, 26-27, 2016. [20] C. Mohan, A. Josephrajkumar, V. Hegde, V. Krishnakumar, P.B. Renjith, A.S. Anjali and P. Chowdappa. Gradient outbreak and bio- suppression of spiralling whitefly in coconut gardens in South India. Indian Coconut Journal, 59(8): 9-12, 2016. [21] L.J. Allen, F. Brauer, P. Van den Driessche and J. Wu, Mathematical epidemiology. Berlin: Springer, 1945: 2019. [22] M. Sivakumar and R. Senthamarai. Mathematical model of epidemics: SEIR model by using homotopy perturbation method. AIP Conference Proceedings, 2112(1), 2019. [23] B. Dhivyadharshini and R. Senthamarai. Mathematical Analysis of a Non Linear Prey Predator System: Analytical Approach By HPM. AIP Conference Proceedings, 2516, 2022. [24] G. Suganya and R. Senthamarai, Mathematical modeling and analysis of the effect of the rugose spiraling whitefly on coconut trees. AIMS Mathematics, 7(7), 13053-13073, 2022. [25] F.A. Basir, A. Banerjee and S. Ray, Role of farming awareness in crop pest management-A mathematical model. Journal of theoretical biology, 461, 59-67, 2019. [26] A. Sharma and A.K. Misra, Modeling the impact of awareness created by media campaigns on vaccination coverage in a variable population. Journal of biological systems, 22(02): 249-270, 2014. [27] F. Al Basir, E. Venturino, S. Ray, P. K. Roy, Effects of awareness program for controlling mosaic disease in Jatropha curcas plantations. Computational and Applied Mathematics, 37:6108–6131, 2018. [28] P. Van den Driessche and J. Watmough. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical biosciences, 180(1-2): 29-48, 2002. [29] M. Sivakumar and R. Senthamarai. Mathematical model of epidemics: Analytical approach to SIRW model using homotopy perturbation method. AIP Conference Proceedings, 2277, (2020). [30] T. Vijayalakshmi and R. Senthamarai. An analytical approach to top predator interference on the dynamics of a food chain model. Journal of Physics: Conference Series, 1000(1), 2018. [31] T. Vijayalakshmi and R. Senthamarai. Application of homotopy perturbation and variational iteration methods for nonlinear imprecise prey-predator model with stability analysis. The Journal of Supercomputing, 78(2): 2477-2502, 2022. B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (4) (2025), 6691 28 of 28 [32] Y. Kuang. Delay Differential Equations with Applications in Population Dynamics. Academic Press, Inc., New York, 1993. [33] R.V. Culshaw and S. Ruan. A delay-differential equation model of HIV infection of CD4+ T-cells. Mathematical Biosciences, 165:27–39, 2000. [34] R.A. Umana, A. Omame, S.C. Inyama, Deterministic and Stochastic Models of the Dynamics of Drug Resistant Tuberculosis. FUTOJNLS, 2(2): 173-194, 2016. [35] S. Boulaaras, S. I. Araz, A. Alharbi. Radiotherapy Effects in Tumor Treatment: A Piecewise Model with Deterministic and Stochastic Approaches. Fractals, 16(44), 2025. [36] S. I. Araz, M. A. Cetin, A. Atangana. Existence, uniqueness and numerical solution of stochastic fractional differential equations with integer and non-integer orders. Electronic Research Archive, 32(2): 733–761, 2024. [37] İ. A. Arık, S. İ. Araz. Crossover behaviors via piecewise concept: A model of tumor growth and its response to radiotherapy. Results in Physics, 41: 105894, 2022. [38] K. S. Kim, S. Kim, I. H. Jung, Dynamics of tumor virotherapy: A deterministic and stochastic model approach. Stochastic Analysis and Applications, 34(3): 483-495, 2016. [39] F.A. Basir, E. Venturino, S. Ray and P.K. Roy, Impact of farming awareness and delay on the dynamics of mosaic disease in Jatropha curcas plantations. Computational and Applied Mathematics, 37(5), 6108-6131, 2018. [40] G. Suganya and R. Senthamarai. Impact of Awareness on the Dynamics of Pest Control in Coconut Trees - A Mathematical Model, Engineering Letters, 30:4, 30(4) 2022. [41] J. Hale. Theory of functional differential equations. Springer, Heidelberg, 1977. [42] M. Bodnar. The nonnegativity of solutions of delay differential equations. Appl Math Lett 13(6):91–5, 2000. [43] X. Yang, L. Chen, J. Chen. Permanence and positive periodic solution for the single species nonautonomus delay diffusive model. Comput Math Appl, 32:109–116, 1996. [44] O. Diekmann, Heesterbeek JAP, and Metz JAJ. On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. Journal of Mathematical Biology, 28(4):365–382, 1990. [45] H.I. Freedman, V.S.H Rao, The trade-off between mutual interference and time lags in predator–prey systems. Bull Math Biol, 45(6):991–1004, 1983. [46] Y. Yuan, L.J.S. Allen, Stochastic models for virus and immune system dynamics. Math. Biosci., 234:84–94, 2011. [47] E. J. Allen, L.J.S. Allen, A. Arciniega, P. E. Greenwood, Construction of equivalent stochastic differential equation models. Stoch. Anal. Appl., 26(2): 274 – 297, 2008. [48] B. Dhivyadharshini, R. Senthamarai. Modeling Rugose Spiraling Whitefly Infestation on Coconut Trees Using Delay Differential Equations: Analysis via HPM. European Journal of Pure and Applied Mathematics, 17(3): 1908 – 1936, 2024.