Anopheles gambiae Trehalase Inhibitors for Malaria Vector Control: A Molecular Docking and Molecular Dynamics Study

Simple Summary Anopheles gambiae is a major malaria vector. This vector is controlled mainly using insecticide, thereby reducing malaria transmission. Since many insecticides are toxic to non-target species, and there is increasing insecticide resistance, there is a need to identify other safer molecules that can be used to control the vector. Molecules inhibiting trehalase have been proposed as safer options for insecticide and fungicide development since humans do not produce trehalose. This work screens several trehalase inhibitors against trehalase of Anopheles gambiae (AgTre) and assesses their safety in humans using in silico methods. Four trehalase inhibitors had a high affinity for AgTre and were predicted to be relatively safe. The compounds also interacted well with trehalase of Aedes aegypti, another important mosquito vector. The results suggest these molecules are safer options for insecticide development to control mosquito-borne diseases. Abstract Trehalase inhibitors are considered safe alternatives for insecticides and fungicides. However, there are no studies testing these compounds on Anopheles gambiae, a major vector of human malaria. This study predicted the three-dimensional structure of Anopheles gambiae trehalase (AgTre) and identified potential inhibitors using molecular docking and molecular dynamics methods. Robetta server, C-I-TASSER, and I-TASSER were used to predict the protein structure, while the structural assessment was carried out using SWISS-MODEL, ERRAT, and VERIFY3D. Molecular docking and screening of 3022 compounds was carried out using AutoDock Vina in PyRx, and MD simulation was carried out using NAMD. The Robetta model outperformed all other models and was used for docking and simulation studies. After a post-screening analysis and ADMET studies, uniflorine, 67837201, 10406567, and Compound 2 were considered the best hits with binding energies of −6.9, −8.9, −9, and −8.4 kcal/mol, respectively, better than validamycin A standard (−5.4 kcal/mol). These four compounds were predicted to have no eco-toxicity, Brenk, or PAINS alerts. Similarly, they were predicted to be non-mutagenic, carcinogenic, or hepatoxic. 67837201, 10406567, and Compound 2 showed excellent stability during simulation. The study highlights uniflorine, 67837201, 10406567, and Compound 2 as good inhibitors of AgTre and possible compounds for malaria vector control.


