In Silico Studies on Compounds Derived from Calceolaria: Phenylethanoid Glycosides as Potential Multitarget Inhibitors for the Development of Pesticides

An increasing occurrence of resistance in insect pests and high mammal toxicity exhibited by common pesticides increase the need for new alternative molecules. Among these alternatives, bioinsecticides are considered to be environmentally friendly and safer than synthetic insecticides. Particularly, plant extracts have shown great potential in laboratory conditions. However, the lack of studies that confirm their mechanisms of action diminishes their potential applications on a large scale. Previously, we have reported the insect growth regulator and insecticidal activities of secondary metabolites isolated from plants of the Calceolaria genus. Herein, we report an in silico study of compounds isolated from Calceolaria against acetylcholinesterase, prophenoloxidase, and ecdysone receptor. The molecular docking results are consistent with the previously reported experimental results, which were obtained during the bioevaluation of Calceolaria extracts. Among the compounds, phenylethanoid glycosides, such as verbascoside, exhibited good theoretical affinity to all the analyzed targets. In light of these results, we developed an index to evaluate potential multitarget insecticides based on docking scores.


Introduction
The continuous growth of the world population has created an enormous pressure to satisfy the global demand for agricultural products. The challenges include the depletion of soil fertility, the constant depredation of natural soils to convert them into agricultural ecosystems, and the ability of the arthropods to obtain resistance against traditional insecticidal controls.
Insecticides have been used for combating insect pests, mainly to increase the yield of food production among other agricultural products. From ancient times, there are records that describe the use of different types of products to combat insect pests [1]. It is known that various nonspecific agents have been used, such as sulfur and poisonous natural extracts, then organochlorines, organophosphates, carbamates, pyrethroids, and rotenoids, among others, and finally compounds, which are specifically designed and synthesized against enzymatic systems of arthropods.
The increasing need for agricultural goods has resulted in misutilization of insecticides, and this has led to the use of a higher concentration of insecticides or to the need for more toxic products. This has resulted in increased toxic effects on other beneficial organisms that coexist with pests in agroecosystems and on the bioaccumulation of higher concentrations of toxic insecticides in the bodies of predators or the final consumers, including humans. Despite these problems, the use of insecticides Despite these limitations, DBVS has been successfully used in the identification of potential templates for new drug development [16][17][18]. Recently, some examples of the use of virtual screening and other computational chemistry tools in natural products research have been described [19][20][21][22].
With this in mind, we wanted to determine the potential use of compounds isolated from Calceloraia as leads in the discovery of multitarget insecticides using DBVS on some proteins recognized as targets for pesticides [12,23,24] and that could be targets for compounds present in Calceolaria extracts based on experimental results [9,25,26]: acetylcholinesterase (AChE), prophenoloxidase (PPO), and ecdysone receptor (EcR). Construction of a ligand library was based on compounds isolated and present in organic extracts with experimentally demonstrated pesticide activity aiming to identify potential structural templates that could be used in the development of new pesticides.

Ligand Construction
All of the ligands were chosen from a previously published review [27], which includes several compounds from different chemical families, including diterpenes, triterpenes, and naphthoquinones with a potential pesticidal activity, and some bioactive flavonoids and phenylethanoid glycosides as well. All of the structures were constructed using Spartan '10 for Windows, and these geometries were optimized using the MMFF force field. Then, these structures were exported to Molegro Virtual Docker 6.0.1 [28]; assignments of charges and ionization were based on standard templates as part of the Molegro software. A complete list of all ligands and their structures are presented in Supplementary Information.

