Larvicidal Activity of Cinnamic Acid Derivatives: Investigating Alternative Products for Aedes aegypti L. Control

The mosquito Aedes aegypti transmits the virus that causes dengue, yellow fever, Zika and Chikungunya viruses, and in several regions of the planet represents a vector of great clinical importance. In terms of mortality and morbidity, infections caused by Ae. aegypti are among the most serious arthropod transmitted viral diseases. The present study investigated the larvicidal potential of seventeen cinnamic acid derivatives against fourth stage Ae. aegypti larvae. The larvicide assays were performed using larval mortality rates to determine lethal concentration (LC50). Compounds containing the medium alkyl chains butyl cinnamate (7) and pentyl cinnamate (8) presented excellent larvicidal activity with LC50 values of around 0.21–0.17 mM, respectively. While among the derivatives with aryl substituents, the best LC50 result was 0.55 mM for benzyl cinnamate (13). The tested derivatives were natural compounds and in pharmacology and antiparasitic studies, many have been evaluated using biological models for environmental and toxicological safety. Molecular modeling analyses suggest that the larvicidal activity of these compounds might be due to a multi-target mechanism of action involving inhibition of a carbonic anhydrase (CA), a histone deacetylase (HDAC2), and two sodium-dependent cation-chloride co-transporters (CCC2 e CCC3).


