Synthesis, Insecticidal Activity and Computational Studies of Eugenol-Based Insecticides †

: Eugenol, a natural phenolic allyl benzene that has been used as an active lead compound showing signiﬁcant biological activities, including insecticidal effects on a wide variety of domestic arthropod pests, was used as the main reagent in the present work. Ester eugenol derivatives were synthesized and evaluated for their insecticidal activities against the Spodoptera frugiperda cell line. Studies of structure-based inverted virtual screening were carried out in order to identify the potential targets associated with the obtained insecticidal activity. The results indicate that the insecticide activity observed is most likely a result of the interaction of these molecules with the odorant-binding proteins and/or with acetylcholinesterase.


Introduction
The global population is increasing at an exponential rate, so it is necessary to ensure agricultural production that meets the actual food requirements. In crop protection, the reduction in damage caused by pathogens and pests in agricultural fields is mainly achieved through the extensive use of synthetic pesticides. To mitigate the environmental problems caused by the intensive use of conventional synthetic pesticides, biopesticides and semisynthetic pesticides based on natural plant products are alternatives as pest-management agents [1,2].
Considering these facts, in the present work, eugenol derivatives were synthesized through an esterification reaction and evaluated for their effect on the viability of Sf9 cells. A structure-based inverted virtual screening protocol was employed to identify the potential proteins associated with the observed insecticidal activity.

Biological Activity of Compounds 3a-c
The impact of ester eugenol derivatives 3a-c on the viability of Sf9 cells was evaluated at 100 μg/mL, following 24 h of exposure. As shown in Figure 1, compounds 3a and 3b, resulting from the eugenol esterification with chlorobenzoic acids, displayed similar toxicity to the starting material eugenol 1. On the other hand, the eugenol esterification

Biological Activity of Compounds 3a-c
The impact of ester eugenol derivatives 3a-c on the viability of Sf9 cells was evaluated at 100 µg/mL, following 24 h of exposure. As shown in Figure 1, compounds 3a and 3b, resulting from the eugenol esterification with chlorobenzoic acids, displayed similar toxicity to the starting material eugenol 1. On the other hand, the eugenol esterification with 3-amino-2-methylbenzoic acid led to increased activity (compound 3c), causing nearly 50% of viability loss. with 3-amino-2-methylbenzoic acid led to increased activity (compound 3c), causing nearly 50% of viability loss.

Inverted Virtual Screening Results
In Table 1, the average score for all the eugenol derivatives is presented for each potential target using each scoring function. In each set of targets, the structure with the highest score was selected and ranked from best to worst, based on the docking programs/scoring functions' predictions. For this study, the four scoring functions (SFs) of GOLD were used. The score of each of the GOLD SFs is dimensionless, with a higher value indicating a better binding affinity.

Inverted Virtual Screening Results
In Table 1, the average score for all the eugenol derivatives is presented for each potential target using each scoring function. In each set of targets, the structure with the highest score was selected and ranked from best to worst, based on the docking programs/scoring functions' predictions. For this study, the four scoring functions (SFs) of GOLD were used. The score of each of the GOLD SFs is dimensionless, with a higher value indicating a better binding affinity. There is a high degree of consistency across all SFs, with odorant-binding proteins (OBPs) and acetylcholinesterases (AChEs) suggesting more likely binding. The octopamine Chem. Proc. 2022, 12, 46 4 of 9 receptor, voltage-gated sodium channel, polyphenol oxidate and N-acetylglucosamine-1phosphate uridyltransferase (GlmU) consistently exhibit lower scores.

Molecular Dynamics Simulations and Free Energy Calculation Results
Simulations were carried out for both groups of targets predicted at the inverted VS stage odorant binding proteins and acetylcholinesterases, in complex with the three most potent eugenol derivatives. The structures chosen were the ones that presented the best score from each group (3K1E for OBP and 1QON for acetylcholinesterases-AChE). Molecular dynamics simulations were used to confirm and validate the inverted screening predictions, allowing for a more detailed analysis. Additionally, the interactions formed between the protein and ligand were further analyzed, and the most determinant residues were identified. The results are presented in Table 2.
The protein RMSD value of OBP compared to the docking pose had an average value of about 2.3 Å. Interestingly, the AChE complexes presented a higher RMSD value. However, in these cases, the standard deviation was low. It is possible that the complexes of AChE and the eugenol derivatives were optimized at the beginning of the simulation to obtain a more stable conformation. Throughout the simulation, all molecules remained bound to their targets and an induced-fit adjustment was observed. A study was also conducted to determine how buried the eugenol derivatives were in the binding pockets, measuring the solvent accessible surface area (SASA) and the percentage of potential SASA. An increase in the percentage of ligand SASA with a lower SASA indicates that the molecule is buried in the target pocket, making it less exposed to solvent. Compound 3b is the most buried in the OBP pocket (with a percentage of ligand SASA of 97% and an average SASA of 18.5 Å 2 ). Regarding AChE, compound 3a is the most buried in the active site (with a percentage of ligand SASA of 94% and a SASA of 31.7 Å 2 ).
The three eugenol derivatives exhibit slightly better binding affinities toward OBP, with compound 3a showing a ∆G bind of −35.6 kcal/mol when in complex with OBP and −33.5 kcal/mol when bound to AChE. When compared to the other compounds, compound 3c is the weakest binder for both OBP and AChE (with ∆G bind of −30.1 kcal/mol when bound to OBP and −29.1 kcal/mol when bound to AChE).
When bound to OBP, the compounds are stabilized primarily by electrostatic interactions with Trp114, Leu76, Gly92, Phe15, Leu80 and Tyr122. From all the compounds studied, the results seem to suggest that compounds 3a and 3b may be good candidates to be used as repellents, having OBP as their main target. Regarding AChE, the main interacting residues are Tyr71, Trp83 and Tyr374.