Molecular Docking Studies
The docking studies were carried out based on the crystal structures of Drosophila melanogaster acetylcholinesterase (DmAchE, PDB codes: 1DX4 [29] and 1QON [29]), Heliothis virescens ecdysone receptor (EcR, PDB code: 3IXP [30] and 2R40 [31]), and Manduca sexta prophenoloxidase (PPO, PDB code: 3HHS [32]). Two different structures from the Protein Data Bank (PDB) were selected to analyze repeatability of results, independent of the PDB structure selected. This was not possible for PPO as no other PDB structure has been reported. In addition, docking studies were carried out in human acetylcholinesterase (hAChE, PDB code: 4EY7 [33] and 4M0E [34]) to determine whether some of the compounds exhibited theoretical preference to the Drosophila/human enzyme. All structures were retrieved from the Protein Data Bank [35]. Docking studies were carried out using a previously reported methodology [36,37]. Briefly, all of the solvent molecules and cocrystallized ligands were removed from the structures. Molecular docking calculations for all of the compounds with each of the proteins were performed using Molegro Virtual Docker v. 6.0.1 [28]. Active sites of each enzyme were chosen as the binding sites and delimited with a 15 Å radius sphere centered on the cocrystallized ligand, except for the PPO structure, which has no cocrystallized ligand, and the sphere was centered on the active Cu 2+ ions. Standard software procedure was used. The assignments of charges on each protein were based on standard templates as part of the Molegro Virtual Docker program, and no other charges were necessary to be set. The Root Mean Square Deviation (RMSD) threshold for multiple cluster poses was set to <1.00 Å. The docking algorithm was set to 5000 maximum iterations with a simplex evolution population size of 50 and a minimum of 25 runs for each ligand. After docking, MolDock Score was calculated as the theoretical binding affinity. To assess the efficacy of this procedure, cocrystallized ligands were also docked to their respective receptors (except for PPO), the top-ranking score was recorded, and the RMSD of that pose from the corresponding crystal coordinates was computed. In all the cases, the RMSD was lower than 2 Å. For each enzyme, the 10 compounds with lower MolDock scores were selected for analyzing their docking poses to identify potential structural requirements for enzyme binding.

Molecular Dynamics Simulations
Molecular Dynamics (MD) simulations were carried out to observe differences that could account for potential selectivity of phenylpropanoids for DmAChE over hAChE. Simulations were performed in YASARA Dynamics v.18.4.24 [38,39] using AMBER14 force field [40]. The initial structures for the MD simulation were obtained from the docking complexes of compound 87 with DmAChE (PDB code: 1QON) and with hAChE (PDB code: 4M0E). Compound 87 (verbascoside) was selected, as it is the most studied compound of the phenylpropanoids derived from Calceolaria. Each complex was positioned into a water box with a size of 100 Å × 100 Å × 100 Å, with periodic boundary conditions. Temperature was set at 298 K, water density to 0.997 g/cm 3 , and pH to 7.4. Sodium (Na + ) and chlorine (Cl − ) ions were included to provide conditions that simulate a physiological solution (NaCl 0.9%). Particle mesh Ewald algorithm was applied with a cut-off radius of 8 Å. A timestep of 2.5 fs was set. The simulation snapshots were recorded at intervals of 100 ps until a total simulation time of 30 ns. Results were analyzed with a script included as part of YASARA software and included RMSD, ligand binding energy variations (using MM-PBSA calculations), and distance of ligand 87 to Ser 283 (DmAChE) or Ser 203 (hAChE), as these residues play a key role in acetylcholinesterase enzymatic activity. For DmAChE, we considered the interaction between the primary alcohol group of the central glucopyranose ring with Ser 238, and for hAChE, the interaction between the hydroxyl group of the ferouyl residue with Ser 203. A similar procedure has been recently reported for the simulation of complexes of drugs with some proteins [41][42][43].

