Exploring the Potentiality of Natural Products from Essential Oils as Inhibitors of Odorant-Binding Proteins: A Structure- and Ligand-Based Virtual Screening Approach To Find Novel Mosquito Repellents.

Odorant-binding proteins (OBPs) are the main olfactory proteins of mosquitoes, and their structures have been widely explored to develop new repellents. In the present study, we combined ligand- and structure-based virtual screening approaches using as a starting point 1633 compounds from 71 botanical families obtained from the Essential Oil Database (EssOilDB). Using as reference the crystallographic structure of N,N-diethyl-meta-toluamide interacting with the OBP1 homodimer of Anopheles gambiae (AgamOBP1), we performed a structural and pharmacophoric similarity search to select potential natural products from the library. Thymol acetate, 4-(4-methyl phenyl)-pentanal, thymyl isovalerate, and p-cymen-8-yl demonstrated a favorable chemical correlation with DEET and also had high-affinity interactions with the OBP binding pocket that molecular dynamics simulations showed to be stable. To the best of our knowledge, this is the first study to evaluate on a large scale the potentiality of NPs from essential oils as inhibitors of the mosquito OBP1 using in silico approaches. Our results could facilitate the design of novel repellents with improved selectivity and affinity to the protein binding pocket and can shed light on the mechanism of action of these compounds against insect olfactory recognition.

Introduction

Mosquitoes are the main agents of vector-borne diseases in the tropical regions of the world, causing a high social, economic, and public health impact on these affected regions.[1,2] The disruption of mosquito–human interaction remains one of the most efficient prophylaxis methods against these diseases, and research on chemical repellents against mosquitoes has advanced due to the understanding of mosquito behavior and chemical olfactory receptors.[3−5]

The olfactory system of insects involves diverse transmembrane odorant receptor proteins located in olfactory membrane neurons, which are expressed in different parts of the insect body.[6] These odorant receptors evolved to respond to several functions in the mosquito life cycle, such as identification of pheromones for reproduction and chemical signals for host recognition.[6−8] The odorant-binding protein 1 (OBP1) is the main olfactory protein involved in the host-seeking mechanism of mosquitoes. The 3D structure of OBP1 is well-conserved across different mosquito species that are vectors of human diseases, such as Aedes aegypti (Protein Data Bank ID: 3K1E),[9]Anopheles gambiae (3V2L and 3R1O),[10] and Culex quinquefasciatus (3OGN),[11] and its structure has been widely investigated for structure-oriented development of novel repellents.[12−14] DEET (N,N-diethyl-meta-toluamide) is one of the most effective, commercially available Food and Drug Administration (FDA)-approved repellents, providing good residual protection against a broad spectrum of insects.[15,16] Its action against OBPs has been extensively investigated as an attractive approach to developing novel bio-inspired repellent compounds by combining ligand- and structure-based approaches.[17,18] However, different studies have reported behavioral insensitivity of different insect species to DEET, which could implicate some inefficiency in its repellency activity.[19−21]

With the increased interest in developing new mosquito repellents from natural products (NPs), essential oils are considered an interesting source due to their widely diverse class of volatile and low-molecular-weight compounds and also to their ovicidal, larvicidal, and repellent activities against human disease vectors.[22−26] With the resurgence of NPs in the development of new bioactive compounds by the pharmaceutical and cosmetic industries[27−30] and due to the cutting-edge technologies of combinatory chemistry, cheminformatics, and molecular modeling, new studies have focused on essential oils to explore their potentiality as mosquito repellents.[23,31,32] Recently, we have used different computational approaches to investigate biomolecular systems with emphasis on enzymatic reaction and inhibition.[33−37] In the present study, using an in silico approach, we performed a comprehensive analysis of the potentiality of 1633 compounds from the essential oils of 71 botanical families deposited in the Essential Oil Database (EssOilDB)[38] by combining a structure- and ligand-based virtual screening, using as reference the structure of DEET complexed to the OBP1 homodimer of A. gambiae (AgamOBP1). We also investigated the affinity and selectivity of these compounds against AgamOBP1 through docking techniques allied with molecular dynamics (MD) simulation and binding free energy calculations.

Results and Discussion