Introduction
Mosquitoes such as Anopheles and Aedes spp. are important vectors of several human diseases [1]. These mosquito-borne diseases, such as malaria, are controlled principally by employing vector-control strategies [2]. These strategies aim to reduce or eliminate the contact of the vector with humans, thus preventing parasite transmission. The primary vector control strategy for malaria involves the use of insecticides either in insecticide-treated nets (ITNs) or for indoor residual spraying (IRS) [3]. These strategies have contributed considerably to averting clinical cases of the disease, with about sixty-eight per cent of the 663 million averted cases between 2000 and 2015 attributed to the use of ITNs [4]. This highlights the importance of vector control in curbing malaria. However, increasing reports of insecticide resistance, the toxicity of insecticides to non-target species, and the persistence of insecticides in the environment have necessitated the search for safer alternatives in vector control [5,6]. Larval source management (LSM) strategies, which prevent the emergence of adult mosquitoes, are one of the complementary strategies employed in vector control [7]. These methods reduce the number of house-entering as well as outdoor biting mosquitoes and have been observed to reduce malaria transmission when combined with the use of ITNs [8]. Amidst this, increasing research efforts are being made to identify and develop new classes of safer compounds that could be used in vector control.
Trehalase inhibitors have been proposed as safer alternatives for insecticides as they are non-toxic and do not persist in the environment [9]. Trehalase (EC 3.2.1.28) is an enzyme that catalyses the reduction of trehalose into two molecules of glucose [10,11]. This reaction is critical in insects, as trehalose is the major hemolymph sugar in insects and is also the fuel required for flight [12]. Trehalose is also needed in larvae to resist stress factors [13]. Trehalase levels have been reported to be upregulated in the salivary gland of Anopheles stephensi during Plasmodium vivax infection, suggesting that it might also play a crucial role in parasite development in the mosquito [14]. Similarly, inhibition of trehalose transporter in Anopheles gambiae was observed to lead to decreased survival of the mosquito as well as decreased Plasmodium development in the mosquito [15]. Inhibition of trehalase by validamycin has been proposed as promising for the development of new insecticides [16]. Validamycin A has been tested in different insect species. For example, it has been observed to cause developmental defects and decrease fecundity in Helicoverpa armigera [17] and cause mortality in Diaphorina citri [18]. Marten et al. [19] reported that validamycin A, a trehalase inhibitor, reduced egg hatching and pupation, skewed the sex ratio, and prevented flight in Aedes aegypti, the vector for dengue fever, chikungunya, Zika fever, Mayaro, and yellow fever. Minimal larva death in Ae. aegypti (7% mortality at 100 and 200 ppm, and 9% mortality at 500 ppm validamycin A administration) was observed in their study. An earlier report by Logan [20] revealed that validoxylamine A, another trehalase inhibitor, did not cause significant mortality of Ae. Aegypti larvae but prevented flight of emerging adults and decreased the number of eggs laid by adult Ae. aegypti mosquitoes. This suggests that trehalase inhibitors can serve as insect growth regulators, sterilants, and flight inhibitors. This is important as the flight of disease vectors to access humans is required for disease transmission. Similarly, sterilants can help to reduce mosquito populations. There is growing research on synthesising other potential trehalase inhibitors and testing their activities in different species as fungicides or insecticides [9]. Although these potential inhibitors have been tested in different organisms, no study has tested these compounds in the malaria vector, Anopheles spp., to elucidate their possible use for malaria vector control. While validamycin A and validoxylamine A have been tested in Ae. aegypti, validamycin A is considered a weaker binder of trehalase, and the safety of validoxylamine A is of concern, since it is a good binder of porcine trehalase. Hence, there is a need to identify more potent and safe inhibitors of trehalase aside from validamycin A in these disease vectors, which can be developed for vector control. These can be carried out using traditional drug design methods or computer-aided drug design methods (CADD).
Traditional drug design methods are costly as they involve the synthesis of numerous compounds followed by in vitro or in vivo testing of their activities. For example, Insects 2022, 13, 1070 3 of 21 D'Adamio et al. [13] synthesised more than 20 compounds and tested their activity against Chironomus riparius trehalase. Apart from the cost involved in synthesising numerous compounds, in vitro or in vivo testing of compounds depend on the availability of recombinant enzymes or the organisms, which come with a cost that multiplies with an increasing number of compounds. However, CADD methods, which include molecular docking and molecular dynamics (MD) simulation, offer cheaper and faster means of screening and identifying potential binders of a protein from a wide library of compounds [21]. These methods rely on the availability of the three-dimensional (3D) structure of the protein of interest and a compound library [22]. Protein structure can be experimentally determined or computationally derived [23]. Potential hits, obtained based on docking scores from molecular docking, can be optimised and re-screened, leading to the guided synthesis of compounds. This helps to reduce the costs involved in insecticide/drug development. Performing molecular dynamics simulation after molecular docking helps to reveal the stability of protein-ligand interactions and may separate false positive from true positive hits [21]. As described above, trehalase inhibitors are promising vector-control compounds. Unfortunately, there is no report about their use against An. gambiae, an important malaria vector. In this study, the 3D-structure of An. gambiae trehalase (AgTre) was determined and molecular docking was employed to screen a library of compounds in order to identify potential lead structures against AgTre that could be developed for malaria vector control. Considerations were given to the interaction of these compounds to trehalase of Ae. aegypti, which is also an important mosquito vector of human diseases. Since trehalase from both organisms share >70% sequence identity, inhibitors could be used to control both vectors, preventing transmission of other mosquito-borne diseases aside from malaria, e.g., dengue fever. In addition, the medicinal chemistry and toxicity potential of these compounds were predicted, and the stability of their complexes with trehalase was determined. Four trehalase inhibitors were observed to be better binders of trehalase of An. gambiae and Ae. aegypti than validamycin A. The compounds were predicted to be non-toxic when compared to validoxylamine. Hence, the study proposes these inhibitors as potential molecules for vector control.

Prediction of 3D Structure of AgTre
Most insects have at least two copies of trehalase genes; however, it has been reported that dipteran insects, such as Anopheles spp., have only one copy of trehalase [24]. The single copy gene of trehalase in An. gambiae is AGAP012053, which codes for its trehalase protein AgTre. The experimental 3D structure of AgTre has not been determined; hence it is not available in the database of the Research Collaboratory for Structural Bioinformatics Protein Data Bank (RCSB PDB) [25]. However, a predicted structure based on homology modelling is available in the SWISS-MODEL Repository [26]. The structure was built using trehalase of Arabidopsis thaliana (7E9X_A) as a template which has the highest percentage identity (33%) to AgTre. Considering the low percentage identity and family distance of the two proteins (plant to insect), ab initio modelling was implemented in predicting the 3D structure of AgTre. The FastA file containing the protein sequence of the AgTre (Trehalase from An. gambiae Pest strain) was retrieved from UniProt Knowledgebase (UniProtKB) [27] with the accession number Q7PZS4. The protein sequence of AgTre was submitted to three ab initio modelling web servers: Iterative Threading Assembly Refinement (I-TASSER), Contact-guided Iterative Threading ASSEmbly Refinement (C-I-TASSER) [28,29], and the Robetta web server. On the Robetta server, a deep learning-based modelling method, RoseTTAFold [30], was used to predict the 3D structure of AgTre. RoseTTAFold gives confidence scores for models (ranging from 0 to 1), with models closer to 1 being better structures. This confidence score corresponds to predicted IDDT score using DeepAccNet [31]. For I-TASSER and C-I-TASSER, default parameters were used in predicting the 3D structure for AgTre. I-TASSER and C-I-TASSER output C-scores within the range between −5 and 2, with a higher C-score signifying higher model confidence.