Biological Assays
The potential of compounds 3a-c as biopesticides was evaluated in assays using the Sf9 (Spodoptera frugiperda) insect cell line. Cells were maintained at 28 • C and cultivated in Grace's medium with 10% FBS. For the evaluation of viability, cells were plated at 3.0 × 10 4 cells/well and exposed to the molecules, after which resazurin was added, and the fluorescence was read at 560/590 nm after 60 min of incubation.

Docking and Inverted Virtual Screening Studies
The Scopus database was searched for papers reporting virtual screening studies involving targets and molecules with insecticidal activity. The publication year and relevance of the target were considered in the selection process. A total of thirteen targets were identified in the eighteen studies found, as shown in Table 3. After preparing each PDB structure (removing water and other molecules from the crystallization process), the crystallographic ligands were extracted and saved in separate files, to provide a reference site for docking coordinates and to posteriorly perform redocking. In the absence of crystallographic ligands, active site residues were selected based on their importance for activity. In order to evaluate the quality of the docking protocol, as well as the capability of the docking software to reproduce geometry and orientation of the crystallographic pose, re-docking was used. The docking program/scoring functions used was GOLD [29] (PLP, ASP, ChemScore and GoldScore). By comparing the predicted docking pose with the crystallographic one through RMSD calculation, the ability of the SF to correctly dock the ligand is evaluated and that is a valuable step in the validation process. A lower RMSD means a better docking prediction.
The three eugenol derivatives were prepared for the study using DataWarrior [30] and OpenBabel [31], and subsequently docked into each PDB structure, using all the four SFs as soon as the protocol was optimized. The parameters that were optimized for each program/scoring function were: the docking coordinates, the docking box dimension or radius, the exhaustiveness, the search efficiency and the number of runs. Lastly, a ranked list was prepared based on the average scores of each target.

Molecular Dynamics Simulations and Free Energy Calculations
The three eugenol derivatives (compounds 3a, 3b and 3c) in complex with the two most promising targets identified from the inverted virtual screening study (odorantbinding protein 1-3KIE and acetylcholinesterase-1QON) were further studied through 100 ns molecular dynamics simulations with the Amber18 software [32].
The pose predicted at the IVS stage (with the PLP scoring function) was selected as the starting point of each MD simulation. ANTECHAMBER with RESP HF/6-31G(d) charges (calculated with Gaussian16 [33] and the general Amber force field (GAFF) [34]) were used to assign parameters. To describe the protein targets, the ff14SB force field [35] was applied. Posteriorly, the protein-ligand complexes were placed in TIP3P water boxes with a minimum distance of 12 Å between the protein surface and the side of the box, and periodic boundary conditions were used. By adding counter-ions (Na + ), the overall charge on the system was neutralized. Ewald's particle-mesh summation method was used to calculate long-range electrostatic interactions, while for short-range electrostatic and Lennard-Jones interactions, a cut-off value of 10.0 Å was used. A time step of 2 fs was used and all bonds involving hydrogen atoms were constrained using the SHAKE algorithm.
Then, the systems were submitted to four consecutive minimization stages to remove clashes, followed by an equilibration and production run, with each minimization having a maximum of 2500 cycles. In the equilibration run, the procedure was divided into two phases; NVT ensemble, where the systems were gradually heated to 298 K using a Langevin thermostat at constant volume (50 ps), and equilibration of the system density at 298 K (subsequent 50 ps). Lastly, the production run was performed during 100 ns with an NPT ensemble at constant temperature (298 K, Langevin thermostat) and pressure (1 bar, Berendsen barostat). The last 70 ns of the simulation were considered for SASA and hydrogen bonding analysis.
In this study, the molecular mechanics/generalized Born surface area method [36] was used in conjunction with the MM/PBSA.py script [37] from Amber. Each simulation's last 70 ns were analyzed with an interval of 100 ps and salt concentration of 0.100 mol dm −3 . Additionally, the contribution of the amino acid residues was evaluated using the energy decomposition method. The MM-GBSA calculations considered 1400 conformations from each MD trajectory taken from the last 70 ns of simulation.

Conclusions
The eugenol derivatives obtained were fully characterized and their biological evaluation as insecticides using the Sf9 (Spodoptera frugiperda) insect cell line has shown that it was possible to obtain a molecule (compound 3c) that was more toxic by tuning its structure.