In the present study, we applied a structure- and ligand-based virtual screening approach starting with 1633 compounds from 71 botanical families obtained from the EssOilDB (see the complete list in Table S1). We used as reference the crystallographic structure of DEET interacting with the AgamOBP1 homodimer, a relevant odorant protein involved in host-seeking recognition by insects. Using a pharmacophoric prediction and structural similarity search, we analyzed the similarity of these compounds with the reference structure. Then, using docking techniques, MD simulation, and binding free energy calculations, we investigated the potentiality of these compounds as inhibitors of AgamOBP1. Our results are further discussed in the following sections.

Structure-Based Filtering and Similarity Search

A total of 121 compounds from essential oils were obtained after 3D structural similarity filtering using a Tanimoto cutoff of 0.7. After pharmacophoric filtering, 29 compounds were obtained. Two pharmacophoric models were selected: (1) The first model considered the interatomic interactions of the crystallographic structure together with the findings of a previous 3D quantitative structure–activity study of DEET.[4,39] This model comprises one aromatic hydrophobic function located in the aromatic ring, two aliphatic hydrophobic functions at the methyl groups of the diethylamine chain, and one hydrogen-bond acceptor formed by the carboxylic group (Figure A). (2) The second model considers the relevance of apolar interactions of DEET in the hydrophobic tunnel and contains an additional aliphatic hydrophobic function formed by the methyl group adjacent to the aromatic ring (Figure B). The molecular structure of DEET is shown in Figure C.

Figure 1

(A, B) Selected pharmacophoric models used to filter the NPs from essential oils. The green circles highlight the aromatic and aliphatic hydrophobic functions, and the orange ones are for the hydrogen-bond acceptor function. (C) Molecular conformer of DEET complexed to AgamOBP1.

Based on the physicochemical and structural correlation of the NPs with commercial repellents obtained with PCA, we selected six compounds to determine their selectivity and affinity to the AgamOBP1 binding pocket using molecular docking, MD simulation, and binding free energy calculations. The carvacryl acetate and thymol acetate were selected due to their well-known repellent activity against different insect species;[40−43] thus, they could be used for a comparative analysis with the other natural products. The compounds thymyl isovalerate, 4-(4-methylphenyl)-pentanal, and p-cymen-8-yl were selected due to their closely physicochemical correlation with commercial repellents DEET and DEPA, as shown by the PCA plot, and the p-anisyl hexanoate was selected due to its structural and pharmacophoric similarities to crystallographic DEET. The common name, chemical class, structural similarity (RMSD and Tanimoto 3D), and some origin species of the selected compounds are shown in Table .

Table 1

Selected NPs from Essential Oils To Investigate Their Affinity against AgamOBP1 Binding Pocketa

The compounds are identified by molecular structure, common name, chemical class, structural similarity to DEET (Tanimoto 3D value and RMSD), and some origin species.

The compounds p-cymen-8-yl, thymol acetate, carvacryl acetate, 4-(4-methylphenyl)-pentanal, thymyl isovalerate, and p-anisyl hexanoate showed satisfactory structural similarity with DEET according to the following Tanimoto 3D values: 0.80, 0.79, 0.71, 0.80, 0.78, and 0.81, respectively. Supporting this structural similarity, p-cymen-8-yl, thymol acetate, carvacryl acetate, 4-(4-methylphenyl)-pentanal, and thymyl isovalerate also exhibited a favorable physicochemical correlation with commercial repellents. Figure depicts the PCA scatter plot; the first circle highlights the correspondence between the NPs 2-methoxy-4,5-methylenedioxypropiophenone, methyl N-methylanthranilate, and (S)-1-(4-acetoxyphenyl)propyl acetate with the repellents dimethyl phthalate (DBP), methylanthranilate, and 3-cyclohexyl propanoic acid. The second circle highlights the chemical correspondence between NPs thymol acetate, carvacryl acetate, 4-(4-methyl phenyl)-pentanal, benzyl (2S)-2-methylbutanoate, thymyl isovalerate, and p-cymen-8-yl with the commercial repellents N,N-diethyl phenylacetamide (DEPA) and DEET. We obtained the following variance percentage for the principal components: 44.22% (PC1), 26.10% (PC2), and 16.79% (PC3). The raw data of compound properties (NPs and repellents) used to calculate the PCA is available in Table S2.

