Advanced Formulation of Ecological Biopesticides Based on Citrus Limonum in Clayey Matrices: Optimization of Diffusive Dynamics

1. Introduction

The preservation of stored grains is a critical global challenge, essential for ensuring food security amid increasing demand. Pests like Sitophilus granarius (granary weevil) pose significant threats by infesting and degrading grain reserves, leading to both quantity and quality losses. Traditional management of these pests has relied heavily on synthetic chemical insecticides, which, despite their effectiveness, have led to issues such as pest resistance, health risks from chemical residues, and environmental damage [1,2,3,4,5]. Consequently, there is a growing need for safer, more sustainable pest control alternatives. Essential oils, particularly lemon essential oil, have emerged as promising candidates due to their insecticidal properties, low toxicity, and environmental safety [6]. Lemon essential oil, rich in limonene, is effective against Sitophilus granarius and is biodegradable, reducing long-term ecological risks [7,8,9]. This study aims to investigate the potential of lemon essential oil as a natural insecticide, focusing on its interaction with porous clay supports to enhance its residual activity and effectiveness. The study will evaluate the oil’s efficacy under various conditions using advanced statistical and modeling techniques to optimize its application in real-world storage environments [10]. The research also explores the use of natural porous clays from Morocco to encapsulate the oil, improving its stability and controlled release, which are critical for maximizing insecticidal efficacy [11,12,13]. The clays will be characterized using advanced analytical techniques, and the oil's diffusion through the clay matrix will be modeled to better understand and optimize its release [14,15]. This integrated approach, combining material characterization, theoretical modeling, and bioinformatics analysis, could lead to innovative pest control solutions and a deeper understanding of the mechanisms involved in the interactions between essential oils and porous matrices, paving the way for new strategies in natural product-based insect control.

2. Materials and Methods

2.1. Materials

2.1.1. Porous Clay