Introduction
Arbovirus infections are caused by viruses transmitted to people through the bite of an infected arthropod, such as the mosquito. Dengue, yellow fever, Zika and Chikungunya are all transmitted by Aedes aegypti L., which is of great epidemiological importance as the main vector of arboviral infections in tropical and subtropical regions [1,2]. Of the above-mentioned diseases, dengue is the most prevalent, being present on five continents, and affecting about 50 million people per year around the world [3]. Dengue presents great potential to assume severe and lethal forms [4].
Most of the derivatives tested (except compounds 14, 15, 16 and 17) demonstrated greater larvicidal potency than cinnamic acid (1) the starting material. The derivatives with alky substituents (compounds 2 through 11) presented LC 50 values ranging from 0.17 to 3.54 mM (p < 0.05). Methyl cinnamate (2) presented high potency compared to cinnamic acid (1); and ethyl cinnamate (3) with an ethyl group in its carbon chain presented increased larvicidal effect (LC 50 : 0.74 mM). However, these results are at odds with data from the literature in which it is reported that the presence of one or more alkyl groups in the chemical structure of cinnamic acid ester radicals lowers the intensity of larvicidal activity [27,28]. Certain compounds presented stronger larvicidal action, such as isopropyl cinnamate (5), methoxyethyl cinnamate (6), butyl cinnamate (7), pentyl cinnamate (8), isopentyl cinnamate (9) and hexyl cinnamate (10) (LC 50 : 0.52, 0.58, 0.21, 0.17, 0.83 and 0.99 mM, as shown respectively in Figure 1. These compounds present medium substituted alkyl chains such as compounds 7, 8, and 10, a heteroatom in the carbon chain (compound 6), and small to median carbon chain, for example compounds 5 and 9; being the pentyl cinnamate (8) (the most active of all the derivatives tested) presented an LC 50  Most of the derivatives tested (except compounds 14, 15, 16 and 17) demonstrated greater larvicidal potency than cinnamic acid (1) the starting material. The derivatives with alky substituents (compounds 2 through 11) presented LC50 values ranging from 0.17 to 3.54 mM (p < 0.05). Methyl cinnamate (2) presented high potency compared to cinnamic acid (1); and ethyl cinnamate (3) with an ethyl group in its carbon chain presented increased larvicidal effect (LC50: 0.74 mM). However, these results are at odds with data from the literature in which it is reported that the presence of one or more alkyl groups in the chemical structure of cinnamic acid ester radicals lowers the intensity of larvicidal activity [27,28]. Certain compounds presented stronger larvicidal action, such as isopropyl cinnamate (5), methoxyethyl cinnamate (6), butyl cinnamate (7), pentyl cinnamate (8), isopentyl cinnamate (9) and hexyl cinnamate (10) (LC50: 0.52, 0.58, 0.21, 0.17, 0.83 and 0.99 mM, as shown respectively in Figure 1. These compounds present medium substituted alkyl chains such as compounds 7, 8, and 10, a heteroatom in the carbon chain (compound 6), and small to median carbon chain, for example compounds 5 and 9; being the pentyl cinnamate (8) (the most active of all the derivatives tested) presented an LC50 of less than 0.2 mM.
The compounds presenting an aryl radical (compounds 12 to 18) also presented less larvicidal effect than the alkyl cinnamates (compounds 2 to 11); of these compounds, benzyl cinnamate (13), 4chlorobenzyl cinnamate (12), and 3-metoxylbenzyl cinnamate (18) all presented better larvicidal activity than the starting material (1), as shown Figure 1. Compound 13 was the most potent aromatic derivative, with LC50: 0.55 mM [17]; yet the presence of ring substituents resulted in decreased larvicidal potency. Comparing compounds 14 and 15 suggests that both bulky alkyl substituents or non-bulky substituents attached to the aromatic ring in para position decrease larvicidal potency (LC50: 12.06 and 11.21 mM, respectively). The isomers 17 and 18, present a methoxy substituent in different positions of the benzyl ring, respectively para or meta. Biological tests revealed that changing the position of this substituent for compound 18 increased larvicidal potency (LC50: 0.82 mM, Figure  1). Among the other aryl derivatives with electron withdrawing substituents, compounds 12 and 16, only compound 12 showed greater larvicidal activity than the starting material (1), (LC50: 3.01 mM, Figure 1).  Thus, the SAR analysis shows certain structural features do increase larvicidal activity, such as cinnamate alkyl side chain, (especially the median linear carbon chain); or aryl radicals without substituents; or aryl radicals with meta position methoxy substituents. Figure 2 summarizes the study of cinnamic acid derivatives for larvicidal activity. The results are encouraging, and they show that esters of cinnamic acid are structurally simpler than the insecticide temephos, exhibit larvicidal  The less potent compounds (though better than the starting material), presented larger alkyl chains, for example dodecyl cinnamate (11); small alkyl chains such as methyl cinnamate (2) and propyl cinnamate (4).
The compounds presenting an aryl radical (compounds 12 to 18) also presented less larvicidal effect than the alkyl cinnamates (compounds 2 to 11); of these compounds, benzyl cinnamate (13), 4-chlorobenzyl cinnamate (12), and 3-metoxylbenzyl cinnamate (18) all presented better larvicidal activity than the starting material (1), as shown Figure 1. Compound 13 was the most potent aromatic derivative, with LC 50 : 0.55 mM [17]; yet the presence of ring substituents resulted in decreased larvicidal potency. Comparing compounds 14 and 15 suggests that both bulky alkyl substituents or non-bulky substituents attached to the aromatic ring in para position decrease larvicidal potency (LC 50 : 12.06 and 11.21 mM, respectively). The isomers 17 and 18, present a methoxy substituent in different positions of the benzyl ring, respectively para or meta. Biological tests revealed that changing the position of this substituent for compound 18 increased larvicidal potency (LC 50 : 0.82 mM, Figure 1). Among the other aryl derivatives with electron withdrawing substituents, compounds 12 and 16, only compound 12 showed greater larvicidal activity than the starting material (1), (LC 50 : 3.01 mM, Figure 1).
Thus, the SAR analysis shows certain structural features do increase larvicidal activity, such as cinnamate alkyl side chain, (especially the median linear carbon chain); or aryl radicals without substituents; or aryl radicals with meta position methoxy substituents. Figure 2 summarizes the study of cinnamic acid derivatives for larvicidal activity. The results are encouraging, and they show that esters of cinnamic acid are structurally simpler than the insecticide temephos, exhibit larvicidal potential, and are either found in nature or analogous to natural cinnamic derivatives.  Thus, the SAR analysis shows certain structural features do increase larvicidal activity, such as cinnamate alkyl side chain, (especially the median linear carbon chain); or aryl radicals without substituents; or aryl radicals with meta position methoxy substituents. Figure 2 summarizes the study of cinnamic acid derivatives for larvicidal activity. The results are encouraging, and they show that esters of cinnamic acid are structurally simpler than the insecticide temephos, exhibit larvicidal potential, and are either found in nature or analogous to natural cinnamic derivatives.

Computational Methods
To study the potential mechanism of action of the compounds herein presented, modeling studies were performed as described in the Materials and Methods section. The list of the potential targets of the studied compounds is presented in Table 2 and it includes the UniProt [29] ID of the identified targets, the proteins descriptions, the genes ID of the target in the VectorBase [30] database and the ID assigned to each potential receptor throughout the rest of this report. Among the proteins listed in Table 2, GSTE2, GSTE4, GSTD2, GSTT1, GSTX2, GSTT2 and ACE were not predicted during the computational target fishing calculations. Regarding the Glutathione transferase enzymes, these were included in our study because there is evidence suggesting that cinnamic acid can inhibit their activity [31]. ACE was included because, despite there is no evidence pointing to its inhibition by cinnamic acid and its derivatives, it is the most widely studied target of insecticidal compounds in Ae. aegypti [32][33][34].

Computational Methods
To study the potential mechanism of action of the compounds herein presented, modeling studies were performed as described in the Section 3. The list of the potential targets of the studied compounds is presented in Table 2 and it includes the UniProt [29] ID of the identified targets, the proteins descriptions, the genes ID of the target in the VectorBase [30] database and the ID assigned to each potential receptor throughout the rest of this report. Among the proteins listed in Table 2, GSTE2, GSTE4, GSTD2, GSTT1, GSTX2, GSTT2 and ACE were not predicted during the computational target fishing calculations. Regarding the Glutathione transferase enzymes, these were included in our study because there is evidence suggesting that cinnamic acid can inhibit their activity [31]. ACE was included because, despite there is no evidence pointing to its inhibition by cinnamic acid and its derivatives, it is the most widely studied target of insecticidal compounds in Ae. aegypti [32][33][34]. Molecular docking of compound 8 to the 27 receptors listed in Table 2 was performed as described in the Materials and Methods section. The detailed results of this step are provided in Table 3. Docking scores provided by Gold are dimensionless and the higher a score value is, the better the predicted pose is. The visual inspection of the predicted ligandreceptor complexes reveals meaningful complexes showing favorable ligand-receptor interactions and complementarity. This is further supported by the obtained docking scores. Following the same consensus scoring protocol used to select the probable binding modes of compound 8 to its potential receptors, it can be seen that the best consensus docking scores are obtained for the ACE, LGUL, HDAC2, RAC1 and GSTE2 targets. On the other hand, the CCC3, GLUR, GSTT1 and AKR6 targets are the ones with the poorest molecular docking scores. It is also worth noting that more than one possible binding mode is predicted for 22 out of the 27 proteins studied, leading to 60 compound 8-receptor complexes to be analyzed.  Due to the computational complexity of molecular recognition and the need of algorithms capable of processing large amounts of compounds in a reasonable time, many factors relevant for this process are simplified or neglected in molecular docking. It has also been shown that the estimation of the free energy of binding from Molecular Dynamics (MD) simulations snapshots during the post processing of molecular docking models can lead to a better selection of molecular targets of bioactive compounds [35,36]. Taking this into account, MD simulations and Molecular Mechanics-Poisson-Boltzmann Surface Area (MM-PBSA) calculations were performed for the 60 predicted compound 8-receptor complexes as described in the Materials and Methods section. The detailed results of these calculations are provided in Table 4. The results of the MM-PBSA calculations are summarized in Figure 3.   The energetic and conformational stability of the complexes predicted for each target were analyzed. The plots of the total energy of the systems and of the ligand Root-Mean-Squared Deviation (RMSD) along the 20 production runs performed for each of these complexes are provided as Supplementary Materials in Figures S1 and S2, respectively. The RMSD calculations took place for the ligand, taking as reference its docking predicted pose. Overall, all systems are energetically stable during the production runs. Moreover, the ligand RMSD relative to the reference structure fluctuates around 2 Å or less for all complexes. These RMSD fluctuations are evidence of the conformational changes on the ligand bound to the receptors, indicating that diverse conformations are obtained for MM-PBSA calculations.
The predicted free energies of binding show that the complexes with ACE, TRXE, GSTE2, GLUR and Cdc42 are predicted as unfeasible due to their positive ΔG values. ACE was not predicted by the target fishing approach and there is no evidence that it could be the target of the compounds herein studied. It was included in our investigations because it is the most widely studied target of insecticidal compounds. The fact that the worst (highest) ΔG value among all targets is obtained for ACE is consistent with these observations. Notably, the ACE and GSTE2 receptors that were among the top scored ones according to the molecular docking calculations are discarded as potential targets of compound 8 after the MM-PBSA calculations. These results highlight the importance of evaluating The energetic and conformational stability of the complexes predicted for each target were analyzed. The plots of the total energy of the systems and of the ligand Root-Mean-Squared Deviation (RMSD) along the 20 production runs performed for each of these complexes are provided as Supplementary Materials in Figures S1 and S2, respectively. The RMSD calculations took place for the ligand, taking as reference its docking predicted pose. Overall, all systems are energetically stable during the production runs. Moreover, the ligand RMSD relative to the reference structure fluctuates around 2 Å or less for all complexes. These RMSD fluctuations are evidence of the conformational changes on the ligand bound to the receptors, indicating that diverse conformations are obtained for MM-PBSA calculations.
The predicted free energies of binding show that the complexes with ACE, TRXE, GSTE2, GLUR and Cdc42 are predicted as unfeasible due to their positive ∆G values. ACE was not predicted by the target fishing approach and there is no evidence that it could be the target of the compounds herein studied. It was included in our investigations because it is the most widely studied target of insecticidal compounds. The fact that the worst (highest) ∆G value among all targets is obtained for ACE is consistent with these observations. Notably, the ACE and GSTE2 receptors that were among the top scored ones according to the molecular docking calculations are discarded as potential targets of compound 8 after the MM-PBSA calculations. These results highlight the importance of evaluating molecular complexes derived from docking calculations with more accurate methods such as MM-PBSA for the identification of potential targets for chemical compounds.
Compound 8 is predicted to have very low affinity for a group of proteins that includes LGUL, AKR2, AKR5, RHO1, RAC, GSTD2, RAC1. AKR1, HDAC1 and GSTT1 with predicted ∆G of binding ranging from −0.07 to −2.51 kcal/mol. A third set of targets formed by GSTX2, AKR4, GSTE4, RHO, GSTT2, AKR6, KCC1A and AKR3 is predicted with low to moderate ∆Gs of binding. Interestingly, this set is enriched with proteins annotated with aldo-keto reductase and glutathione transferase functions, both involved in maintaining the redox equilibrium in the cell. Finally, the best values of predicted ∆G of binding are obtained for the HDAC2, CCC2, CA and CCC3 targets.
Overall, the obtained free energies of binding point to HDAC2, CCC2, CA and CCC3 as the most probable targets of the cinnamic acid derivatives herein studied. To get more insights into the possible mechanism of action of these compounds, the predicted complexes of compound 8 with them were analyzed in detail. In Figure 4    A feature common to all complexes depicted in Figure 4 is the extensive network of contacts between the ligand and the receptors. Among the most probable targets of compound 8, two are metalloenzymes (CA and HDAC2) and two are ion channels (CCC2 and CCC3). In both metalloenzymes the carbonyl group of the compound points to the Zn 2+ ions, an interaction that produces favorable electrostatic interactions that positively contribute to the stability of these complexes. The complex of compound 8 with CA is mainly stabilized, in addition to the electrostatic interaction with Zn 2+ , by Van der Waals and hydrophobic interactions. Its phenyl moiety occupies a mostly hydrophobic region lined by W5, Y7, Y66, C67, H97 and T202 while its pentyl tail binds to a zone delimited by H120, V122, V144, L200, T201, V209 and W211.
The predicted binding pose of the ligand to HDAC2 shows the ligand's aromatic ring stacked between H139, F203 and P204 at the entrance of the binding cavity. On the other hand, the aliphatic part of compound 8 orientates to the bottom of the receptor pocket and interacts with M28, G136, L137, G147, C149, G298, G299 and Y301. As in the complex with CA, compound 8 forms no hydrogen bond with the receptor and most interactions are of hydrophobic and Van der Waals types.
The predicted binding modes of compound 8 to the CCC2 and CCC3 ions transporters are very similar. In both cases the phenyl ring of the ligand orientates toward the bottom of the ions binding cavity while its pentyl moiety heads to the entrance of the pocket. The root mean square deviation (RMSD) between the conformations adopted by compound 8 in both receptors is 2.59 Å. The highly similar binding modes of the ligand to these ion channels can be a consequence of the 65% overall identity between their transmembrane domains and, specifically, of the 79% identity among the residues directly interacting with compound 8.
The regions of the binding cavities accommodating the phenyl substituent of compound 8 are formed by I157, G159, V160, M161, N241, P361, S362 and T465 in CCC2 and by the corresponding I137, G139, V140, M141, M220, N221, P341, S342 and T447 residues in CCC3. On the other side of the cavity, the aliphatic moiety of the ligand interacts with F230, N234, T364, G365, Q367, A368, N469 and Y533; and with N136, F210, N214, N217, T344, S348, N451 and F515 in CCC2 and CCC3, respectively. In addition, compound 8 is predicted to hydrogen bond the side chains of N469 and Q237 of CCC2 through its carbonyl group and the main chain oxygen, respectively. The first of these hydrogen bonds is observed in more than 97% of the analyzed MD snapshots while the second one is less frequent (22% of the analyzed snapshots). In contrast, the network of hydrogen bonds predicted for CCC3 is more extensive and it includes the side chains of N217, N451 and N221. All these three residues can serve as hydrogen bond donors for the carbonyl group of compound 8, while the main chain oxygen of the ligand can accept a hydrogen bond from N217. Overall N217, N451 and N221 of CCC3 are observed to hydrogen bond the ligand in 82%, 73% and 42% of the analyzed MD snapshots, respectively.
The slightly better ∆G of binding predicted for CCC3 compared to CCC2, despite the high identity of their binding pockets, can be explained from the different networks of hydrogen bonds that they form with compound 8. This difference arises from the N217 (CCC3) to Q237 (CCC2) mutation. In both receptors several hydrogen bonds are formed between the residues pointing to the binding cavity. One of these is formed between the side chains of N241 and Q237 of CCC2, mediated by their amine and carbonyl groups, respectively. The same hydrogen bond is observed between N221 and N217 of CCC3. Notably, these hydrogen bonds are observed during all the 200 MD snapshots analyzed in each system and are not disrupted by the ligand. The side chain of glutamine is one carbon longer than that of asparagine, making the amine and carboxyl groups of the glutamine side chain locating furthest from the protein backbone than in asparagine. This leads to the movement of the side chain of N241 in CCC2, pulled by the hydrogen bond with Q237, to a position unfavorable for hydrogen bonding the carbonyl group of compound 8. This is not the case of CCC3, where the hydrogen bond between the side chains of the corresponding N221 and N217 favors the formation of hydrogen bonds between them and the ligand.
All these results suggest a multi-target mechanism of action of the cinnamic acid derivatives herein studied. This hypothesis is supported by literature reports where cinnamic acid has been evaluated for its larvicidal activity against Ae. aegypti. For instance, the low and moderate ∆Gs of binding predicted for GSTD2, GSTT1, GSTX2, GSTE4 and GSTT2 agree with the inhibitory activity of the glutathione s-transferase activity produced by cinnamic acid results reported by Nobsathian and collaborators (2018), [31].
On the light of our results, we postulate that the principal targets of the cinnamic acid derivatives are CA, HDAC2, CCC2 and CCC3.
Histone deacetylase activity is related to maintaining the stability of the chromosomes in the cell and in proteins expression. However, we found no study in the scientific literature focusing on the larvicidal effect of histone deacetylase inhibitors. On the other hand, carbonic anhydrases has been explored as the targets of Ae. aegypti larvicidal compounds [37]. Specifically, CA has been proposed to play an important role in bicarbonate metabolism, pH regulation and ions transportation in the alimentary canal of Ae. aegypti [38]. Thus, its inhibition is a strategy to explore for the development of larvicidal compounds.
CCC2 and CCC3 have been characterized as Na + -coupled cation chloride cotransporters (CCCs) and associated to ions transport [39,40]. However, the exact mechanisms by which this process occurs remains unclear. It will require more investigations to test the hypothesis that the studied cinnamic acid derivatives could exert their larvicidal effect through the blockage of CCC2 and CCC3. However, we consider that this particular mechanism of action deserves additional attention in future research efforts.

Chemical Characterization and Reagents
All of the chemical products used during synthesis were from Sigma-Aldrich. 1 H-NMR (200 and 50 MHz), (400 and 100 MHz) and 13 C-NMR (500 and 125 MHz) spectra were respectively recorded on VARIAN MERCURY, BRUKER-ASCEND and VARIAN-RMN-SYSTEM spectrometers. Chemical shifts (δ) are expressed in parts per million (ppm) using TMS as an internal standard. Spin multiplicities are given as s (singlet), brs (broad singlet), d (doublet), t (triplet), q (quartet), quint (quintet), sex (sextet), sept (septet) and m (multiplet). High Resolution Mass Spectrometry was performed using an Ultraflex II TOF/TOF spectrometer with a high performance solid state laser (λ = 355 nm), and a reflector using the MALDI technique-matrix assisted laser desorption ionization. Column adsorption chromatography (CC) was performed on silica gel (Merck 60, 230-400 mesh); analytical TLC was performed on pre-coated silica gel plates (Merck 60 F 254 ). FTIR spectra were recorded in a Bruker FTIR spectrometer, Vertex 70 model, using KBr pellets.

Synthesis of Compounds 2-9; Fischer's Esterification
To a 100 mL flask, cinnamic acid (0.25 g, 1.69 mmol) and alcohol (50 mL) were added in the presence of sulfuric acid (0.4 mL), and were heated under reflux until completion of the reaction (5-24 h), being verified by single spot TLC. The alcohol was then removed under reduced pressure until about half and the solution was diluted with 20 mL of water. The product was extracted using ethyl acetate (3 × 15 mL), and the organic phase was neutralized successively with 5% sodium bicarbonate, washed with water, then dried over anhydrous sodium sulfate, and filtered. After evaporation under reduced pressure, certain compounds required a column for purification, this phase yielded the ester derivatives, shown in Scheme 1 [41].

Synthesis of Compound 10; Reaction with Halide
A mixture of cinnamic acid (0.2 g, 1.35 mmol) in acetone (16.4 mL) was heated under reflux in the presence of triethylamine (0.73 mL) and halide (1.39 mmol) until a complete reaction (24 h), which was verified by single spot TLC. The solvent was then removed under reduced pressure and the solution was diluted with 20 mL of water. The product was extracted with ethyl acetate (3 × 20 mL). The organic phase was treated with water, dried over anhydrous sodium sulfate, filtered, and concentrated under reduced pressure. The residue was purified on a silica gel column (eluent: hexane-ethyl acetate, 9:1), as shown in Scheme 1 [42].

Synthesis of Compounds 11-18; the Mitsunobu Reaction
In a 10 mL flash we added cinnamic acid (0.1 g, 0.68 mmol) and alcohol (0.68 mmol); being dissolved in 2.25 mL tetrahydrofuran. The reaction mixture was stirred under magnetic stirring at 0 • C for about 30 min. Afterwards, diisopropyl azodicarboxylate (0.11 mL, 0.68 mmol) and triphenylphosphine (0.15 g, 0.68 mmol) were added as esterification agents, with continuous stirring at room temperature for about 44-68 h, which was verified by a single spot in TLC. The solvent was removed under reduced pressure and the solution was diluted with 10 mL of water. The product was then extracted with ethyl acetate (3 × 10 mL), and the organic phase was neutralized successively with saturated sodium bicarbonate, washed with brine, dried over anhydrous sodium sulfate, filtered, and evaporated under reduced pressure. The residue was then purified on a silica gel column (eluent: hexane-ethyl acetate, 9:1), as shown in Scheme 1 [43].
The compounds were diluted in distilled water (10 mL) at different concentrations (0.0125 to 10 mg/mL). Some substances were diluted in 4% DMSO, while others in 4% Tween80, (due to the compound's structural characteristics), for later solubilization in distilled water. For a dose response curve and LC 50 determination, fifteen L4-stage larvae from Ae. aegypti were transferred into Falcon tubes containing the 10 mL solution of compounds at different concentrations. The negative control group was prepared using distilled water and a solubilizing agent. The positive control group was prepared using a solution of the insecticides Imiprothrin 0.02%, Permethrin 0.05% and Esbiothrin 0.1% [3]. The tubes were incubated for 24 h at 27 ± 2 • C, over 12 h of natural light and 12 h of darkness [53,54]. Larvae mortality was verified after 24 h of incubation and all tests were carried out in triplicate.

Larvae Mortality Rate
The larval mortality rate was verified 24 h after exposure to the compounds. The larvae were transferred to petri dishes, being considered lifeless when a total absence of movement was observed. Larvicidal activity was calculated according to Equation (1).

Statistical Analysis
Statistical analysis was performed using GraphPad Prism software for Windows version 5.0 (GraphPad Software, San Diego, CA, USA). Significant differences between groups were analyzed by ANOVA followed by the Tukey post-test when appropriate (p < 0.05).

Molecular Modeling
The overall molecular modeling workflow used in this research bases on that employed in our previous publications [35,36]. The potential molecular targets of the studied series of compounds in Ae. aegypti was identified through a consensus computational target fishing approach. Then, the most active compound in the series (8) was docked to its potential targets to obtain potential ligand-receptor complexes. The free energies of binding of the compound to its potential targets were estimated from molecular dynamics simulations. Then, the most probable targets of the studied chemicals were selected as those providing the lowest free energies of binding. Figures containing molecular structures were obtained with UCSF Chimera [55], the frequencies of interactions were visualized with Cytoscape [56] and the interaction diagrams were obtained with LigPlot+ 2.2 [57] from the representative MD snapshot. Otherwise noted, the parameters for modeling were kept to their default values in all used software. All molecular modeling steps are described in detail below.

Targets Selection
Compounds 2, 3, 5, 7, 8 and 13 were independently used as input to the Similarity Ensemble Approach (SEA) web server [58] for the identification of their potential molecular targets. From the lists of predicted targets for these compounds, those predicted for all the six chemicals were selected as potential targets of the congeneric series of compounds. The consensus list of potential targets was used as input to a Blast [59], search, through its NCBI web implementation, against the Ae. aegypti (taxid:7159) proteome on the Reference proteins (refseq_protein) database. Ae. aegypti proteins identical in at least 40% to any query sequence and covered in at least 80% by the Blast alignment were selected as potential targets of the studied compounds in this organism. Based on literature reports, six glutathione s-transferase and the acetylcholinesterase enzymes of Ae. aegypti were also studied as potential targets of the assayed compounds [31,32].

Molecular Docking
The initial three-dimensional structure of compound 8 was obtained with OpenEye's Omega [60] and partial atomic charges were added to it with Molcharge [61]. For molecular docking, the structure of GSTE2 (code 5FT3) was obtained from the Protein Data Bank [62]. The rest of the evaluated Ae. aegypti proteins had no three-dimensional structure solved and homology models were built for them using the SWISS-MODEL web server [63]. Any catalytically relevant cofactor was added to the models of the receptors prior to molecular docking. Furthermore, the binding of compound 8 was evaluated on the glutamate binding site of the GLUR receptor only.
Docking calculations proceeded with the Gold program [64] as previously described [35,36]. Briefly, compound 8 was docked to the selected receptors using the CHEMPL scoring function for primary docking. The search efficiency of Gold was set to 200% and 30 different ligand conformations were predicted for each target. The obtained ligand poses were rescored with the ASP, ChemScore and GoldScore scoring functions of Gold. The obtained binding modes were ranked in descending order according to the consensus score defined by Equation (2): where, S i,j is the score of pose i according to scoring function j and S j and std(S j ) represent the average of scoring function j and its standard deviation, respectively. During molecular docking, the side chains of the amino acids directly pointing to the binding cavity were considered flexible. Any predicted pose of compound 8 with a consensus score higher than 1 was selected for further analyses.

Molecular Dynamics Simulations and Free Energies of Binding Estimation
Molecular dynamics (MD) simulations and the estimation of the free energies of binding were carried out with Amber 2018 [65]. The Amber forcefields ff14SB, lipid17 and gaff were used for proteins, lipids and non-aminoacidic residues, respectively. All MD simulations proceeded in explicit solvent.
Compound 8-receptor complexes were divided into two classes for MD simulations. The first group consisted in water-soluble proteins and the second one was formed by membrane receptors. The later class included the CCC2, CCC3 and KCC1A proteins while the rest of the targets formed the first group. The complexes containing soluble proteins were parameterized with the tleap tool of Amber. These systems were enclosed in truncated octahedron boxes, solvated with TIP3P water molecules and any excess charges in them were neutralized by the addition of either Cl − or Na + ions. Long range electrostatic interactions were treated with the Particle Mesh Ewald (PME) method in all MD simulations steps. The solvated systems underwent a two steps energy minimization protocol, the first of which consisted in 500 steps of the steepest descent method followed by 500 cycles of conjugate gradient at constant volume. During this stage, all atoms except the solvent and counterions were restrained with a force constant of 500 kcal/mol. Å 2 and the PME distance cutoff was set to 12 Å.
The second part of energy minimization took place at constant volume too and included 1500 iterations of the steepest descent method and 1000 cycles of conjugate gradient, with no restraints and a PME distance cutoff of 10 Å. Afterward, the energy minimized systems were gradually heated from 0 K to 300 K for 20 ps with a time step of 2 fs, keeping all atoms except solvent and counterions restrained with a force constant of 10 kcal/mol.Å 2 . From this step on, the PME distance cutoff was set to 10 Å, all bonds involving hydrogen atoms were constrained with the SHAKE algorithm and temperature was controlled using a Langevin thermostat with a collision frequency of 1.0 ps −1 .
Next, the systems were equilibrated for 100 ps using a time step of 2 fs at constant pressure and constant temperature. At this step pressure was set to 1 bar and controlled with isotropic position scaling setting the relaxation time to 2 ps. Temperature was set to 300 K. The last snapshot of the equilibration process was used as input to 20 MD production runs. Each of these lasted 2 ns and were initialized with velocities to allow for a better sampling of the complex conformational space.
The complexes of compound 8 with the membrane proteins were prepared with the CHARMM-GUI server [66,67]. These complexes were integrated in a bilayer membrane containing 50 1-palmitoyl−2-oleoylphosphatidylcholine (POPC), 50 1-palmitoyl−2oleoylphosphatidylethanolamine (POPE) and 25 cholesterol (CHL) molecules on each side. The systems were enclosed in rectangular boxes that were later solvated with TIP3P water molecules and neutralized by adding K + and Cl − counterions. The solvated systems were then energy minimized, heated and equilibrated in Amber using the input files provided by CHARMM-GUI. Production runs were run as for the soluble proteins.
The free energies of binding of compound 8 to all receptors were estimated with the MM-PBSA method as implemented in Amber [67]. MD snapshots for MM-PBSA calculations were extracted from the 40 ns simulations of each target every 200 ps. This lead to 200 snapshots being selected for the free energy of binding calculation on each complex. The ionic strength for these calculations was set to 0.100 mM. Membrane proteins were treated under the implicit membrane model of MMPBSA.py with the membrane dielectric constant set to 4.

Conclusions
The present study investigated the larvicidal activity of seventeen cinnamic acid esters derivatives against Ae. aegypti L. L4 larvae. Based on the findings of the investigation, we concluded that the butyl cinnamate (7), pentyl cinnamate (8) and benzyl cinnamate (13) derivatives present the highest larvicidal activities. The study revealed that certain structural features are important for larvicidal action, such as: alkyl substituents with medium chain; or aryl groups without substituents. Extensive molecular modeling studies on 60 possible protein-inhibitor complexes that involved a total of 2.4 µs MD simulations lead to the proposal of a multi-target mechanism of action for these compounds. Further, the modeling results may guide future experiments focusing on optimizing the larvicidal activity of cinnamic acid derivatives against Ae. aegypti, while elucidation their mechanism of action, involving CA, HDAC2, CCC2 and CCC3 proteins, and motivating additional in vitro studies to optimize the larvicidal activity of this class of compounds.
Supplementary Materials: The following are available online, Figure S1: Total energy, in kcal/mol, of the studied complexes along the MD simulations, Figure S2: Ligand RMSD, in Å, relative to the starting docking conformation along the MD simulations.