Figure 2

Scatter plot showing the chemical space of the NPs and the commercially used repellents. The x, y, and z axes exhibit the contribution of each principal component for the chemical profile of the compounds. The selected compounds are identified by color arrows: p-cymen-8-yl (blue), thymol acetate (red), carvacryl acetate (yellow), and 4-(4-methylphenyl)-pentanal (green), thymyl isovalerate (black), and p-anisyl hexanoate (orange).

Interestingly, several studies have highlighted the repellent activity of thymol and thymol acetate against different insects, which are vectors of diseases such as Ixodes ricinus (Acari: Ixodidae),[40]Anopheles subpictus Grassi, Anopheles stephensi (Diptera: Culicidae),[42,44]Culex pipiens (Diptera: Culicidae),[41] and crop pests such as Meligethes aeneus (Fabricius) (Coleoptera: Nitidulidae).[31] Similarly, carvacryl acetate, a monoterpene derivative of carvacrol found in a high percentage in the essential oils of some Lamiaceae species such as Clinopodium sp.[45] and Thymus sp.,[46] has a well-reported insecticidal and oviposition deterrence against different insect species.[54,55] In addition, p-cymen-8-ol has a high repellent activity against Lasioderma serricorne Fabricius (Coleoptera: Anobiidae),[47,48] and its derivative molecule p-cymene was used against A. gambiae.(49)

Several studies have reported the specific inhibitory activity of some NPs from essential oil against insect OBPs with satisfactory binding affinity, such as sesquiterpenes (e.g., benzaldehyde, β-myrcene, and α-copaene)[13] and alcohols (e.g., n-decanol and n-dodecanol).[12] In the present study, we identified that predominantly monoterpenoids, such as p-cymen-8-yl, thymol acetate, and carvacryl acetate, mimic the binding mode of DEET, exhibiting similar pharmacophoric groups and intermolecular interactions with the protein pocket, which sheds light on the molecular mechanism of action of these compounds against olfactory recognition of insects, reinforcing the previous experimental repellency studies for this chemical class.[40,41,48,50]

Interactions of NPs with the AgamOBP1 Binding Pocket

The AgamOBP1 structure is composed of two monomers, each consisting of six α-helices, with the odorant-binding pocket located at the center of a hydrophobic tunnel through the dimeric interface.[39,51] The crystallographic structure of AgamOBP1 bound to DEET reveals that the binding pocket cavity is formed by the residues Leu80, Leu73, Leu76, and His77 (α-helix 4); Met91, Ala88, Met89, and Gly92 (α-helix 5); and Trp114 (α-helix 6); also, the residues Leu96, Lys930, Arg940, and Leu960 and the two molecules of DEET bound to the dimeric interface of AgamOBP1 interact with each other by the methyl carbon atoms. In addition, studies have demonstrated that odorant molecules could interact in several locations of the OBP binding pockets with partial occupancies, that is, in the central cavity of the monomeric subunit as well as to the OBP dimeric interface.[52−54] Based on these results, we selected for docking analyses the dimeric interface of the DEET-binding pocket located between the two monomers, which is formed by the residues of the helices α4, α5, and α6 (intermolecular interactions and atomic distances obtained in docking are available in Table S3).

The noncovalent interactions formed between the AgamOBP1 binding pocket and the NPs were analyzed over 10 ns of MD simulation. The selected NPs formed numerous hydrophobic interactions with the AgamOBP1 binding pocket, with residues Leu76, His77, Met89, Leu96, and Trp114 repeating interactions similar to those observed with the crystallographic structure of DEET. AgamOBP1 residues Ala88 and Trp114 formed π–alkyl interactions with the aromatic ring of the selected compounds, and some residues, such as Gly92, Leu73, and Met89 (chain A) as well as Lys93 and Leu96 (chain B), formed hydrogen bonds with the oxygenated groups (see Table S4). We also noted that water molecules interact with DEET, with the residues located in the binding pocket, such as Trp114, Asp78, Gly92, and Ser79, and with the main chain of residues Cys95 and Gly92.