The porous clay used in this study was sourced from the Bejaad region of Morocco (32°47'01.7"N 6°13'52.3"W) and characterized using several advanced techniques. X-ray diffraction (XRD) was performed with a Bruker CCD-Apex device, while the oxide content was analyzed using a Siemens SRX 3000 X-ray fluorescence spectrometer (XRF) [16]. Fourier-transform infrared spectroscopy (FT-IR) was conducted using a Bruker Alpha spectrometer [17], and the morphology was examined with a scanning electron microscope (SEM-EDX) [18]. Thermogravimetric analysis (TGA) was performed using a Du Pont analyzer [19].

2.1.2. Lemon Essential Oil

Lemon essential oil (Citrus limonum) was extracted through hydrodistillation and analyzed using gas chromatography-mass spectrometry (GC-MS) with a Hewlett Packard 5971A instrument [20]. The oil's molecular profile was further determined using gas chromatography with flame ionization detection (GC-FID) [21].

2.1.2. Preparation of the Clay/Essential Oil Mixture

The RC and GC clays were prepared by grinding and heating, then mixed with lemon essential oil in various concentrations. The mixture was placed in metal cylinders for further experimentation [22].

2.2. Theoretical Considerations

2.2.1. Diffusion Model

The desorption of lemon essential oil within the porous clay matrix can be modeled using Fick second law of diffusion. This law describes how the concentration of a diffusing substance changes over time and space within a medium. For this study, the general form of Fick second law, applied to a one-dimensional system, can be expressed as:

$${\displaystyle \frac{\partial C\left(x,t\right)}{\partial t}}=D{\displaystyle \frac{{\partial }^{2}C\left(x,t\right)}{\partial x^2}}$$

where:

-$C\left(x,t\right)$ is the concentration of the essential oil at position x and time t,

- D is the diffusion coefficient of the essential oil in the clay matrix,

- ${\displaystyle \frac{\partial C\left(x,t\right)}{\partial t}}$ is the rate of change of concentration with time,

- ${\displaystyle \frac{{\partial }^{2}C\left(x,t\right)}{\partial x^2}}$ is the second derivative of concentration with respect to position, representing the spatial variation in concentration.

In this model, the diffusion coefficient D is a critical parameter that depends on the physical properties of the lemon essential oil, the characteristics of the clay, and the interaction between the oil and the clay matrix [23]. The boundary conditions and initial conditions are defined based on the experimental setup, where the concentration of the oil at the surface of the clay and within the matrix at the initial time point is considered.

To accurately describe the desorption behavior, the model must account for the specific properties of lemon essential oil, such as its volatility, viscosity, and interaction with the clay pore structure. The heterogeneous nature of the porous clay, including its varying pore size distribution and surface area, may affect the diffusion process, potentially leading to non-uniform desorption rates [24].

2.2.2. Advanced Models

Non-linear diffusion models were considered to account for complex interactions and external factors like temperature and humidity, providing a more accurate representation of the oil's desorption behavior [25,26].

2.2.3. Assumptions

Key assumptions included one-dimensional diffusion, constant diffusivity, and uniform initial concentration within the clay matrix [27].

2.2.4. Model of Mass Transfer by Diffusion

The diffusion process was further described using Fick's second law in cylindrical coordinates, with simplified equations to model the mass transfer of the essential oil [28].

2.2.5. Mathematical Framework

Fick second law for diffusion in cylindrical coordinates (r, θ, z) is expressed as:

$${\displaystyle \frac{\partial C}{\partial t}}={\displaystyle \frac{\left(\frac{1}{r}\right)\partial }{\partial r}}\left(rDr{\displaystyle \frac{\partial C}{\partial r}}\right)+{\displaystyle \frac{\left(\frac{1}{r}\right)\partial }{\partial \theta }}\left({\displaystyle \frac{\frac{D\theta }{r}\partial C}{\partial \theta }}\right)+{\displaystyle \frac{\partial }{\partial z}}\left(Dz{\displaystyle \frac{\partial C}{\partial z}}\right)$$

where:

- $C\left(r,\theta ,z,t\right)$ is the concentration of EOCL at a position $\left(r,\theta ,z\right)$ and time t,- $Dr,D\theta ,Dz$ are the diffusion coefficients in the $r,\theta ,andz$ directions, respectively.

Given that the height of the cylinder is much smaller relative to its diameter and lateral surface, the diffusion can be assumed to primarily occur in the axial direction (along z). This simplifies the model to one-dimensional diffusion, reducing Fick second law to:

$${\displaystyle \frac{\partial C}{\partial t}}=Dz{\displaystyle \frac{{\partial }^{2}C}{\partial z^2}}$$

2.2.6. Boundary and Initial Conditions

The solution to the above equation requires appropriate boundary and initial conditions:

Initial Condition: At $t=0$, the concentration of EOCL is assumed to be uniform throughout the clay:

$$C\left(z,0\right)=C_0\mathrm{for}\mathrm{all}0\le z\le h$$

Boundary Conditions:

- For $t>0$ at the surface z = h, the concentration of EOCL at the boundary becomes:

$$C\left(z=h,t\right)=C_t$$

where Ct is the concentration at the boundary at time t.

2.2.7. Analytical Solution

The analytical solution to the simplified diffusion equation has been derived by Crank (1979) for these boundary conditions and is given by:

$$C\left(z,t\right)=\left({\displaystyle \frac{4}{\pi }}\right)\sum _{n=0}^{\mathrm{\infty }}\mathit{exp}\left[-Dz{\displaystyle \frac{{\left(2n+1\right)}^2{\pi }^{2}t}{h^2}}\right]\cos \left[\left(2n+1\right)\pi {\displaystyle \frac{z}{h}}\right]$$

The total mass Mt of EOCL in the clay RC at any instant t is obtained by integrating the concentration over the thickness of the material and the surface S:

$$M_t={\int }_0^{h}S·C\left(z,t\right)dz$$

Under the assumption that the desorption process tends towards equilibrium over an infinite desorption time, the analytical solution for the mass at any time t relative to the initial mass M₀ is given by:

$${\displaystyle \frac{M_t-M_{\mathrm{\infty }}}{M_0-M_{\mathrm{\infty }}}}=\left({\displaystyle \frac{8}{{\pi }^2}}\right)\sum _{n=0}^{\mathrm{\infty }}\left({\displaystyle \frac{1}{{\left(2n+1\right)}^2}}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left[-\mathit{Dz}{\displaystyle \frac{{\left(2n+1\right)}^2{\pi }^{2}t}{h^2}}\right]$$

where:- M₀ is the initial mass of EOCL absorbed,

- M_∞ is the mass of EOCL at equilibrium.

For practical purposes, the equation can be simplified by considering only the first term of the series and assuming M∞ = 0:

$$\left({\displaystyle \frac{M_t}{M_0}}\right)=\left({\displaystyle \frac{8}{{\pi }^2}}\right)\exp \left(-Dz{\displaystyle \frac{{\pi }^{2}t}{h^2}}\right)$$

Given that (8/π²) ≈ 1, the equation further simplifies to:

$${\displaystyle \frac{M_t}{M_0}}=\exp \left(-Dz{\displaystyle \frac{{\pi }^{2}t}{h^2}}\right)$$

The diffusivity Dz can be determined from the slope of a linear or logarithmic extrapolation of the simplified equation. The evaporation constant K can be approximated by :

$$K={\displaystyle \frac{{\pi }^{2}Dz}{h^2}}$$

The initial evaporation rate F₀ is obtained from the initial slope of the desorption curve as a function of time:

$${\mathrm{F}}_0=\left({\displaystyle \frac{1}{\mathrm{S}}}\right)\underset{\left(\mathrm{t}\to 0\right)}{\lim }{\displaystyle \frac{dM}{dt}}$$

These parameters, diffusivity Dz, evaporation constant K, and evaporation rate F are critical in describing the desorption process and are central to the modeling and understanding of EOCL behavior in the porous clays RC and GC [29].

2.3. Experimentation and Modeling - Essentials

2.3.1. Insecticidal Activity Test

$$M\mathrm{\%}={\displaystyle \frac{M_{test}-M_{control}}{M_{control}}}\times 100$$

where:

$M\mathit{\%}$ : is the corrected mortality,

$M_{test}$ : is the observed mortality during the test,

$M_{control}$ :​ is the mortality observed in the control.

Lethal Dose ( $\mathit{L}{\mathit{D}}_{50}$ ) Determination: The lethal dose required to kill 50% of the insect population $L{D}_{50}$ was determined by linear interpolation of curves that plot the percentage of mortality as a function of the logarithm of the tested concentration [34]. This provides a quantitative measure of the toxicity of EOCL, EOCL + RC, and EOCL + GC under the conditions tested.

This procedure allowed for a detailed assessment of EOCL insecticidal efficacy, both alone and when mixed with RC and GC clays, under controlled conditions, providing valuable data for optimizing its use in pest control.

2.3.2. Optimization of Transfer Conditions

To optimize the experimental conditions for the diffusion transfer of EOCL (lemon essential oil) in both RC and GC porous clay media, a statistical approach based on experimental designs was employed. The transfer process was evaluated by measuring the evaporation rates, denoted as F.

A full factorial design was used to structure the experiments, considering four factors ($X_i$​):

Factor 1 = C: Concentration of essential oil (0.01 mL/cm³ and 0.02 mL/cm³).

Factor 2 = T: Temperature (25 °C and 30 °C).

Factor 3 = D: Cylinder diameter (1 cm and 2 cm).

Factor 4 = M: Mass of clay rock (0.05 g and 0.10 g).

The total number of tests required was calculated using the formula:

Number of tests =$2^{\mathrm{k}}$where 2 represents the two levels (-1 and +1) and k is the number of factors studied. In this study, with $\mathrm{k}=4\mathrm{k}$Number of tests= $2^4$=16

This results in 16 tests being conducted for each type of clay (RC and GC). Table 1 presents the design matrix of the experiments and the corresponding evaporation rates obtained for both RC and GC clays. Each row in the Table 1 represents a specific test, with the factors coded according to their high (+1) or low (-1) levels [35].

The polynomial model used to describe the relationship between the factors and the evaporation rates F for each clay type is expressed as:

$$F=b_0+\sum _{i=1}^n{b}_{iXI}+\sum _{i=1}^n{\sum }_{j=1}^{n-1}{b}_{ij}{X}_i{X}_j+\sum _{i=1}^n{\sum }_{j=1}^{n-1}\sum _{k=1}^{n-2}{b}_{ijk}{X}_i{X}_j{X}_k+b_{ijk}{X}_i{X}_j{X}_k{X}_l$$

where:

$b_0$​: is the mean,

$b_i$: represents the main effects (4 main effects for the factors),

$b_{ij}$ : represents the second-order interaction effects (6 interactions),

$b_{ijk}$: represents the third-order interaction effects (4 interactions),

$b_{ijkl}$​: represents the fourth-order interaction effect (1 interaction).

This model allows for the analysis of how each factor, and their interactions, influence the evaporation rates, enabling the optimization of the transfer conditions for EOCL in both RC and GC clays [36].

2.3.3. Study of Insecticidal Activity Behavior of Essential Oil Mixture in Clay Porous Media

The interaction between the behavior of EOCL in the porous clay media (RC and GC) and its insecticidal activity was studied by analyzing the relationship between the parameters of mass transfer by diffusion and the mortality rate in insecticidal activity. This analysis was conducted using Principal Component Analysis (PCA), a powerful statistical tool for compressing and synthesizing information. PCA is particularly useful when dealing with a large amount of quantitative data, as it helps to interpret and identify correlations between different variables, providing insights into the effectiveness of EOCL in different clay matrices [37,38].

Table 2. LD₅₀ ​ of the Insecticidal Activities of EOCL Fixed on Porous Media (RC and GC).

Diameter (cm)	Temperature (°C)	LD50 (a) for RC	(b) for RC	LD50 (a) for GC	(b) for GC
D1 = 1	T1 = 25	0.228	0.110	0.229	0.115
	T2 = 30	0.350	0.109	0.349	0.111
	T3 = 35	0.035	0.017	0.036	0.018
D2 = 2	T1 = 25	0.529	0.032	0.530	0.034
	T2 = 30	0.758	0.008	0.759	0.009
	T3 = 35	0.021	0.007	0.022	0.007
D3 = 3	T1 = 25	≫ 0.056	0.016	≫ 0.057	0.017
	T2 = 30	0.061	0.017	0.062	0.018
	T3 = 35	0.013	0.005	0.014	0.006

≫: LD50 is greater than 1 mL/cm3. a: EOCL alone. b: EOCL + (RC or GC) clay.

Table 2. LD₅₀ ​ of the Insecticidal Activities of EOCL Fixed on Porous Media (RC and GC).

Diameter (cm)	Temperature (°C)	LD50 (a) for RC	(b) for RC	LD50 (a) for GC	(b) for GC
D1 = 1	T1 = 25	0.228	0.110	0.229	0.115
	T2 = 30	0.350	0.109	0.349	0.111
	T3 = 35	0.035	0.017	0.036	0.018
D2 = 2	T1 = 25	0.529	0.032	0.530	0.034
	T2 = 30	0.758	0.008	0.759	0.009
	T3 = 35	0.021	0.007	0.022	0.007
D3 = 3	T1 = 25	≫ 0.056	0.016	≫ 0.057	0.017
	T2 = 30	0.061	0.017	0.062	0.018
	T3 = 35	0.013	0.005	0.014	0.006

≫: LD50 is greater than 1 mL/cm3. a: EOCL alone. b: EOCL + (RC or GC) clay.

2.4. The Bioinformatics Approach

Computational research was integrated as an additional step in this study to provide a comprehensive theoretical understanding of limonene, the primary molecule found in Citrus limonum, and its interaction with key proteins 2NXX and 1C3Z. These interactions offer unique insights into the biological mechanisms underlying development regulation and survival under extreme conditions, while also opening new avenues for protein docking applications [39]. Two distinct approaches were employed for this analysis. The first approach involved modeling the molecular interactions of limonene, with a particular focus on determining essential energetic descriptors. The second approach involved the docking study of limonene with selected proteins. For molecular modeling, the MMFF94 (Merck Molecular Force Field 94) method was employed, which is a set of force field parameters used in computational chemistry to simulate molecular interactions within the framework of molecular mechanics [40]. Several energetic descriptors were calculated, including the energy levels of the highest occupied (E_Homo) and lowest unoccupied (E_Lumo) molecular orbitals, as well as derived parameters such as the energy gap (EGAP), chemical hardness (η), electronegativity (χ), electrophilicity index (ω), and molecular flexibility (s). The formulas used to calculate these descriptors are listed in Table 3 [41].

In this study, the proteins encoded by the PDB structures 2NXX and 1C3Z were analyzed to understand their biological functions and interactions with limonene, a predominant molecule in Citrus limonum. These proteins exhibit distinctive structural features that influence their ligand-binding capabilities and present valuable opportunities for protein docking research [42].

The 2NXX protein, identified as the ultraspiracle receptor (USP) from the insect Tribolium castaneum, is part of the nuclear receptor family and is crucial in regulating developmental processes such as molting and metamorphosis. Remarkably, this receptor can maintain a functional conformation without requiring a specific ligand, setting it apart from its vertebrate counterparts. The protein surface, particularly the binding domain, plays a vital role in its molecular interactions. Analysis of the solvent-accessible surface (SAS) revealed the protein unusual ability to interact effectively without a ligand, highlighting its unique properties [43]. The 1C3Z protein, an antifreeze protein (THP12) from the yellow mealworm Tenebrio molitor, is vital for the survival of organisms in cold environments by preventing ice crystal formation and growth. THP12 helical structure and hydrophilic surface allow it to bind effectively to small hydrophobic ligands, enhancing its protective function against freezing. The stability of this protein under extreme conditions is significantly supported by electrostatic interactions on its surface [44].

The surface properties of 2NXX and 1C3Z, including their solvent accessibility and interaction potential with limonene, offer a robust foundation for protein docking studies. The ability of 2NXX to function without a ligand, combined with 1C3Z affinity for hydrophobic ligands such as limonene, facilitates advanced modeling of protein-ligand interactions. These models can be instrumental in designing targeted molecules and exploring the mechanisms of action of these proteins under different biological conditions [45]. Following the optimization of the protein structures using Swiss-PDB Viewer V4.1, several critical steps were undertaken for docking preparation. These included the removal of water molecules and heteroatoms, the addition of polar hydrogens, and the assignment of charges using the Kollman method. The geometry of limonene was optimized with density functional theory (DFT) calculations using the B3LYP/6-311G basis set. Docking simulations were conducted with Autodock VINA 1.1.2 (Table 4), using specific grid dimensions for each protein. The Biovia Discovery Studio (Studio, 2021) software was then utilized to visualize the protein-ligand interactions, allowing for a detailed analysis of how limonene interacts with the target proteins and providing valuable insights into their potential as biological agents [46].

2.5. Computer and Statistical Processing

All computational processing was performed using MATLAB software. The numerical modeling results obtained were compared with the experimental results to identify key parameters and maximize interpretation. In addition to MATLAB, other software tools were employed for various calculations and graphical representations, including EXCEL, ORIGIN LAB, DESIGN PLAN, and XLSTAT. These tools supported the comprehensive analysis and helped ensure that the results were accurately interpreted and effectively communicated. Furthermore, principal component analysis (PCA) has proven to be an indispensable mathematical tool for exploring the correlations between various calculated parameters. The goal was to visualize and synthesize data related to the energetic and qualitative properties of the molecular docking of limonene with the studied proteins [47,48].

3. Results and Analysis

3.1. Chemical Characterization of RC and GC Clay Media

X-ray Diffraction (XRD) Analysis: The XRD analysis revealed the mineralogical complexity of the RC and GC samples, which are primarily composed of phyllosilicates, especially kaolinite and illite. The distinct kaolinite peaks observed at 7.75 Å and 14.5 Å suggest that this mineral formed under conditions of low ionic charge and intense chemical weathering, which are typical of tropical or subtropical climates. Illite, found in both samples, points to diagenetic processes, likely linked to higher temperature and pressure conditions, indicating a complex geological history. The significant presence of quartz, with peaks at 19.8 Å and 26 Å, suggests a detrital input, probably of alluvial or aeolian origin. Specific reflections from carbonates like calcite and dolomite suggest formation in carbonate-rich environments, such as alkaline lakes or shallow marine zones. Finally, the detection of hematite in the RC sample, marked by reflections at 23.8 Å and 29.7 Å, indicates oxidation conditions, probably due to prolonged exposure to air, supporting the hypothesis of a pedogenic origin (Figure 2) [49,50].

X-ray Fluorescence (XRF) Chemical Composition Analysis: The XRF chemical analysis showed that silicon dioxide (SiO2) and aluminum oxide (Al2O3) are predominant in the samples, with notable variations between RC and GC. The GC sample, with its high SiO2 content, reflects a strong presence of silicate minerals such as quartz and silica-rich phyllosilicates, whereas the RC sample, with a more modest SiO2 content, is dominated by iron oxides, primarily Fe2O3, associated with hematite. The higher concentration of Al2O3 in GC suggests an abundance of aluminous clay minerals like kaolinite, while the lower content in RC may indicate a composition richer in non-aluminous minerals. Other oxides present, such as Na2O, K2O, MgO, CaO, and TiO2, are in moderate proportions and reflect the presence of accessory minerals like feldspars, dolomite, and calcite, indicating distinct sedimentary and diagenetic processes (Table 5) [51,52].

Transform Infrared Spectroscopy (FTIR) Analysis: FTIR analysis allowed us to identify the functional groups present in the samples, confirming the presence of clay minerals like kaolinite and illite, as well as carbonates and iron oxides. The absorption bands at 1025 cm⁻¹, associated with Si-O vibrations, reveal a well-ordered crystalline structure typical of aluminosilicates. The band at 920 cm⁻¹, linked to the Al2OH bond, further supports the identification of kaolinite, while the band at 3620 cm⁻¹ indicates the presence of structural water in the clay minerals. The bands at 475 cm⁻¹ and 540 cm⁻¹, representing Si-O and Si-O-Al bonds, confirm the presence of phyllosilicates. The band at 1517 cm⁻¹, corresponding to C-O vibration, indicates the presence of calcite, while the band at 1455 cm⁻¹ in RC, associated with hematite, explains the red color of this sample, revealing geochemical conditions favorable to iron oxidation (Figure 3) [53,54].

Scanning Electron Microscopy (SEM) and Energy Dispersive X-ray (EDX) Analysis: SEM analysis revealed significant morphological differences between the samples. The RC sample exhibits an irregular and microporous structure, suggesting intense weathering, potentially due to wetting and drying cycles or chemical interactions. This morphology is typical of clays subjected to pedogenetic processes. In contrast, the GC sample shows a more regular morphology, suggesting a more stable depositional environment. EDX analysis confirmed a high iron content in RC, mainly in the form of iron oxide (Fe2O3), corroborating the presence of hematite detected by XRD and FTIR. The GC sample is dominated by silicon and aluminum, reinforcing its richness in phyllosilicates like kaolinite and illite. These differences between RC and GC reflect distinct geological and environmental conditions (Table 6) (Figure 4 and Figure 5) [55,56].

Thermogravimetric Analysis (TGA): Thermogravimetric analysis identified three main phases of mass loss in the samples, reflecting dehydration and thermal decomposition processes. The first phase, below 200 °C, corresponds to the loss of surface water and light organic matter. The second phase, between 430 °C and 670 °C, corresponds to the dehydroxylation of clay minerals, particularly kaolinite, forming metakaolinite. This phase indicates a significant kaolinite content in both samples. The third phase, between 700 °C and 800 °C, corresponds to the decomposition of illites and carbonates like calcite and dolomite. The significant mass loss at this stage confirms the presence of these minerals, suggesting formation conditions that favor their stable precipitation. These thermogravimetric results reveal the complexity of the samples and their thermal behavior, providing insights into the geochemical conditions that influenced their formation (Figure 6) [57,58].

3.2. Chemical characterization of essential oil of Citrus limonum:

The analysis of Citrus limonum essential oil, conducted using gas chromatography coupled with mass spectrometry, reveals a detailed composition of the compounds present, each identified by its retention time (RT), retention index (RI), and area percentage (Area %), indicating their relative proportion in the sample. Limonene, which dominates with 51.41% of the total composition, is a major component typical of lemon essential oils. This monoterpene is not only responsible for the characteristic citrus aroma but is also recognized for its insecticidal properties, as well as its antioxidant, anti-inflammatory, and antimicrobial activities [59]. Additionally, β-Pinene (6.46%) and α-Pinene (4.84%), also common monoterpenes in citrus oils, contribute to the insecticidal efficacy and the aromatic complexity of the oil. γ-Terpinene (4.60%) and α-Terpineol (4.13%), present in notable amounts, also play a role in the biological activities, including antimicrobial and soothing effects [60]. The analysis also reveals several minor compounds, each representing less than 1% of the total, such as cis-Citral (1.27%), trans-Caryophyllene (1.06%), and cis-α-Bergamotene (1.97%). Although present in smaller quantities, these compounds can enhance the biological activity of the oil, particularly in the antimicrobial and insecticidal domains, while adding to the aromatic complexity [61]. Citrus limonum essential oil is predominantly composed of monoterpenes, notably limonene, β-Pinene, and α-Pinene, which are typical of citrus oils and responsible for its distinctive lemony aroma. Sesquiterpenes, such as β-Elemene and trans-Caryophyllene, though present in smaller amounts, can also play an important role in the overall aroma and therapeutic properties, particularly in enhancing insecticidal efficacy (Table 7) [62].

The analysis of Citrus limonum essential oil by gas chromatography coupled with mass spectrometry (GC-MS) shows significant variability in the chemical composition depending on the geographical origin of the sample. In several studies, limonene is consistently identified as the major compound, but its proportion varies considerably. Najwa Nasser AL-Jabri conducted a comparative analysis of lemon essential oils from Turkey and India. Limonene was the main compound in both varieties, representing 78.93% in the Turkish essential oil, followed by 5.08% β-pinene, while in the Indian variety, limonene represented 53.57% and α-terpineol reached 15.15% [59]. Similarly, Njoku et al. analyzed the essential oil of lemon from Nigeria, finding 54.15% limonene, followed by 6.21% 2-cyclohexan-1-ol and 3.82% β-pinen [60]. These variations indicate that the chemical composition of Citrus limonum essential oil is influenced by climate, soil, and other environmental factors specific to each region. In comparison, Bertuzzi et al. (2013) reported that lemon essential oil contains 71.87% limonene as the main constituent [61]. In another study, Himd et al. (2016) identified a limonene content of 64.19% in the essential oil of Lisbon lemon peel (Portugal), followed by β-pinene (7.76%) and α-terpinene (5.45%) [62]. However, a study by Moufida & Marzouk (2003) on lemon essential oil shows that geranial, at 3.65%, is the main compound, which is quite different from other studies [63]. These differences highlight the impact of local conditions and processing methods on the composition of Citrus limonum essential oils, emphasizing the importance of precise chemical characterization for each sample based on its geographical origin.

3.3. Insecticidal Activities in Porous Clay Media

The results of the insecticidal tests for EOCL, both alone and in combination with the porous clay media RC and GC, against Sitophilus granarius are summarized in Table 2, which presents the lethal doses of 50% LD₅₀ as a function of the studied parameters, such as cylinder diameter (D) and incubation temperature (T).

From these results, it is evident that the LD₅₀ ​ values decrease when EOCL is fixed on the porous clay media RC and GC, indicating an enhancement in insecticidal efficacy [64]. This suggests that the fixation of EOCL on the porous media alters its insecticidal potency. Additionally, the factors of cylinder diameter and temperature continue to play significant roles in maintaining the insecticidal properties, with larger cylinder diameters and higher temperatures generally promoting more effective insecticidal activity [65]. Furthermore, the analysis indicates that the insecticidal activity of EOCL when fixed on porous media, as well as its persistence, depends on the nature of the interactions between the essential oil compounds and the porous media. The specific interactions between the EOCL and the clay components likely contribute to the observed changes in insecticidal performance, underlining the importance of the media in modulating the efficacy of EOCL. This conclusion is consistent for both RC and GC clay types, suggesting a similar mechanism of action for the EOCL when used with these different porous media (Table 8) [66].

The experimental results reveal that the desorption of lemon essential oil in porous clays is significantly influenced by factors such as oil concentration, temperature, clay particle diameter, and the mass of clay used. These parameters determine the amount of essential oil available to act as an insecticide against Sitophilus granarius. Temperature plays a crucial role in this process, with desorption accelerating as the temperature increases, particularly between 25°C and 35°C [67]. For example, red clay exhibits faster mass loss at 35°C, indicating enhanced diffusion and potentially improved insecticidal efficacy at higher temperatures. Green clay demonstrates a slightly faster desorption rate under the same conditions, with desorption rates increasing from 30% to 60%, compared to 25% to 55% for red clay [68]. The concentration of the essential oil also affects desorption. Higher concentrations tend to slow initial desorption, likely due to increased clay saturation, which extends the insecticidal effect as more oil remains available over a longer period. For instance, at 25°C and with a particle diameter of 1 mm, the residual mass is higher for a stronger concentration, indicating slower but more sustained desorption. For red clay, the residual mass decreased from 75% to 45% as the temperature increased, while for green clay, it decreased from 70% to 40% [69].

Clay particle diameter influences desorption as well. Smaller particles, due to their higher specific surface area, promote increased interaction with the essential oil, leading to faster desorption and a more immediate insecticidal effect. The results show that green clay with finer particles has a higher desorption rate than red clay, highlighting the importance of particle size [70]. Lastly, a larger mass of clay leads to slower desorption of the essential oil, as more clay can adsorb and gradually release the oil, thus prolonging the insecticidal effect. This effect is more pronounced in red clay than in green clay, suggesting that red clay may be more effective for applications requiring prolonged release [71]. To maximize the insecticidal efficacy of lemon essential oil against Sitophilus granarius, it is crucial to consider the interaction between temperature, oil concentration, clay particle size, and clay mass. Faster desorption, facilitated by higher temperatures and finer particles, could provide a more immediate insecticidal effect, while slower desorption associated with higher concentration and greater clay mass could prolong the duration of the effect [72]. The analysis of desorption data in red and green clays allows for the identification of optimal conditions to maximize insecticidal activity, which is essential for optimizing the storage and application conditions of lemon essential oil. Regarding the estimation of LD₅₀, it was determined that reliable values require well-distributed data covering an adequate mortality range. Through adjustments to the logistic model parameters and linear interpolation, the LD₅₀ for lemon essential oil against Sitophilus granarius was estimated to be 0.02 µl/cm³, highlighting the importance of proper data distribution for accurate analysis (Table 9) (Figure 7) [73].

3.4. Optimization of Essential Oil Transfer Conditions in Porous Media

The optimization of the evaporation flux parameter for essential oil transfer in porous media was conducted using a full factorial design involving four factors: concentration of essential oil (C), temperature (T), cylinder diameter (D), and mass of clay medium (M). Each factor was tested at two levels, coded as +1 for the high level and -1 for the low level, resulting in a total of 16 simulations. The mathematical model associated with this first-degree polynomial accounts for the main effects of these factors, as well as their second, third, and fourth-order interactions [74]. The coefficients for the flux model of essential oil in the red clay (RC) medium indicate that the mean flux achieved was 2.6061×10^-4 Similarly, the coefficients for the flux model of essential oil in the green clay (GC) medium provide insight into the behavior of the oil in this different porous medium. In both RC and GC, an increase in the concentration of essential oil enhances the flux, while a rise in temperature also contributes to an increased flow [75]. Conversely, increasing the cylinder diameter and the mass of the clay medium results in a reduction of the flux in both types of clay. Graphical representations of the flux as a function of two interacting factors in 3D for both RC and GC demonstrate that the interaction between concentration and temperature significantly increases the flux. However, interactions between concentration and diameter, concentration and mass, mass and temperature, and temperature and diameter generally result in a decrease in flux. The interaction between diameter and mass also leads to a reduction in the flow rate for both red and green clays [76]. Overall, the determination coefficients R² for the models developed for both RC and GC media indicate that the optimized models are well-explained by the data. The Fisher test value (F_test = 2.42) suggests that the model is not significantly influenced by noise or uncertainties. However, the predictive coefficient of determination (R²_pred ​) yielded a negative value, implying that the total mean flux might be a better predictor of the model response [77]. Despite this, the response adequacy, measured by the flux-to-uncertainty ratio, was found to be 6.68, exceeding the desirable threshold of 4, which confirms the adequacy of the response. This analysis underscores the complex interactions between various factors affecting the flux of essential oil in porous media, particularly in red and green clays, and highlights the importance of optimizing these conditions to ensure efficient transfer of essential oil (Figure 8 and Figure 9) [78].

3.5. Modeling the Diffusion Process of Citrus Limonum Essential Oil (EOCL) in Porous Media RC and GC

To better understand how Citrus Limonum essential oil (EOCL) diffuses in porous media like RC and GC, two complementary approaches were taken: studying the kinetics of desorption and simulating concentration profiles within a cylindrical model. These steps are crucial for identifying key parameters such as the diffusion coefficient Dz, evaporation flux F, evaporation constant K, and activation energy Ea, all of which are essential for accurately describing the oil behavior in these media. The Fickian diffusion model was chosen to analyze desorption, based on the mass percentage of oil over time [79]. For the RC medium, the diffusion coefficient Dz was estimated at 9.81×10^-6 cm²/h. under optimal conditions (C = 0.01 mL/cm³, T = 30°C, M = 0.10 g, D = 2 cm). The results showed a significant decrease in oil mass, confirming the validity of the model used [80]. This simulation is not only precise but also useful for predicting the remaining amount of oil over time, which is crucial for applications such as insecticidal efficacy. Studies conducted on the GC medium produced similar results, though some adjustments are necessary due to the specific properties of this medium [81]. Overall, the Fickian diffusion model has proven to be an effective tool for understanding and predicting EOCL diffusion, providing a solid foundation for optimizing the use of essential oils in various practical applications (Figure 10 and Figure 11), [82].

3.6. Interaction of Essential Oils in Porous Media and Insecticidal Efficacy

The study of the behavior of Citrus limonum essential oil (EOCL) in porous media, particularly in red clay (RC) and green clay (GC), is based on several key parameters: evaporation flux (F), diffusivity (Dz), evaporation constant (K), activation energy (Ea), and insecticidal activity, measured by mortality percentage (M%). These parameters, detailed in Table 8 and Table 9, provide a complex but valuable insight into how EOCL interacts with these porous media. The activation energy values, while consistent and low across both media, suggest similar energetic behavior between RC and GC. The observed negative values indicate a dominant adsorption process, highlighting the strong capacity of these clays to retain essential oils. This retention capability is crucial for understanding the desorption dynamics and the insecticidal efficacy of EOCL [83]. To further analyze these interactions, Principal Component Analysis (PCA) was conducted, focusing on evaporation flux (F), diffusivity (Dz), evaporation constant (K), and mortality percentage (M%) across 16 different cases. The PCA identified two principal components (F1 and F2), which together explained 82% of the total variability, with eigenvalues of 1.89 and 1.38, respectively [84]. The mapping of these components revealed that the diffusion coefficient (Dz) and the evaporation constant (K) are closely related, as indicated by their proximity on the F1 axis. This correlation suggests a dependency between these two parameters, aligning with theoretical expectations. On the F2 axis, the evaporation flux (F) and mortality percentage (M%) also show a strong correlation, implying that insecticidal efficacy, as measured by mortality, is largely influenced by the rate of EOCL evaporation. Interestingly, the near-orthogonal relationship between the pairs (Dz, K) and (F, M%) suggests independence between diffusion-related parameters and evaporation flux-related parameters [85]. This distinction underscores the complexity of the interactions between the physical diffusion processes within the clay matrices and the resulting biological insecticidal activity. The PCA not only clarifies the interactions between EOCL and the clay media but also highlights the key factors influencing its insecticidal effectiveness. Understanding these relationships is essential for optimizing the use of essential oils in pest control applications, where the choice of porous medium can significantly impact the efficacy of the treatment (Table 10, Table 11 and Table 12), (Figure 12) [86].

The results provide a detailed analysis of the correlation between several key parameters involved in the desorption of EOCL essential oil and its insecticidal activity in two types of clay media: red clay (RC) and green clay (GC). The parameters examined include evaporation flux (F), diffusivity (Dz), evaporation constant (K), and mortality percentage (M%). In the RC medium, the F1 axis, representing 47.31% of the total variability, is the most significant in explaining the desorption behavior of EOCL. A strong correlation between diffusivity (Dz) and the evaporation constant (K) is observed, suggesting that changes in diffusivity directly influence the evaporation constant. This highlights the crucial role of diffusion in the desorption dynamics within this medium [87]. Simultaneously, the F2 axis, accounting for 34.70% of the variability, reveals a close relationship between the evaporation flux (F) and the mortality percentage (M%). This indicates that the evaporation rate of EOCL has a direct impact on its insecticidal efficacy, with faster evaporation being associated with higher mortality rates. The orthogonal relationship between the pairs (Dz, K) and (F, M%) suggests that diffusion processes are independent of evaporation flux and the resulting insecticidal activity, revealing a complex interplay between physical diffusion processes and biological effectiveness [88]. In the GC medium, the F1 axis similarly accounts for 47.31% of the total variability and shows a strong correlation between diffusivity (Dz) and the evaporation constant (K). This underscores the critical role of diffusion in shaping the evaporation characteristics of EOCL in this medium. Furthermore, the F2 axis, which explains 34.70% of the variability, shows a significant correlation between the evaporation flux (F) and the mortality percentage (M%), echoing the findings in the RC medium. This highlights that the evaporation rate of EOCL is a key factor in its insecticidal efficacy. The orthogonal relationship between the pairs (Dz, K) and (F, M%) in the GC medium, as in the RC medium, underscores the independence of diffusion processes from the factors controlling evaporation rate and insecticidal activity, emphasizing the complexity of these interactions in green clay [89]. Overall, in both RC and GC media, the strong correlation between diffusivity (Dz) and the evaporation constant (K) underscores the critical importance of diffusion in the EOCL desorption process. The evaporation flux (F) emerges as a key determinant of insecticidal activity (M%), with faster evaporation improving insecticidal efficacy. The independence between diffusion-related and evaporation-related parameters indicates that these factors operate through distinct mechanisms. Therefore, it is essential to consider them separately when optimizing EOCL desorption and insecticidal effectiveness [90].

3.7. Bioinformatic Studies

The bioinformatics study relies on two distinct computational approaches focused on analyzing the limonene molecule and its interaction with two key proteins, 2NXX and 1C3Z. The first approach involves the energetic analysis of limonene, while the second explores its molecular docking with these proteins, which are crucial for understanding biological processes such as development and survival under extreme conditions. The results of the energetic analysis are presented in (Figure 13) and (Table 13). The Figure 13 displays the energy level diagram of the molecular orbitals of limonene, optimized using the MMFF94 method, while Table 13 summarizes additional calculations and relevant energetic parameters. These data reveal that limonene exhibits significant electronic stability, characterized by an energy gap (EGAP) of 6.677 eV and a chemical hardness of 3.338 eV, indicating a strong resistance to reactivity [91]. An electronegativity of 7.634 eV suggests a robust electron-attracting ability, contributing to its stability. The electrophilicity index of 8.729 eV indicates a moderate capacity of the molecule to accept electrons during chemical reactions, while a molecular flexibility value of 0.149 eV points to moderate reactivity. In summary, limonene demonstrates high electronic stability, low reactivity, and a moderate ability to act as an electrophile in chemical reactions. These results are crucial for understanding the electronic behavior and potential reactivity of limonene, providing a solid foundation for further investigation of its interactions with the 2NXX and 1C3Z proteins. These interactions are particularly relevant given the biological roles of these proteins, and the findings from the molecular docking studies will shed light on the potential applications of limonene in biological and therapeutic contexts [92].

The docking results of limonene with the proteins 2NXX and 1C3Z reveal subtle yet significant differences in the interaction between the ligand and these two proteins. The binding free energy of -6.5 kcal/mol for 2NXX indicates a slightly more stable interaction with limonene compared to the binding free energy of -6.3 kcal/mol observed with 1C3Z. Although this difference is modest, it suggests that the binding site of 2NXX is better suited to accommodate limonene, providing a more favorable environment for interaction [93]. Furthermore, the pKi of 4.77 for 2NXX, slightly higher than the 4.62 for 1C3Z, confirms a stronger binding affinity of limonene for 2NXX. This higher pKi reflects greater stability of the formed complex, reinforcing the idea that 2NXX is a more suitable partner for limonene [94]. The ligand efficiency, measured at 0.65 kcal/mol/non-H atom for 2NXX compared to 0.63 for 1C3Z, also indicates more efficient binding of limonene to 2NXX. This result could be attributed to better steric compatibility or a more favorable alignment of binding surfaces in 2NXX. The torsional energy, which is identical for both complexes at 0.3113, suggests that the conformational flexibility of limonene remains constant regardless of which protein it binds to. This indicates that the energetic cost for limonene to adapt to the binding site conformation is similar in both cases, highlighting good conformational compatibility with both proteins [95]. In summary, while both 2NXX and 1C3Z offer binding sites that are well-suited to limonene, the data indicate a slight preference for 2NXX. This preference is supported by a more negative binding free energy, a higher pKi, and a slightly superior ligand efficiency. These results suggest that 2NXX may provide a more optimal binding site for limonene, which could have implications for its biological function in the presence of similar compounds. However, the close similarity of the results for both proteins indicates that limonene is an effective ligand for each, with a slightly stronger affinity for 2NXX (Figure 14), (Table 14) [96].

The principal component analysis (PCA) presented in the graph highlights the complex relationships between docking parameters and the studied proteins. The F1 axis, as the first principal component, accounts for a significant portion of the variance in the data. Parameters that deviate most from the origin along this axis contribute the most to the overall variability of the results. In contrast, the F2 axis, the second principal component, explains a smaller portion of the variance, as indicated by the annotation "F2 (0.00%)." Parameters near this axis show little variation compared to F1. In this analysis, Binding Free Energy and Ligand Efficiency stand out with thicker arrows, emphasizing their importance. Their proximity suggests a positive correlation: an increase in binding affinity is generally associated with better ligand efficiency [97]. On the other hand, pKi, represented by an arrow pointing in a different direction, appears to capture a distinct aspect of ligand-protein interactions, likely related to the strength of inhibition or the ligand affinity for the protein. The vertical orientation of Binding Free Energy and Ligand Efficiency indicates that these parameters are primarily influenced by F1, confirming the minimal contribution of F2 [98]. The black dots, representing the proteins, are positioned based on their coordinates in the principal component space. Their proximity to certain arrows reveals which proteins are most influenced by these parameters. A protein located near Binding Free Energy or Ligand Efficiency might indicate that the interaction between the ligand and this protein is strongly determined by these parameters, whereas a protein far from these arrows could be influenced by other factors. These findings suggest that certain proteins exhibit varying binding affinities or efficiencies depending on the characteristics of the ligand. This understanding is crucial for drug development and for deepening knowledge of the biological mechanisms underlying molecular interactions [99]. The correlation circle provides a visual representation of these correlations: arrows reaching the edge of the circle are strongly correlated with the F1 or F2 axes, while those near the center show weaker correlations. The balanced distribution of parameters around the circle indicates that they capture diverse dimensions of ligand-protein interactions, thereby enriching the overall analysis. This PCA offers a clear and in-depth visualization of the complex relationships between multiple docking parameters and the studied proteins, providing valuable insights into the molecular mechanisms underlying these interactions (Figure 15) [100].

4. Conclusions

The study highlights the promising effectiveness of encapsulating lemon essential oil (Citrus limonum) within natural porous clay matrices to combat the pest Sitophilus granarius, a significant concern for the preservation of stored grains. The research demonstrated that using porous clays, sourced from the Bejaad region in Morocco, as carriers for the essential oil helps to overcome the limitations related to the volatility and rapid degradation of these active compounds while optimizing their insecticidal potential. Analyses revealed that the fixation of the essential oil on clay matrices significantly enhances its insecticidal efficacy, reducing the lethal dose required to achieve the desired effect and allowing for a more controlled and prolonged release of the active compound. This increased effectiveness is attributed to specific interactions between the essential oil and the mineralogical components of the clays, which favorably influence the diffusion and persistence of limonene, the main bioactive compound in the essential oil. Furthermore, mathematical modeling of the diffusion process, combined with rigorous statistical approaches, allowed for a detailed understanding of the mechanisms of essential oil release from the clay matrices. These analyses showed that parameters such as temperature, oil concentration, clay particle size, and mass play a crucial role in the desorption rate and, consequently, in the insecticidal efficacy. A fine understanding of these mechanisms enables the adjustment of application conditions to maximize effectiveness while minimizing risks to the environment and human health. The study also opened up new perspectives through bioinformatic analysis of molecular interactions between limonene and certain target proteins in insects, suggesting that these interactions could enhance the effectiveness of the essential oil as a natural biopesticide. This paves the way for innovative and sustainable strategies for the development of natural biopesticides, offering viable alternatives to traditional chemical pesticides. In conclusion, this work makes a significant contribution to the understanding and optimization of essential oil use in biological control applications. It paves the way for innovative strategies for the development of natural biopesticides, offering sustainable alternatives to conventional chemical pesticides while maintaining the necessary effectiveness for crop protection. This integrated approach, combining material characterization, theoretical modeling, and bioinformatic analyses, could serve as a model for future research in the field of biological control using essential oils.