Construction of the Virtual Multitarget Index and the Weighed Multi-Target Index
The virtual Multitarget index of each compound was determined. To compare the multitarget index of the analyzed compounds, we propose a virtual multitarget index (vMTi), which was calculated for the three insect targets (EcR, PPO, and DmAChE) using formula (1): where MDi corresponds to the MolDock score of the molecule in a specific target and MDr is the MolDock score of the reference ligand; we considered the compound with the lowest MolDock score (which has the highest theoretical affinity) for each target as the reference. Compounds with higher values have a higher multitarget index. However, binding to hAChE is an undesirable condition.
To take this into consideration, we propose a weighed MTi (wMTi), which was calculated using formula 2, using an external coefficient n, which represents the desirability of binding to a specific target: To calculate wMTi, we gave values of n = 0.3 to desirable targets (EcR, PPO, and DmAChE) and n = −0.3 to hAChE.
In addition to these calculations, a contour plot was built with Minitab using PPO, DmAChE, and EcR docking scores. This plot can help identify potential multitarget compounds because those compounds would appear in valleys since they would have lower MolDock scores.

Docking Studies Results
Tables 1-3 show data for the 10 compounds with a lower average MolDock score in the ecdysone receptor (EcR), prophenoloxidase (PPO), and acetylcholinesterase (both DmAchE and hAChE) docking study. A table with complete docking results is presented in Supplementary Information.  We also wanted to check if the distance of compound 87 (verbascoside) to the catalytic site of acetylcholinesterase variates during the simulation time. We selected the potential interactions to Ser 203 in hAChE or Ser 238 in DmAChE predicted by molecular docking as described in Methodology. Figure 1b shows the variation of distances of compound 87 to these key serine residues along the simulation time. As seen in this figure, the distance to catalytic site diminished from 4.8 Å to 3.0 Å after 3 ns of simulation in the case of DmAChE (average distance = 3.6 Å), while it maintained almost the same in hAChE (average distance = 5.21 Å). To estimate the difference in binding energy of compound 87 in its complex with DmAChE and hAChE, MM-PBSA methods were applied. Ligand binding energy suggests better binding of compound 87 to DmAChE (E = −131.5 kJ/mol) versus hAChE (E = −134.0 kJ/mol) as binding energy calculations implemented in YASARA Dynamics indicates that the higher the energy value, the better the binding. Figure 2 shows a contour graphic that compares docking scores of the analyzed compounds on PPO, DmAChE, and EcR. The zone in red corresponds to those molecules with high theoretical affinity against all the three molecular targets. Figure 3 shows the structure of the phenylpropanoids which resulted with the highest values of vMTi and wMTi. Table 4 shows a list of compounds with higher vMti and wMTi values; a full list is included in Supplementary Information.