An overview of the intermolecular interactions formed between the residues of the AgamOBP1 binging pocket and the selected NPs and docked DEET is shown in Figure and compared with a structural superposition with the crystallographic DEET. The diagram showing the interactions of the six NPs with the AgamOBP1 binding pocket over 10 ns of MD simulation is shown in Figure . Thus, we then performed binding energy calculations to analyze the binding affinities of the selected NPs complexed to the AgamOBP1 pocket.

Figure 3

Superposition of the average structure obtained from MD simulation of each NP and the docked DEET with the crystallographic structure: (A) p-cymen-8-yl, (B) 4-(4-methylphenyl)-pentanal, (C) p-anisyl hexanoate, (D) docked DEET, (E) thymyl isovalerate, (F) carvacryl acetate, (G) thymol acetate. Crystallographic DEET (PDB: 3N7H) is shown in light blue, the residues are in green, and the docked ligands are in gray. Residues from chain B are indicated with a prime (′).

Figure 4

Noncovalent interactions of the selected NPs against the AgamOBP1 binding pocket: (A) 4-(4-methylphenyl)-pentanal, (B) thymol acetate, (C) thymyl isovalerate, (D) carvacryl acetate, (E) p-cymen-8-yl, (F) p-anisyl hexanoate. Residues from chain B are indicated with a prime (′).

Conformational Stability of the AgamOBP1–Ligand Complexes and the Binding Affinities of the Natural Products

The conformational stability of AgamOBP1 bound and unbound to the ligands over 100 ns of MD simulation is shown by the RMSD plots (Figure ). We observed that the AgamOBP1 heterodimer complexed with the ligands reached equilibrium at 70 ns of MD simulation, exhibiting the following average RMSD values: 1.45 ± 0.21 Å (complexed with DEET), 1.74 ± 0.29 Å (complexed with thymol acetate), 1.74 ± 0.32 Å (thymyl isovalerate), 1.61 ± 0.23 Å (carvacryl acetate), 1.43 ± 0.20 Å (p-cymen-8-yl), 1.93 ± 0.43 Å (p-anisyl hexanoate), and 1.49 ± 0.21 Å (4-(4-methylphenyl)-pentanal). Corroborating these results, the RMSD plot of the NP structures exhibited stability over the simulations and indicated a favorable interaction with the AgamOBP1 binding pocket, which is also observed for the DEET structure (Figure S1).

Figure 5

RMSD plots of the AgamOBP1 structure bound and unbound to the ligands obtained over 100 ns of MD simulation. (A) AgamOBP1 unbound to the ligands (black), OBP1 complexed with thymyl isovalerate (yellow), p-anisyl hexanoate (blue), and carvacryl acetate (green). (B) AgamOBP1 complexed with DEET (red), thymol acetate (yellow), p-cymen-8-yl (blue), and 4-(4-methyl phenyl)-pentanal (purple).

During the first MD simulation of the AgamOBP1-DEET complex (positive control), we noted that the protein structure undergoes a conformational change at 60 ns, which can be seen by the RMSD values with a deviation of 3.0 Å (see Figure S2 for triplicate MD simulation of the AgamOBP1–DEET complex). These increased RMSD values could be explained by the movement of a loop segment located in the N-terminal region. Using structural superposition (Figure S3A,B), we analyzed the conformations of the AgamOBP1–DEET complex before (60 ns, yellow color) and after (70 ns, blue color) the increased RMSD values, and we noted that the protein structure maintains slight deviations in its conformation, except for a loop segment (L1) located at the N-terminal region. The loop segment L1 moves toward the α-helix, leading to the formation of intermolecular interactions that stabilize the whole protein structure (Figure S3 C). The conformational changes in the loop do not alter the AgamOBP1 active site, and the DEET structure maintains a stable interaction with residues of the binding site.

