From the ChEMBL database, we select a set of chemical structures for the construction of a predictive model. The set was composed by 162 diverse chemical structures that have been studied and inhibited larvae of Ae. aegypti.

The compounds were classified using values of −logIC₅₀ (mol/L) = pIC_50, which led us to assign 85 active (pIC₅₀ ≥ 4.15) and 76 inactive compounds (pIC₅₀ < 4.15). In this case, the IC₅₀ represented the concentration required for 50% inhibition of larvae of Aedes aegypti. Through a bibliographical search, 11 flavonoids (Gikonyo et al., 1998; Rajkumar and Jebanesan, 2008; Perumalsamy et al., 2015; Guarda et al., 2016) have been added to this bank against Aedes aegypti larvae, with a total of 173 molecules. Another dataset of molecules isolated from S. paniculatum was built from the phytochemical study of this species. For all structures, SMILES codes were used as input data in a Marvin 18.10, 2018, ChemAxon (http://www.chemaxon.com) (ChemAxon, 2018). We used (Standardizer software 2018) [JChem 18.10, 2018; ChemAxon (http://www.chemaxon.com)] to canonize structures, add hydrogens, perform aromatic form conversions, clean the molecular graph in three dimensions, and save compounds in sdf. format. Three dimensions structures (3D) were used as input data in the software Dragon 7.0 to generate descriptors (Talete, 2016) where it is possible to predict the biological and physical-chemical properties of the molecules. This calculus was realized for the set of chemical structures with knowledge activity against Ae. aegypti larvae.

The Dragon 7.0 software (Talete, 2016) generated 1,232 descriptors for 173 molecules of known activity against Ae. aegypti larvae. These descriptors were used as input date no Knime software for generation the predictive model. The Dragon 7.0 software can calculate 5,270 molecular descriptors, covering several approaches. These descriptors are organized in 30 logical blocks. The descriptors list includes the simplest types of atoms, functional groups and fragment counting, topological and geometric descriptors, tridimensional descriptors, but also several of them estimates properties such as logP, and Lipinski (Mauri et al., 2006).

The Knime 3.4.0 software (Knime 3.4.0 the Konstanz Information Miner Copyright, 2003–2014, www.knime.org) (Berthold et al., 2006) was used to perform the analyzes and created the model in silico. The descriptors and class variables were imported from the software Dragon 7.0, and for each one the data were divided using the “Partitioning” node with the “stratified sample” option to create a training set and a test set, and encompassing 80 and 20% of the compounds, respectively.

Although the compounds were selected randomly, the same proportion of active and inactive samples was maintained in set. For internal validation, we employed cross-validation using 10 randomly selected, stratified groups, and the distributions according to activity class variables were found to be maintained in all validation groups and in the training set.

Descriptors were selected, and a model was generated using the training set and the Random Forest algorithm (RF) (Salzberg, 1994), using the WEKA nodes (Hall et al., 2009). The parameters selected for RF included the following settings: number of trees to build = 50, seed for random number generator = 1, for the model. The internal and external performances of the selected models were analyzed for sensitivity (true positive rate, i.e., active rate), specificity (true negative rate, i.e., inactive rate), and accuracy (overall predictability).

Additionally, the sensitivity and specificity of the Receiver Operating Characteristic (ROC) curve was find to describe true performance with more clarity than accuracy. It was also used the domain of applicability based on Euclidean distances in order to signal compounds in the test set for which the predictions may be unreliable. Similarity measurements are used to define the applicability domain of the model based on Euclidean distances between all training, test and virtual screening compounds.

The distance from a compound of a test compound to its nearest neighbor in the training set is compared to the predefined applicability domain limit, if the similarity is beyond this limit, the prediction is considered unreliable (Zhang et al., 2006).

The target proteins of Ae. aegypti 1YIY (Han et al., 2005), 1PZ4 (Dyer et al., 2003), with their respective inhibitor ligands were downloaded from Protein Data Bank (http://www.rcsb.org/pdb/home/home.do), details of each enzyme can be seen in Table 1.

The sulfated flavonoids isolated from S. paniculatum were submitted to molecular docking using the Molegro Virtual Docker, v. 6.0.1 (MVD). All the water compounds were deleted from the enzyme structure, and the enzyme and compound structures were prepared using the default parameter settings in the software package (ligand evaluation: Internal ES, Internal H-Bond, Sp2-Sp2 Torsions, all checked; number of runs: 10 runs; algorithm: MolDock SE; maximum interactions: 1,500; Max. population size: 50; Max. steps: 300; neighbor distance factor: 1.00; Max. number of poses returned: 5).

The docking procedure was performed using a GRID of 15 Å in radius and 0.30 in resolution to cover the ligand-biding site of the enzyme’s structures. Template docking was used to focus the search on poses similar to the ligand interactions and conformation. The Moldock score algorithm was used as the score function (Thomsen and Christensen, 2006). The energy of the crystallized ligand of each protein is calculated automatically by the Moldock Score. Based on the pose of the crystallized ligand, the Moldock score algorithm automatically converts the ligand’s energy to the energy scale used by it.

Molecular dynamics simulations were performed to estimate the stability of interactions between proteins and ligands using Gromacs 5.0 software (Berendsen et al., 1995; Abraham et al., 2015). The topology of the ligands was prepared using the PRODRG topology generator (http://davapc1.bioch.dundee.ac.uk/cgi-bin/prodrg/submit.html) (Schuttelkopf and Aalten, 2004) applying the GROMOS43a1 force field. The protein topology was also prepared using the GROMOS43a1 force field. The molecular dynamics simulation was performed using the SPC water model of point load extended in a dodecahedral box (Bondi, 1964). The system was neutralized by adding ions (Cl⁻ and Na⁺) and minimized to remove bad contacts between complex molecules and solvent. The system was also balanced at 300 K using the V-rescale algorithm at 100 ps as NVT set (constant number of particles, volume, and temperature) followed by equilibrium at 1 atm of pressure using the Parrinello-Rahman algorithm as set NPT (constant number particles, pressure and temperature) to 100 ps. The DM simulations were performed in 50,000 cycles at 1 ns. In order to determine the stability of the structure and if the complex is stable near the experimental framework, the mean square root displacement (RMSD) of all heavy atoms was calculated relative to the starting structures. Residual fluctuations (RMSF) were also analyzed to understand the role of residues near the receptor binding site. To verify the interaction energy of the protein-ligand complex, we calculated the short-range Coulombic interaction energy and the Lennard-Jones short-range energy for 10 ns. The RMSD and RMSF plots were generated in the Grace (http://plasma-gate.weizmann.ac.il/Grace/) software and the protein and ligands were visualized at UCSF Chimera (Pettersen et al., 2004).

Analyzing the Ae. aegypti model, the internal cross validation and the external test demonstrated similar statistical performance, with accuracy higher than 81%, and showing to be a model with great performance. The Table 3 summarizes the statistical rates of the RF model.

With these results it was possible to calculate the Matthews correlation coefficient (MCC) for generation evaluation of the Ae. aegypti model. The MCC correlations observed and predictive binary classifications, resulted in a value between −1 and +1, where +1 is a perfect prediction, and −1 indicates a complete disagreement between prediction and observation. For the external test the value obtained was 0.70 and 0.72 for the internal cross validation, revealing a good prediction, and robustness of the model.

The graphic of ROC curve that analyze the model performance was generated for the external test set and obtained an area under the curve of 0.93, and 0.91 for the internal cross validation (Figure 2). Considering that a perfect model has an area under the curve equal to 1, it is possible to state that the model can perform a high classification rate for this RF method.

The Ae. aegypti model was used in order to triage the sulfates flavonoids isolated from S. paniculatum for the identification of possible bioactive molecules against Ae. aegypti larvae. The six flavonoids tested showed a potential activity of more than 53% (Figure 3).

It is possible to observe in Figure 3 that the position of the sulfate group influences the probability of potentially active of the flavonoids. When the sulfate group is in position 7, the probability of active potential is greater than when it is positioned in position 8. Another fact that we can observe is that the substitution of a methoxyl in position 4′ also positively influences the probability of active potential flavonoids, a fact best observed in flavonoids 9a and 9b.

The six sulfates flavonoids isolated from S. paniculatum were tested with the two Ae. aegypti target proteins (1YIY and 1PZ4) in order to screen the compounds with the highest potential larvicidal activity, thus a virtual screening based on the receptor structure. Then, calculations were made to obtain the molecules with the best probabilities of active potential against Ae. aegypti larvae. These calculations were carried out using the following formula: Prob= E MT E M , if  E MT < E L Where, E_MT is the energy of the test molecule, E_M is the lowest energy obtained from the molecules tested and E_L is the energy of the ligand obtained from the crystallography of the test protein, the fraction is conditioned to the energy of the test molecule being smaller than the energy of the ligand, that is, and only molecules which have obtained binding energy lower than the binding energy of the crystallographic binder will be considered potentially active (Acevedo et al., 2018; Barros et al., 2019).

Molecular docking was validated using the redocking technique, in which the inhibitor ligand of each protein obtained the same pose after the docking was performed, and thus validating the docking performed in the study. The energy offered by Moldock is not determinant to influence the validation of the docking, because in the literature the methodology used was another one, and the crystallography.

Flavonoid 6 obtained the best probability of interaction at the active site of the 1PZ4 protein but did not obtain a probability of interaction at the active site of the 1YIY protein. All other flavonoids were likely to fit and interact at the active site of the two studied proteins. These results can be seen in Table 4.

Calculations were also made to select multitarget molecules through the docking, combined with the virtual screening model, that is, to select molecules that are potentially active for the two enzymes, and which were also predicted to be active in the created virtual screening model. For these calculations, the following formula was used: Pro b Comb =  (Prob+ ( 1+ Es p Cross ) x  P RF 2+ Es p Cross  , If Pro b Comb  >0,5 Where, Prob is the potential probability, Esp_Cross is the value of the specificity in the cross-validation of the virtual screening model, and P_RF is the probability of each flavonoids being active in the RF model. And this combined probability is conditioned, that is, only molecules with values above 0.5 will be considered potentially active (Acevedo et al., 2018; Barros et al., 2019). Already eliminating also, the two molecules that had results considered not reliable by the model. The results are presented in Table 4.

Only flavonoid 6 was not multitarget since it was not likely to interact at the active site of the 1YIY protein. Also analyzing the results of the combined probability between the molecular docking and the RF model, we can see that the flavonoids that have the sulfate group in position 7 were shown to be more potentially active in both the model and in the 1PZ4 protein, a protein that is present in the larvae intestine. Ae. aegypti mosquito suggests that the sulfate in this position favors the fitting of the flavonoid in the active site of the protein.

While flavonoids that have the sulfate group in position 8 were more likely to have an active potential in 1YIY, a protein present in the head of the Ae. aegypti adult mosquito, and thus it is suggested that in this position the flavonoid presents a better fit into the active site of the protein.

To characterize the interactions and to evaluate the interaction energies between the proteins and the ligands, simulations of molecular dynamics (MD) were performed. The analysis of the mean square deviation (RMSD) of the sterol protein (PDB ID 1PZ4) showed that the enzyme reached conformations ranging from 1,5 to 2,4 nm in size in 1 ns. While Kinurenine aminotransferase (PDB ID 1YIY) protein showed a RMSD with conformations ranging from 1.5 to 3.5 nm in 1 ns. The results show that the protein Kinurenine aminotransferase is more unstable, requiring greater mobility and interaction with its ligands (Figures 4, 5).

The mean square root fluctuations (RMSF) of the individual residues were also calculated. It was verified that residues 1, 23, 25, 105, 106, and 110 of the Sterol -sulfated flavonoid complexes (Mol 6, 7, 8, 9a, 9b, and 10) presented higher fluctuation values, indicating the key residues in the conformational change of the protein structure (Table 5). Only the first residue of the Sterol-PLM complex (PDB binder) showed a high fluctuation value. Kinurenin aminotransferase complexed to flavonoids also showed several residues with high RMSF values (above 0.3); corresponding to residues 12, 13, 14, 15, 16, 16, 17, 18, 22, 64, 65, 75, 159, 160, 262, 308, 357, 423, 427, 428, and 429. While complexed Kinurenine to the PDB binder (PMP) showed few high fluctuations, occurring in residues 12, 13, 14, 73, 308, and 428. This indicates that flavonoids contribute more to the flexibility of the protein, being able to interact with a greater number of amino acids.

The Coulombic and Lennard-Jones interaction energies of the protein-ligand complexes were calculated and the flavonoid linkers (6, 7, 8, and 9b) were shown to have a higher interaction stability with the sterol enzyme when compared to the ligand PDB PLM (Figure 6). As for Kinurenin, all flavonoid binders showed greater interaction stability of Lennard-Jones than PDB PMP ligand (Table 5).

Although the RMSD analysis of the sterol-PLM complex is better when compared to the flavonoids, we can suggest that through RMSF and interaction energy calculations, flavonoids interact better with the two proteins because they allow greater flexibility and interaction with amino acids, and interaction stability.