Docking Studies on Ecdysone Receptor
Induction of molting in Arthropods coincides with a release of 20-hydroxyecdysone (20-E), a steroidal-type hormone. Prior to each of the larval molts, at pupariation, at pupation, and during metamorphosis, the hormone is released in carefully timed spurts, coinciding with major morphological transitions. Ecdysone receptor (EcR) exists in three isoforms. Each requires a partner during the heterodimerization, a Drosophila homolog of vertebrate RXR protein named ultraspiracle (USP) protein. Although ecdysone can bind to EcR on its own, binding is significantly augmented by the participation of USP. The interaction between EcR and ecdysone is a crucial event for the development of insects, which is why it represents an interesting molecular target against pest insects [44]. Table 1 shows data for the 10 compounds with a lower average MolDock score in the ecdysone receptor docking study. Phenylethanoid glycoside derivatives have better theoretical binding to EcR among the evaluated compounds. Among them, the results indicate that binding to EcR occurs through hydroxyl groups of caffeoyl and phenyl ethyl residues in similar fashion to 20-E (Figure 4a). Analysis of Figure 4a,b reveals that compound 87 (verbascoside, which is a major phenylethanoid glycoside present in Calceolaria extracts) adopted a J-shaped conformation in the binding site of EcR, which is similar to the conformation adopted by natural ligand of EcR, 20-E [31] and the interaction pattern is very similar between these compounds: phenylethanoid residue interacts through hydrogen bonds with Arg383 and Glu309 (in magenta in Figure 4a,b) like 2β and 3β hydroxyl groups of ring-A in 20-E, rhamnose residue interacts with Ala398 (in blue in Figure 3a,b) like C-6 ketone moiety in ring-B in 20-E, feruoyl residue interacts with Thr 343 and rhamnose residue with Thr 346 (also in blue in set of Figure 4) like 14-α hydroxyl group, and additional interactions between feruoyl residue and central glucose ring of verbascoside with Tyr 408 and Asn 504 (shown in yellow) are seen in a similar fashion for 25-OH group of side chain of 20-E. This interaction is notable because 20-E interacts with this residue via a water linkage [31], hence the interaction of verbascoside with Asn504 could increase affinity, as some studies have demonstrated that ligands designed to displace the water molecules exhibit higher affinity [45]. Though there are no previous studies on phenylethanoid glycosides as EcR ligands, there is some experimental evidence that can support docking results. It has been previously reported that the ethyl acetate extract of C. talcana, of which verbascoside is its major component (compound 87), caused a developmental disruption of D. melanogaster and S. frugiperda larvae. In addition, the authors of the study proposed verbascoside as a disruptor of ecdysteroid metabolism [26]. In another study, the incorporation of verbascoside in the artificial diet of the pest Agrilus planipennis caused a 100% mortality at 45 mg/g of artificial diet [46] when testing the toxic effect of verbascoside against at least three different pest insect species. In addition, Harmatha and Dinan have reported that some polyhydroxylated stilbenoids have an antagonist EcR activity [47], and a previously reported pharmacophore model indicated that the presence of hydroxyl groups in an ecdysteroid template is important for EcR binding [48]. On the other hand, the presence of some other phenylethanoid glycosides, such as calceolariosides A, B, and C, in active extracts of Calceolaria as well as in Fraxinus spp. can be related to the strong larval molting disruption observed when the larvae of different species were exposed to extracts with high amounts of phenylethanoid glycosides or directly to different amounts of the isolated compounds [46]. All these experimental data strongly suggest that phenylpropanoid glycosides could be EcR ligands.
As described here, the effect of verbascoside and related phenylethanoid glycosides against the ecdysone receptor can explain only one of the multiple effects exerted by these compounds. The experimental evidence remarkably indicates that molting disruption exerted by phenylethanoid glycosides cannot be the only mechanism that explains the strong exerted insecticidal properties. The above information suggests a possible antagonist and multienzymatic inhibitory mechanism that phenylethanoid glycosides can exert, causing larvae disruption activity by acting as EcR antagonist in addition to other mechanisms.