In the present study, we calculate the binding free energy using SIE method.[55] As seen from computational binding results, thymyl isovalerate (−7.34 ± 0.09 kcal mol–1), p-anisyl hexanoate (−6.85 ± 0.10 kcal mol–1), and 4-(4-methylphenyl)-pentanal (−6.81 ± 0.08 kcal mol–1) are the most efficient inhibitors of the OBP1 odorant pocket followed by p-cymen-8-yl, carvacryl acetate, and thymol acetate (Table and Table S5). Note that we have used PCA for analyzing the convergent trajectories and selected the time interval for the binding free energy calculation. The analysis of convergent trajectories for each complex studied over the 100 ns of MD is depicted in Figure S4. It is also worth noting that a previous study obtained similar results using docking energies for DEET (ΔG = −5.86 kcal mol–1 ; Ki = 50.51 μM) and potential repellents of AgamOBP1, such as 2-methyl-1-(1-oxodecyl)piperidine (ΔG = −7.36 kcal mol–1; Ki = 4.04 μM), 1-(1-oxoundecyl)piperidine (ΔG = −7.20 kcal mol–1; Ki = 5.28 μM), and N,N-diethyl-3-phenylpropanamide (ΔG = −6.49 kcal mol–1; Ki = 17.21 μM).[39] Similarly, a previous study performed binding free energy calculations using the MM/GBSA method for β-caryophyllene (sesquiterpene), β-myrcene (monoterpene), and cis-β-ocimene (monoterpene) complexed to different OBP classes of Hydroides elegans, and their energy values approximated our calculated values of OBP1,[13] which is also consistent with the inhibitory activity.

Table 2

Predicted Binding Free Energies Calculated by SIE of the Selected NPs from Essential Oil with the AgamOBP1 Binding Pocket

compounds	ΔGSIE (kcal mol–1)
DEET	–6.85 ± 0.13
carvacryl acetate	–6.76 ± 0.12
thymyl isovalerate	–7.34 ± 0.09
thymol acetate	–6.66 ± 0.07
4-(4-methylphenyl)-pentanal	–6.81 ± 0.08
p-anisyl hexanoate	–6.85 ± 0.10
p-cymen-8-yl	–6.75 ± 0.10

The crystallographic structure of AgamOBP1 complexed with DEET reveals that the majority of residues located at chain A, such as Leu73, Leu76, His77, Ala88, Met91, Gly92, Lys93, and Leu96, form hydrophobic contacts with the ligands.[39] In the present study, we noted that the NPs also showed a similar binding mode, and the residues Leu73, His77, Leu80, Ala88, Gly92, Leu96, and Trp144 exhibited the most energetic contribution for ligand stabilization in the AgamOBP1 binding pocket (Figure ). In contrast, the residues of AgamOBP1 chain B exhibited a minor contribution to ligand stabilization, and we noted a special interaction with Met89, Lys93, and Leu96 that was revealed as the most important to ligand binding.

Figure 6

Ligand pairwise per-residue energy decomposition analysis of AgamOBP1 binding pocket. ΔG values of (A) AgamOBP1 chain A and (B) AgamOBP1 chain B: DEET (blue), thymol acetate (yellow), p-cymen-8-yl (light green), carvacryl acetate (red), 4-(4-methylphenyl)-pentanal (cyan), thymyl isovalerate (green), and p-anisyl hexanoate (dark blue).

Structural waters have been reported as relevant for the binding stability of OBP-odorant molecules and to the OBP recognition process.[39] Our analyses identified that water molecules interact with DEET (DEET-C1 with Wat1815-O; average distance: 4.06 ± 0.20 Å) and with the side chains of residues located in the binding pocket, Trp114 (occupancy: 77.60%, Trp114-NE1 atom with Wat1815-O; average distance: 1.97 ± 0.15 Å), and minor interactions were noted for the residues Asp78 (occupancy: 0.25%), Gly92 (0.15%), and Ser79 (0.20%) and with the main chain of the residues Cys95 (3.65%) and Gly92 (0.10%). These interactions with water molecules are consistent with the crystallographic findings. Considering the NPs, we do not note H-bond interactions with water molecules, but similar interactions occurred with residues of the AgamOBP1 binding pocket when bound to thymol acetate, p-anisyl hexanoate, and thymyl isovalerate. The non-occurrence of interactions between water molecules and ligands could be explained by the different occupancies formed by these NPs in the AgamOBP1 binding pocket, which fill the same space previously occupied by the water.