Structure Assessment of AgTre Models
Structure assessment of all predicted models in this study (RoseTTAFold, I-Tasser, and C-I-Tasser), as well as available models from SWISS-MODEL and AlphaFold, was carried out using the SWISS-MODEL web server structure assessment tool [32]. Protein model structures from AlphaFold predictions have been integrated into the SWISS-MODEL web server [33,34]. For structural assessment, Qualitative Model Energy Analysis (QMEAN) [35] values around 0 showed reliability, while values ≤−4.0 indicated a lowquality model [36]. A normalised QMEAN4 score compares a predicted model structure with a non-redundant set of PDB structures, providing information on how many standard deviations away from the mean is the predicted model score given a score distribution from a large set of experimentally determined structures. The value varies from |Z-score| < 1 to |Z-score| > 2, with |Z-score| < 1 being closer to the mean, hence a better structure. Ramachandran plot and statistics were obtained, with a model having >98% of its residues in the most favoured regions considered to be of good quality. Similarly, Ramachandran outliers should be <0.2%, Rotamer outliers <1%, number of C-Beta deviations and twisted prolines/non-prolines should be zero. Further assessment was carried out using ERRAT [37] and VERIFY3D [38,39] on SAVES servers (v6.0). For VERIFY3D and ERRAT, models with ≥80% of its amino acid residues having average 3D/1D profile scores ≥0.2 and overall quality factor >90%, respectively, were considered to be of good quality.

Active Site Prediction
The active site residues necessary for the catalytic activity of AgTre were predicted using CASTp [40] and PrankWeb servers [41]. Active site prediction for trehalase from Ae. aegypti (AlphaFold model of Q16V81: AF-Q16V81-F1) was carried out using PrankWeb.

Ligand Preparation
Validamycin A, trehalozin, trehazolamine, validoxylamine A, casuarine, and uniflorine have been reported to be inhibitors of trehalase in different species. Hence, they were used as the starting compounds to generate more similar structures [9]. However, validamycin A was employed as the control ligand or reference standard due to its commercial availability and reported activity in different insect species. Compounds with similar structures to the control ligand and the previously reported trehalase inhibitors in different species were searched for and retrieved with a 100% search of the PubChem database [42] using a Tanimoto threshold of 80%. This search resulted in a total of 3022 compounds, which were downloaded in their SDF formats. Using the OpenBabel package in PyRx software, the SDF files were minimised using the Universal Force Field (uff) [43] to obtain appropriate proper bond lengths between the different atoms. They were then converted into the pdbqt format and further used for the virtual screening studies.

Protein Preparation
The best modelled protein structure of AgTre and the downloaded structure of Ae. aegypti (AlphaFold: AF-Q16V81-F1) were prepared using UCSF Chimera [44]. The proteins were minimised using 2000 steepest descent steps and 500 conjugate gradient steps with a conjugate gradient step size of 0.02. AMBER ff14SB charges were added to standard amino acid residues. The prepared PDB file was converted to pdbqt using AutoDock in PyRx.

Virtual Screening
Virtual screening of the 3022 compounds against AgTre was carried out using AutoDock Vina wizard in PyRx [45]. The vina search space was fitted around the predicted active site residues leading to xyz dimensions of 28 30.0916, and xyz centre of −2.1668, −2.0597, and −6.4698 Å. When the virtual screening was completed, binding energies were obtained for the starting compounds, control ligand, and top 9 ligands. Discovery Studio was used to visualise the Trehalase-ligand binding interaction [46].

Toxicity Screening and Ligand Optimisation
The control ligand, previously reported inhibitors, and top 9 hits were screened for possible ecotoxic effects using admetSAR 2.0 [47], toxicity using ADMETlab 2.0 [48], and their medicinal chemistry properties were determined using SwissADME [49]. ADME-Topt [50] was used to optimise the scaffolds of some top hits with the aim of seeing if optimised compounds would have a better binding affinity. A non-toxic constraint (i.e., compounds showing no carcinogenicity, Ames toxicity, acute oral toxicity) was used for optimisation. The new compounds were docked against AgTre and AaTre.