Docking Studies on Prophenoloxidase
Melanization, a process performed by phenoloxidase (PO) and controlled by the prophenoloxidase (PPO) activation cascade, plays an important role in the invertebrate immune system in allowing a rapid response to pathogen infection. The activation of the PPO system, by the specific recognition of microorganisms by pattern-recognition proteins (PRPs), triggers a serine proteinase cascade, which eventually leads to the cleavage of the inactive PPO to the active PO that functions to produce melanin and toxic reactive intermediates. The importance of PPO-PO is due to cuticular sclerotization and defense against pathogens and parasites. PO catalyzes hydroxylation of monophenols to o-diphenols and oxidation of o-diphenols to quinones. Quinones take part in sclerotization and tanning of the cuticle and serve as precursors for synthesis of melanin [49,50]. Therefore, PPO is a very suitable molecular target for designing pesticide compounds. Table 2 shows docking results from the PPO study. Among the compounds with better PPO binding, the best are flavonoids and phenylethanoid glycosides. The results are in agreement with several previous in vitro and in silico studies [10,[51][52][53][54][55] in which most tyrosinase inhibitors possess a phenol moiety as the pharmacophore. Among them, flavonoids appear as effective competitive inhibitors of this enzyme. In addition, Karioti et al. [53] and Muñoz et al. [10] have reported tyrosinase inhibitory activities of some phenylethanoid glycosides, which have lower but comparable inhibitory activities compared with flavonols and flavones.
Hydroxyl groups of caffeoyl or phenylethyl residues have been proposed as essential structural requirements to display inhibitory activities against PPO because of their chelating properties. Figure 5 shows that phenylethanoid glycosides could bind to the PPO catalytic site through interaction with His residues, which are required to form a complex with copper ions. This is in agreement with the previously reported information that indicates verbascoside as a substrate for tyrosinase [10]. Among phenylethanoid glycosides, compounds 86 and 93 could bind better than other analogs. These compounds are monoglycosides, whereas other evaluated phenylethanoid glycosides are diglycosides, which is in agreement with previous reports that indicate the increase in the number of sugar units and the reduction of PPO inhibitory activity [53]. This observation could be explained in terms of the higher molecular volume of diglycosides that prevents them access to the active site as is shown Figure 5; the diglycoside compound 88 is shown in yellow and monoglycoside 86 is in cyan. From this figure, it can be concluded that it is possible that monoglycosides could bind closer to the catalytic site of PPO.

Docking Studies and Molecular Dynamics Simulations on Drosophila and Human Acetylcholinesterase
Acetylcholinesterase is a serine hydrolase that is vital for regulating the neurotransmitter acetylcholine in insects. This enzyme is an excellent molecular target for the development of insecticides [56,57]. The well-known active site has a deep and narrow gorge, with a catalytic site at the bottom and a peripheral site at the entrance. Acetylcholinesterase is a molecular target used to control insects that affect public health (e.g., mosquitoes, flies, cockroaches, among others) as well as those that affect agriculture and gardening (e.g., grasshoppers, aphids, caterpillars, among others) [58]. Current anticholinesterase insecticides work through phosphorylation of a serine residue at the AChE catalytic site, which disables the catalytic function and causes enzyme antagonism. Because this serine residue is also ubiquitous in AChEs of mammals and other species with cholinergic nerves, the use of anticholinesterase insecticides to target the serine residue causes serious off-target toxicity [59]. Therefore, it is necessary to evaluate in silico the affinity of molecules against the insect together with mammalian enzymes to determine whether there are some structural features that lead to design inhibitors specifically against insect enzymes. Table 3 displays the MolDock scores of the top 10 compounds with better average docking scores for DmAChE studies in PDB 1DX4 and 1QON structures. Data obtained during docking studies with hAChE (PDB structures 4EY7 and 4M0E) is also shown. Selectivity ratio versus hAChE was calculated with the average docking score obtained in both enzymes. Although some differences could be appreciated in docking scores values, the same tendency was observed in both DmAChE structures and both hAChE structures as phenylethanoid glycosides are among the compounds with a better theoretical binding in the four docking studies. It has been reported that verbascoside and extracts containing other phenylethanoid glycosides are moderate inhibitors of AchE [25]. Figure 6a,b shows the predicted binding mode of verbascoside (compound 87) in the active site of both DmAchE and hAChE, respectively. In the case of the docking study carried out in DmAChE, verbascoside and the rest of the analyzed PEGs adopted a Y-shaped conformation, with the central sugar core interacting with residues of the catalytic triad (colored in yellow in Figure 6a) and the phenylethoxy chain interacting with other residues within the bottom of the gorge. In the case of hAChE, the analyzed PEGs adopted a similar conformation, but the phenylethoxy chain did not reach the bottom of the gorge. An explanation to this is that PEGs interact with Asp74 through a hydrogen bond in hAChE, whereas this residue is absent in DmAChE [60], limiting the access of PEGs to the catalytic triad, and this could explain the better theoretical affinity of PEGs to DmAChE compared to hAChE.
To analyze additional differences that could account for potential selectivity of PGs to DmAChE over hAChE, MD simulations were performed. The RMSD values could indicate the stability of the protein relative to its conformation. Figure 4 shows the comparison of the RMSD of the protein backbone profile during the 30 ns simulation of the complexes of compound 87 and DmAChE and hAChE. Both complexes have RMSD values around 2Å (RMSD = 2.13 for DmAChE and RMSD = 1.88 for hAChE), suggesting that both complexes are stable.
An important difference that was observed during visual analysis MD simulations was the variation of distance of verbascoside (compound 87) to catalytic site in DmAChE, as it appeared to move closer to key residues Glu 237 and Ser 238, while it seemed that compound 87 maintained a constant distance to equivalent residues Glu 202 and Ser 203 in hAchE complex. We confirmed this by measuring the distance of verbascoside to Ser 238 or Ser 203. As seen in Figure 2b, the distance to catalytic site diminished from 4.8 Å to 3.0 Å after 3 ns of simulation in the case of DmAChE (average distance = 3.6 Å), while it maintained almost the same in hAChE (average distance = 5.21 Å). This could account for the better binding of 87 towards DmAChE. This could be an important difference that could be useful to future design of selective inhibitors to DmAChE based on the phenylethanoid glycoside template. Additionally, MM-PBSA calculations suggest that 87 binds better to DmAchE (E = −131.5 kJ/mol) than to hAChE (E = −134.0 kJ/mol). Overall, we conclude that, based on molecular docking calculations and MD simulations, compound 87 is an interesting starting point for the design of selective DmAchE inhibitors, an important factor to consider in terms of potential toxicity against human beings.

Virtual Multitarget Index and Weighed Multitarget Index
A multitarget drug has been defined as the integration of multiple pharmacophores in one molecule with the purpose that it can have two or more simultaneous mechanisms of action [61]. Though the concept of multitarget drugs is an important research topic [62,63], there are few examples of the study of multitarget insecticides [12].
Computational tools like molecular docking and multitarget quantitative structure-activity relationship models (mt-QSAR) have been recently used for prediction and discovery of multitarget compounds [64][65][66][67]. However, the development of parameters to measure the multitarget index of a ligand is not an easy task. As described during the lapatinib discovery [68], the simple average of biological activity against several targets (or in our case, docking scores) can be misleading, because a compound with a promising average value of bioactivity on two or more targets could correspond to a multitarget compound or to a highly selective compound against only one target. Thus, a reference parameter and weighed coefficients for biological activities of interest should be included. Herein, we propose the use of contour graphics, such as the one shown in Figure 2, and a multitarget index (vMTi) to identify potential multitarget insecticides. In the contour plot of Figure 2, the docking scores in three targets of interest (EcR, PPO, and DmAChE) are shown. Compounds that have shown greater theoretical affinity for the three targets will appear in the bottom left of the plot and inside the red and yellow contour areas. This would be a first criteria to identify potential multitarget compounds, because selective compounds would not appear in this area. As expected, the phenylethanoid glycosides appear in this area, but other compounds, such as abietane 3 and isopimarane 77, can be considered to have potential multitarget properties. Table 4 shows a list of compounds with higher vMti and wMTi values. In this table, phenylethanoid glycosides appear as compounds with a better multitarget profile. In addition, compounds 88, 89, and 87 ( Figure 2) have not only high vMTi values, but also the highest wMTi due to their higher selectivity to DmAChE than hAChE. Thus, these compounds are interesting candidates for the development and evaluation of safer multitarget insecticides.

Conclusions
Our results complement the experimental results obtained during Calceolaria extracts evaluation as biopesticides, and suggest that some of the compounds, such as the phenylethanoid family, can be used for the development of multitarget bioinsecticides. Based on the docking studies, it appears that verbascoside and other phenylethanoid glycosides could exert their bioactivity by modifying the activity of various receptors like EcR and enzymes like PPO and AChE, as was suggested and confirmed in previous experimental assays [25,26]. Theoretical affinity, together with vMTi and wMTi, can be useful for the rational design of multitarget bioinsecticides. Verbascoside appears as a good candidate for the development of a multitarget insecticide due to its prolific natural occurrence and its chemical and biological properties [69].

Acknowledgments:
The authors thank Dirección de Posgrado e Investigación of Universidad La Salle for access to additional computational resources.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.