Our results revealed new structural insights about the mechanism of action of these volatile compounds sourced from aromatic plants against the olfactory recognition of mosquitoes, and the findings indicate that the structures of these NPs could be further validated against AgamOBP1.

Conclusions

We have demonstrated that NPs from essential oils, such as thymyl isovalerate, thymol acetate, p-anisyl hexanoate, and p-cymen-8-yl, mimic the binding mode of DEET, a well-known repellent, forming hydrophobic and H-bond interactions with the AgamOBP1 binding pocket. Our results are supported by computational evidence, such as (1) structural, pharmacophoric, and binding mode similarity of these compounds with that of DEET as verified through structural alignment, pharmacophoric prediction, and molecular docking; (2) high conformational stability of these compounds in the OBP1 binding pocket as analyzed through MD simulations; and (3) high predicted binding affinity of the NPs when compared with the DEET–OBP1 complex as revealed by crystallographic structure. We also found that thymol acetate, 4-(4-methyl phenyl)-pentanal, thymyl isovalerate, and p-cymen-8-yl share similar physicochemical properties with the commercial repellents DEPA and DEET. Considering that computational methods have been a cost-effective and predictive approach for analyzing the binding affinity of olfactory receptors with the odorant molecules[7,56] and that they have consequently been applied to the identification of novel potential repellents against different insect species,[57] our findings reinforce the potentiality of NPs as an interesting source for the development of novel mosquito repellents. Further experimental studies should be performed with the compounds p-anisyl hexanoate and 4-(4-methyl phenyl)-pentanal.

Materials and Methods

An overview of the structure- and ligand-based virtual screening methodology applied in the present study is shown in Figure and described in detail in the sections below.

Figure 7

Overview of the applied computational methodology.

Search for 3D Structural Similarity

First, we analyzed the structural similarity of the NPs from essential oils with the DEET crystallographic conformer using Screen3D (ChemAxon).[58] A Tanimoto coefficient cutoff value of 0.7 was applied as a measure of 3D similarity for selecting similar compounds from the initial library with 1633 compounds from the EssOilDB.[38] Screen3D automatically generated the 3D conformers; thus, we limited the maximum number of conformers per compound to 4. The Tanimoto coefficient cutoff is defined by eq .[59]Here, sA, B denotes the similarity between both compounds (A and B), x means the jth features of compound A, x represents the jth features of compound B, and xx is the feature present in both analyzed compounds.

Pharmacophoric Prediction and Filtering

Next, we used the Pharmit server to screen the NP structures with similar pharmacophoric groups and shape with the DEET conformer complexed to AgamOBP1 (PDB: 3N7H, X-ray structure, resolution: 1.6 Å).[60] We analyzed the superposition of the predicted pharmacophores of the natural compound with two pharmacophoric models of DEET complexed with the AgamOBP1 homodimer. Based on these initial screenings, six compounds were selected for further analysis.

Molecular Docking

The Molegro Virtual Docker (MVD) program, which uses the evolutionary algorithm MolDock,[61] was used to perform molecular docking simulations to analyze the binding mode of the NPs complexed to the AgamOBP1 homodimer. First, to validate the docking protocol, redocking simulations were performed with the DEET structure against the AgamOBP1 homodimer. Root-mean-square deviation (RMSD) values less than or equal to 1.0 Å were considered satisfactory for replicating the ligand binding mode in the crystallographic structure. Next, for the docking simulations, docking grids with radii of 8 to 15 Å, depending on the selected ligand, were positioned in the AgamOBP1 active site between both AgamOBP1 monomers using the spatial coordinates of crystallographic DEET as a reference. Water molecules, ions, and the DEET structure were removed from both OBP1 chains before the docking simulations. Then, we performed the flexible docking protocol of MVD, which consider the flexibility of the ligand and residue side chains.[62] All poses obtained were ranked according to the best superposition with the crystallographic structure of DEET and lower docking energy.

Molecular Dynamics Simulation and Analysis of the AgamOBP1–Ligand Complexes