Molecular Dynamics Simulation
The MD simulation was carried out by using Nanoscale Molecular Dynamics (NAMD) [51] and Visual Molecular Dynamics (VMD) [52]. This helped to determine the stability of the best poses of some of the best hits from the docking studies in the binding site of AgTre. Energy minimisation of 1000 steps was performed to fix the backbone atoms, while the production simulation run was carried out for 1 ns, equivalent to 500,000 steps. The simulation was performed at a constant pressure of 1 atm and temperature of 310 K using Periodic Boundary conditions. The protein structure file (PSF) of the target was generated separately from that of the ligand using VMD, while those of ligands were generated using the Charmm36 forcefield of the Charmm-GUI web server [53]. These were generated to define the bond types, bond angles, atom types, and the number of molecules in the simulation system. The topologies (PDB and PSF) of the protein and ligands were merged, and the complexes were solvated using VMD to generate the cubic water box. The other necessary parameters (time and Periodic Boundary conditions) for the simulation were defined in a script and run using NAMD. The RMSD (Root Mean Square Deviation), RMSF (Root Mean Square Fluctuation), and PCA (Principal Component Analysis) of the simulation results were determined using Bio3D on the Galaxy Europe platform [54], while the hydrogen bond (h-bond) analysis was carried out using VMD software. Graphs were plotted using the plot() function in R-programming software [55].

Structure Prediction and Assessment of AgTre
The 3D structure of AgTre predicted by RoseTTAFold had a confidence score of 0.85, which is close to 1, suggesting a good structure. The best model from I-TASSER had a C-score of −0.22, while that from C-I-TASSER had a C-score of 0.45, suggesting that the model obtained from C-I-TASSER was a better model than that obtained from I-TASSER. I-TASSER and C-I-TASSER models had QMEAN < −4.0, while other models were >−4.0. However, the RoseTTAFold model (−0.10) had QMEAN closest to 0, suggesting it to be of better quality (Table 1). Similarly, it had a normalised QMEAN4 of |Z-score| < 1 as compared to the AlphaFold model with 1 < |Z-score| < 2 ( Figure 1A,B). This further suggests that the RoseTTAFold predicted model is the best comparable model with a nonredundant set of PDB structures. The predicted structure from AlphaFold had 96.13% of its residues in Ramachandran favoured, 0.70% Ramachandran outliers, and 0.82% Rotamer outliers compared to 98.06%, 0.35%, and 0%, respectively, observed for the predicted structure from RoseTTAFold in this current study (Table 1). When the RoseTTAFold model and AlphaFold model are aligned, the RMSD, which gives the average deviation between the corresponding atoms between the two proteins, was 0.896 Å when considering 476/570 residues, suggesting the two proteins had similar folds for most parts ( Figure 1C). Hence, the structure in this study is comparable to that predicted using AlphaFold. Among all the models tested, the model from RoseTTAFold prediction had the highest number of residues in the Ramachandran favoured region (98.06%). It was the only model with >98% residues in the Ramachandran favoured region ( Figure 1D). Based on this parameter, the models from I-TASSER (77.18%) and C-I-TASSER (79.67%) were inferior to all other models. While all models passed ERRAT and VERIFY3D, the RoseTTAFold model had no twisted prolines, cis prolines, rotamer outliers, and C beta deviations as compared to the other models (Table 1). Hence the RoseTTAFold model was observed to be the best model and was used subsequently in this study.

Active Site Prediction
The predicted active site residues of the RoseTTAFold AgTre model and the AaTre model are presented in Table 2. Active site prediction from CASTp yielded a total of 116 pockets, while the prediction from PrankWeb gave 12 pockets for AgTre and 5 pockets for AaTre. The pocket with the highest probability and pocket score was selected for Prank-Web, while the pocket with the largest area and volume was selected for CASTp. The probability for the selected active site pocket for AgTre and AaTre predicted using Prank-Web were 0.948 and 0.855, respectively. For AgTre, 17 amino acids were predicted to be necessary for catalytic activity in the selected pockets from both CASTp and PrankWeb.

Active Site Prediction
The predicted active site residues of the RoseTTAFold AgTre model and the AaTre model are presented in Table 2. Active site prediction from CASTp yielded a total of 116 pockets, while the prediction from PrankWeb gave 12 pockets for AgTre and 5 pockets for AaTre. The pocket with the highest probability and pocket score was selected for PrankWeb, while the pocket with the largest area and volume was selected for CASTp. The probability for the selected active site pocket for AgTre and AaTre predicted using PrankWeb were 0.948 and 0.855, respectively. For AgTre, 17 amino acids were predicted to be necessary for catalytic activity in the selected pockets from both CASTp and PrankWeb. Four additional residues were predicted to be active site residues by CASTp only. Since the predicted binding pocket from the two programs overlapped, the predicted active site residues from CASTp were used to define the grid box used for molecular docking. Twenty-one (21) residues were predicted as active site residues of AaTre. In the two species, 19 predicted active site residues of trehalase were conserved. This suggests the possibility of both organisms interacting with similar compounds.