Using the Amber16 package,[63] MD simulations were performed for the following AgamOBP1 systems: (1) AgamOBP1 unbound to the ligands (ligand-free), (2) AgamOBP1 complexed with DEET (PDB: 3N7H), and (3) the six selected NPs from essential oils. First, the partial atomic charges of the ligands were determined using the restrained electrostatic potential (RESP) protocol[64] through quantum mechanical calculations carried out in the Gaussian 09 program.[65] using the Hartree–Fock method[66] with the 6-31G* basis set.[67] The biomolecular systems were solvated in the tLeap module using an octahedral truncated water box with TIP3P, an explicit solvation model.[68] The distance between the water box wall and the atoms of the solvated system was set to 12 Å. Na+ counterions were added to the water box to maintain electroneutrality. The force fields ff14SB[69] and General Amber Force Field[70] were used to parameterize the protein (AgamOBP1) and the ligand structures (NPs and DEET), respectively. An energy minimization protocol with six steps that included the steepest-descent and conjugated gradient algorithm was performed for all systems. All hydrogen atoms, water, and ions were minimized by 10,000 cycles for each step followed by minimization of the whole AgamOBP1 system with the progressive decrease of restraints. Next, the whole system was heated through 10 heating steps. The first heating step was performed at a constant volume for 20 ps, increasing the temperature to 100 K. From the second to the ninth step, 1 ns was used to raise the temperature gradually from 100 to 275 K. In the 10th heating step, the temperature reached 300 K, which was then maintained using 5 ns of MD simulation to equilibrate the density at constant pressure (1 bar). The temperature was maintained at 300 K by coupling to a Langevin thermostat using a collision frequency of 2 ps–1, and the constant isotropic pressure was maintained at 1 bar by using the Berendsen barostat. All stages of the simulations employed a cutoff of 10 Å for nonbonded interactions, and the particle mesh Ewald (PME) method was used to compute the long-range electrostatic interactions. The SHAKE algorithm was applied for all H bonds during MD simulation, and the time step was set to 2 fs.[71] Finally, MD simulation was performed in 100 ns. The RMSD plot was obtained for each system during the simulation using the heavy atoms of the protein backbone. The noncovalent interactions of the receptor–ligand complex were analyzed over 10 ns of MD simulation using PyContact[72] and visually inspected using the Visual Molecular Dynamics (VMD) program.[73]

Analysis of Conformational States over the MD Simulation

To calculate the binding free energy and to measure the molecular interactions between the ligands and AgamOBP1 binding pocket, we selected the intervals of the MD trajectories based on the convergent conformational states of AgamOBP1–ligand complexes obtained during 100 ns of MD simulation using principal component analysis (PCA). PCA is a transformation technique that converts a series of potentially coordinated observations present in a covariance matrix, reducing the linear correlations among them, thus transforming into a set of orthogonal vectors named principal components (PCs).[74] The first component (PC1) maximizes the variance data in the data set and the rest of the variance is represented by the second (PC2) and third component (PC3), respectively. The PCA combined with MD simulations of proteins has been widely applied to analyze the local and collective movements of protein structures,[75,76] to determine the conformational changes that favor enzyme catalysis,[77] to explore the functional roles of ion binding in protein structure,[78] and to sample the convergence and correspondence between the protein structures over the MD simulation.[79] Our conformational analyses were performed in the CPPTRAJ[80] and Bio3D[81] programs, and the scatter plots were build using PC1 and PC2.

Binding Free Energy Calculations

To compute the binding free energy between the ligands (NPs and DEET) and the AgamOBP1 structure, we selected three intervals with 1000 frames from each MD trajectory. To calculate the binding free energy, we applied the solvated interaction energy (SIE) method[55] available in SIETRAJ.[82] We also performed a ligand pairwise per-residue energy decomposition analysis using the molecular mechanics/generalized Born surface area (MM/GBSA) method[83] available in Amber16.[63]

Analysis of the Chemical Space of Natural Products from Essential Oils

To analyze the chemical space of the filtered NPs obtained from essential oils, we compared them with 18 approved or experimental repellents (Table S6) using PCA. Six physicochemical and structural descriptors were used in PCA: number of rotatable bonds, c log P, molecular weight, number of H bond donors, H bond acceptors, and number of aromatic rings. These descriptors have been applied to determine the herbicide-, fungicide-, pesticide-, and insecticide-likeness of compounds.[84] All properties were calculated using the Instant JChem suite.[58]