Molecular Docking
The results obtained from the virtual screening and post-screening studies, ADMET studies, and lead optimisation using ADMETopt are presented in this subsection. Additionally, results comparing binding energies (obtained from virtual screening) between AgTre and AaTre are reported.

Virtual Screening and Post-Screening Studies
The virtual screening studies revealed that all the previously reported trehalase inhibitors that were screened had lower binding energy (i.e., higher affinity) than the control ligand, validamycin A (−5.4 kcal/mol, indicating lowest affinity) for AgTre (Table 3). Among the previously reported trehalase inhibitors, validoxylamine A had the lowest binding energy (−8.8 kcal/mol), indicating the highest affinity for AgTre compared to validamycin A, followed by trehazolin (−7.5 kcal/mol). Uniflorine, casuarine, and trehazolamine had quite similar binding energies for AgTre (−6.9, −6.6, and −6.4 kcal/mol), suggesting that AgTre might have similar affinity for these compounds. Of all the 3022 compounds screened, only 9 compounds had binding energies <−8.8 kcal/mol observed for validoxylamine A (the previously reported trehalase inhibitor with the lowest binding energy). The binding energies were between −8.9 to −9.1 kcal/mol ( Table 3). The compounds interacted with different residues in the active site of AgTre through conventional hydrogen bonds, carbon-hydrogen bonds, Pi interactions, and van der Waals interactions. The ligands formed conventional hydrogen bonds, Pi-donor hydrogen bonds, or carbon-hydrogen bonds with two or more of the active site residues Glu258, Asn186, Gln197, Asp292, Ser259, Gly290, Gln438, Arg256, Arg142, Tyr192, Ala287, Phe143, Trp149, Lys26. However, the control ligand, validamycin A, did not bind to any predicted active site residue of the protein. Validamycin A rather formed hydrogen bonds with Gly504, Gln506, Asn307, Ser266 (Table 3, Figures S1 and S2). Other ligands further interacted with other active site residues through van der Waal interactions. Validoxylamine, trehazolamine, 10690241, and 67837201 formed Pi-Sigma interactions with Tyr507, while 25023458 formed Pi-Sigma interactions with Phe143. In addition, validoxylamine, 101104782, 13364642, 25023458, 21021639, and 58618560 had Alkyl/Pi-Alkyl interactions with one or more of Tyr147, Phe143, Arg142, Tyr507, or Lys26. However, the only interaction validamycin A had with an active site residue was a van der Waal interaction with Tyr507. This further indicates that validamycin is a poor binder of AgTre.

ADMET Studies and Lead Optimisation Using ADMETopt
ADME screening of the top nine hit compounds, control ligand, and previously reported inhibitors using SwissADME revealed that they were not blood-brain barrier permeants or CYP 450 inhibitors. The compounds were predicted to be soluble and had no structural alert for Pan-Assay INterference compoundS (PAINS), suggesting that they are not promiscuous binders but would be target-specific [56,57]. Of the top 9 hit compounds, only 10406567, 10690241, and 67837201 had no structural alert for Brenk ( Table 4), suggesting that they are non-toxic and metabolically stable [58]. All hit compounds had synthetic accessibility scores lower than validamycin A, which had the highest synthetic accessibility score of 6.34. This suggests that these compounds might be easier to synthesise than validamycin A. The control ligand, starting compounds, and top 9 hits were predicted by admetSAR to be non-toxic to honeybees, crustaceans, and fish. The compounds were predicted to be biodegradable except casuarine and trehazolamine, which had a mild nonbiodegradable probability of 0.575 and 0.55, respectively (Table 4). Likewise, employing ADMETlab, all screened compounds were predicted to be non-carcinogenic, cause no irritation or corrosion to the eyes, and cause no skin sensitisation. Predictions by ADMETlab revealed that only 10406567, 10690241, and 67837201 among the top 9 hits had probabilities <0.35 for DILI and respiratory toxicity. Among these three compounds, 10406567 and 10690241 were predicted to have a low probability of being nephrotoxic by admet-SAR (probability of 0.4811 and 0.5058, respectively), while 67837201 was not nephrotoxic. Uniflorine (9794258) was predicted to be the safest previously reported compound.

ADMET Studies and Lead Optimisation Using ADMETopt
ADME screening of the top nine hit compounds, control ligand, and previously reported inhibitors using SwissADME revealed that they were not blood-brain barrier permeants or CYP 450 inhibitors. The compounds were predicted to be soluble and had no structural alert for Pan-Assay INterference compoundS (PAINS), suggesting that they are not promiscuous binders but would be target-specific [56,57]. Of the top 9 hit compounds, only 10406567, 10690241, and 67837201 had no structural alert for Brenk (Table 4), suggesting that they are non-toxic and metabolically stable [58]. All hit compounds had synthetic accessibility scores lower than validamycin A, which had the highest synthetic accessibility score of 6.34. This suggests that these compounds might be easier to synthesise than validamycin A. The control ligand, starting compounds, and top 9 hits were predicted by admetSAR to be non-toxic to honeybees, crustaceans, and fish. The compounds were predicted to be biodegradable except casuarine and trehazolamine, which had a mild non-biodegradable probability of 0.575 and 0.55, respectively (Table 4). Likewise, employing ADMETlab, all screened compounds were predicted to be non-carcinogenic, cause no irritation or corrosion to the eyes, and cause no skin sensitisation. Predictions by ADMETlab revealed that only 10406567, 10690241, and 67837201 among the top 9 hits had probabilities < 0.35 for DILI and respiratory toxicity. Among these three compounds, 10406567 and 10690241 were predicted to have a low probability of being nephrotoxic by admetSAR (probability of 0.4811 and 0.5058, respectively), while 67837201 was not nephrotoxic. Uniflorine (9794258) was predicted to be the safest previously reported compound.

Molecular Dynamics Simulation
The control compound (validamycin A), the previously reported trehalase inhibitor with the best affinity for trehalase (validoxylamine A), as well as 67837201, Compound 1, Compound 2, and 10406567 in complex with AgTre, were subjected to MD simulation.

Root Mean Square Deviation (RMSD) and Root Mean Square Fluctuation (RMSF)
RMSD and RMSF results from the molecular dynamics simulation studies are shown in Figure 2A-D, and the individual results are in Figures S4-S8. The stability of protein, indicated by the RMSD changes with respect to Frame No, is presented in Figure 2A. Similarly, the average RMSDs of the protein and standard deviations are shown in Table 7. The binding of validamycin A, 67837201, 10406567, and Compound 1 to AgTre resulted in smaller RMSD changes in the protein compared to the binding of validoxylamine A (the previously reported trehalase inhibitor with the best affinity for trehalase). On the contrary, the binding of Compound 2 resulted in more conformational changes to the protein structure. However, all conformation changes were less than 3 Å, suggesting that the conformation of the protein was stable throughout the MD simulation [59]. RMSF describes the flexibility of each amino acid residue in the protein during the MD simulation. The larger the RMSF value is, the more flexible the residue is. Regarding Figure 2B, the amino acid residues with RMSF > 1.5 Å represent the flexible residues in the loop region of the protein. Similar flexible patterns in the protein were observed upon binding of the different ligands. The RMSDs of the ligands (validamycin A, 67837201, 10406567, Compound 2, and validoxylamine A) are shown in Figure 2C, while that of Compound 1 is shown in Figure 2D. Similarly, the average RMSDs of the ligands and standard deviations are shown in Table 7. The average RMSDs of validoxylamine A, validamycin A, 67837201, and 10406567 reveal that these compounds do not significantly change their orientation during the simulation; hence, they form stable complexes with the protein. This is also confirmed in the histogram ( Figure S7). However, Compound 1 was unstable with large changes in its RMSD, having an average RMSD of 5.46 Å with a standard deviation of 1.40 Å (Table 7). Compound 2 showed slight changes early in the simulation but returned to conformations similar to the docking conformation. It had an average RMSD of 1.17 Å and a standard deviation of 0.23 Å. The compound (10406567) exhibited the best stability during the MD simulation with an average RMSD of 0.96 Å and a standard deviation of 0.07 Å. This was followed by Compound 2 and 67837201; 67837201 had an average RMSD of 1.46 Å with a standard deviation of 0.10 Å. The average RMSDs for these compounds was lower than those of validamycin A and validoxylamine A, suggesting that the three compounds maintained relatively stable conformations throughout the simulation period. a standard deviation of 0.23 Å. The compound (10406567) exhibited the best stability during the MD simulation with an average RMSD of 0.96 Å and a standard deviation of 0.07 Å. This was followed by Compound 2 and 67837201; 67837201 had an average RMSD of 1.46 Å with a standard deviation of 0.10 Å. The average RMSDs for these compounds was lower than those of validamycin A and validoxylamine A, suggesting that the three compounds maintained relatively stable conformations throughout the simulation period.    Figure 3. This suggests that the molecular dynamics simulation procedure captured the major or dominant motions rather than the much less dominant ones. Considering the first three principal components for 67837201, 10406567, Compound 1, and Compound 2, the order of variance was 10406567 < Compound 1 < 67837201 < Compound 2. PC1 accounted for 41.9, 36.8, 40.1, 32.7, 49.7, and 42.4% of the variance for validoxylamine A, validamycin A, 67837201, and 10406567, Compounds 2 and 1. When PC1 alone was considered, the order of variance was 10406567 < 67837201 < Compound 1 < Compound 2. Of the hit compounds screened, Compound 2 had the highest variance, while 10406567 had the lowest variance either when PC1 only or the first three principal components were considered. 10406567 < Compound 1 < 67837201 < Compound 2. PC1 accounted for 41.9, 36.8, 40.1, 32.7, 49.7, and 42.4% of the variance for validoxylamine A, validamycin A, 67837201, and 10406567, Compounds 2 and 1. When PC1 alone was considered, the order of variance was 10406567 < 67837201 < Compound 1 < Compound 2. Of the hit compounds screened, Compound 2 had the highest variance, while 10406567 had the lowest variance either when PC1 only or the first three principal components were considered.

Hydrogen Bond Analyses
The percentage occupancy of the h-bond >10% between the compounds and amino acid residues of the protein is presented in Table 8. For validoxylamine A, 67837201, and 10406567, Compound 1, and Compound 2, multiple hydrogen bonds were identified between the active site of the protein and the ligand. However, validamycin A did not form hydrogen bonds with any of the active side residues throughout the simulation. The compounds formed h-bonds with similar residues they interacted with in the molecular docking studies (either through hydrogen bond, Pi, unfavourable donor-donor, or van der Waals interactions). Compounds 1 and 2 formed hydrogen bonds with Arg195 (46.26%) and Asp150 (57.72%), respectively, which were not observed in the molecular docking studies. The hydrogen bond interaction with Asp292 was persistent for validoxylamine A (75.46%) and 10406567 (54.56%), while it had an occupancy of 22.62% for 67837201. However, it had an h bond occupancy of <10% for Compound 1 (5.92%) and 2 (6.87%). Glu258 had an h-bond occupancy of 59.33% for 67837201, Asn187 had an occupancy of 58.67% for Compound 1, and Asp150 had an occupancy of 57.72% for Compound 2. Hydrogen bond interactions for all the compounds with other residues were also observed, but these were not persistent (<10%) and were only present for a minority of the time-lapse of the simulation.

Discussion
In this study, validamycin A was observed to have a lower affinity for AgTre and AaTre than validoxylamine A and other previously reported compounds tested (Tables 3 and 6). These results support previous studies in other organisms suggesting validamycin A to be a poor binder of trehalase compared to validoxylamine A. For example, Asano et al. [60] reported an IC 50 of 72 µM for validamycin A in Rizochtonia solani, compared to an IC 50 of 140 nM for validoxylamine A. Similarly, validoxylamine A has been observed to have a better inhibitory effect on trehalase than validamycin A in Spodoptera litura (IC 50 of 48 nM vs. 370 nM, respectively), porcine intestine (14 nM vs. 120 nM, respectively) [61], porcine kidney (2.4 nM vs. 250 µM) [62], and termites (3.2 µM vs. 402 µM) [63]. In a study by Asano et al. [64], injection of 10 µg of validoxylamine A in the last instar larvae of Spodoptera litura resulted in 70% death in prepupae and 30% abnormal pupae, with no normal pupae or adult being formed. However, injection of validamycin A resulted in 25% death in prepupae, 30% abnormal pupae, and 45% normal pupae, although only 11% of the normal pupae emerged as adults [64]. These studies affirm that validoxylamine A is a better inhibitor of trehalase than validamycin A, corroborating the observation in this study. In this present study, trehazolin (binding energy of −7.5 kcal/mol) was observed to be a better binder of AgTre than validamycin A. Ando et al. [65] reported trehazolin to have an IC 50 of 66 nM in Rizochtonia solani compared to IC 50 of 72 µM for validamycin A, and 140 nM for validoxylamine A reported earlier by Asano et al. [60], thus suggesting trehazolin is a better inhibitor than validamycin A and validoxylamine A. However, in porcine kidney, Kyosseva et al. [62] reported trehazolin to have an IC 50 of 16 nM compared to 2.4 nM for validoxylamine A and 250 µM for validamycin A. This suggests that describing trehazolin as a better trehalase inhibitor when compared to validoxylamine might be organism-dependent. In this current study, validoxylamine was observed to be the best inhibitor of AgTre among the previously reported inhibitors, followed by trehazolin. Uniflorine, trehazolamine, and casuarine had quite similar binding affinities for AgTre, while validamycin A had the least affinity (Table 3). Unlike other compounds which interacted with active site residues of AgTre, validamycin A interacted with other residues (Table 3, Figure S1). Unlike AgTre, in which validamycin A had van der Waal interactions with tyrosine 507 only in the active site, validamycin A had a carbon hydrogen bond with Arg 193, van der Waals interactions with Glu554, Lys77, and Tyr555 in the active site of AaTre. However, conventional hydrogen bonds were formed with the non-active site residues Gly192, Glu313, and Arg312. In both organisms, validamycin A had the highest binding energy among the tested ligands (Table 6). This might explain the low larvicidal activity observed with validamycin treatment in Ae. aegypti by Marten et al. [19]. MD simulation was employed to investigate the stability of the protein-ligand complexes in a flexible system, unlike molecular docking, which assumes that the protein and ligand are in a rigid conformation in most cases [66]. Despite the low affinity for trehalase and not interacting much with active site residues, validamycin A was relatively stable during AgTre-validamycin A MD simulation ( Figure 2C, Table 7). It had an h-bond occupancy of 33.97% for Glu262 and 13.03% for Gly504 (Table 8), while the h-bonds formed with other residues persisted below 10% during the whole simulation. These results further corroborate validamycin A as a poor binder of AgTre. This suggests that these other molecules, aside from validamycin A, could serve as better insecticides for malaria vector control.
Uniflorine, the safest previously reported compound showing no predicted toxicity (Table 4), had a binding energy of −6.9 kcal/mol for AgTre and −7.6 kcal/mol for AaTre (Table 6). Uniflorine has been reported to show >5649 selectivity for C. riparius trehalase (177 ± 18 nM) as compared to trehalase from porcine kidney (>1 mM) [67]. Uniflorine is a weaker binder of trehalase from porcine kidney when compared to validoxylamine A, for which an IC 50 of 2.4 nM has been reported, but a better binder than validamycin A, for which an IC 50 of 250 µM was reported [62]. In this present study, a similar trend was observed with uniflorine having a lower affinity for AgTre and AaTre than validoxylamine A (−8.8 kcal/mol) but a better affinity than validamycin A in both organisms. Considering the high selectivity of uniflorine and its safety, it could be further tested as a possible compound for malaria vector control.
In this study, Compounds 1 and 67837201 were the hit compounds with no toxicity and yet low binding energy for AgTre (Tables 3-5). Their binding energies were higher than that of uniflorine, suggesting that they are better binders of AgTre. RMSD results from further examination of the stability of their protein-ligand complexes by MD simulation suggested that Compound 1 might have a binding energy different from that observed in the molecular docking studies as the compound was unstable during the simulation ( Figure 2D). This highlights the importance of following up molecular docking studies with MD simulation, as MD simulation has been reported to help separate false positive hits from true positive hits obtained from molecular docking studies [21]. 67837201 was stable during the simulation ( Figure 2C); however, unlike in AgTre, where a binding energy of −8.9 kcal/mol was observed during docking studies, it had a higher binding energy for AaTre (−5.9 kcal/mol). Both 10406567 and Compound 2 had good affinity for trehalase of both organisms (Table 6) and displayed stability throughout the simulation ( Figure 2C), suggesting they could be excellent compounds for vector control. Further experimental studies to assess their toxicity are necessary, as both were predicted to have a slight probability of being nephrotoxic (Table 4). Similarly, their selectivity for these target organisms needs to be compared to that of non-target organisms in further experimental investigations. However, it has been reported that the disease arising from trehalase inhibition in humans is a mild disease, which can be reversed by avoiding the consumption of mushrooms [68], suggesting trehalase inhibitors as safer alternatives for vector control compared to, for example, organophosphates. Although trehalase enzymes are present in bacteria, fungi, and arthropods [69], uniflorine, 67837201, 10406567, and Compound 2 were predicted to be non-toxic to honeybees, crustaceans, and aquatic life. While these compounds were screened in silico and simulation performed at 310 K, reports of in vivo testing of validamycin A and validoxylamine A in different studies give confidence that the compounds identified in this study would be effective in vivo despite different climatic conditions in which mosquitoes may exist. Experimental validation remains to be carried out as the next future step.

Conclusions
This study screened trehalase inhibitors against AgTre to identify inhibitors with better affinity than validamycin A, which could be developed further for malaria vector control. The study proposes uniflorine, 67837201, 10406567, and Compound 2 as possible compounds for malaria vector control. These compounds also bind to AaTre; as such, they can be used to control Aedes aegypti, thus preventing transmission of several other diseases such as Mayaro and dengue fever. The trehalase inhibitors can be applied to breeding sites, as they could affect the development of mosquitoes and prevent their flight to humans, consequently preventing the transmission of mosquito-borne diseases.