Natural Products as New Treatment Options for Trichomoniasis: A Molecular Docking Investigation

Trichomoniasis, caused by the parasitic protozoan Trichomonas vaginalis, is the most common non-viral sexually-transmitted disease, and there can be severe complications from trichomoniasis. Antibiotic resistance in T. vaginalis is increasing, but there are currently no alternatives treatment options. There is a need to discover and develop new chemotherapeutic alternatives. Plant-derived natural products have long served as sources for new medicinal agents, as well as new leads for drug discovery and development. In this work, we have carried out an in silico screening of 952 antiprotozoal phytochemicals with specific protein drug targets of T. vaginalis. A total of 42 compounds showed remarkable docking properties to T. vaginalis methionine gamma-lyase (TvMGL) and to T. vaginalis purine nucleoside phosphorylase (TvPNP). The most promising ligands were polyphenolic compounds, and several of these showed docking properties superior to either co-crystallized ligands or synthetic enzyme inhibitors.


Introduction
Trichomoniasis is a sexually-transmitted disease (STD) caused by the parasitic protozoan Trichomonas vaginalis and is the most common non-viral STD with an estimated 3.7 million cases in the United States [1]. Only about 30% of individuals infected with T. vaginalis experience symptoms of genital discomfort, itching, burning or discharge, but there can be severe inflammation, increased risk of HIV infection, cervical cancer, preterm delivery and low birth weight [1]. Trichomoniasis can be treated with antibiotics, usually metronidazole or tinidazole, but there are increasing reports of resistance to these drugs [2]. There are currently no alternative drugs approved for the treatment of refractory cases of trichomoniasis, emphasizing the need for new treatment options. Recent investigations have identified several T. vaginalis proteins that may serve as targets for drug discovery and development [3,4].

Trichomonas vaginalis Protein Targets
Proteases are known to carry out important biological processes in parasitic protozoa and are therefore potential drug targets. Papain-like proteases have been suggested to be involved in nutrition and hemolysis, as well as able to induce apoptosis in human vaginal epithelial cells [3]. More than 40 papain-like cysteine proteases have been identified in the T. vaginalis degradome, which have been implicated as virulence factors [4].
Triosephosphate isomerase (TPI) is a glycolytic enzyme that catalyzes the interconversion of glyceraldehyde 3-phosphate and dihydroxyacetone phosphate and is an essential component of the glycolytic pathway [5]. Because of its importance in glycolysis, TPI has been identified as a good drug target for antiparasitic chemotherapeutics [6].
Lactate dehydrogenase (LDH) catalyzes the interconversion of lactate to pyruvate with concomitant interconversion of NAD + to NADH. LDH is a key enzyme in glycolysis and is found in nearly all living cells. Because T. vaginalis lactate dehydrogenase (TvLDH) is required for parasite survival, but is not similar to human LDH, TvLDH may be regarded as a suitable target for drug discovery [7,8].
Methionine gamma-lyase (MGL) has been characterized in several bacteria species, as well as the parasitic protozoans Entamoeba histolytica and Trichomonas vaginalis [9]. The enzyme degrades sulfur-containing amino acids to α-keto acids, ammonia and thiols and plays a key role in the regulation of sulfur-containing amino acids. Mammals do not have MGL, so this enzyme is a potential drug target for anti-Trichomonas chemotherapy.
Thioredoxin reductase (TrxR) catalyzes the reduction of thioredoxin, and the thioredoxin system provides a defense against oxidative damage due to oxygen metabolism and redox signaling [10]. Mammalian TrxRs and TrxRs from parasitic protozoa are different classes with different mechanisms of activity [11]. Because TrxR is a strong antioxidant that protects T. vaginalis from oxidative stress, the parasite lacks glutathione or catalase and T. vaginalis thioredoxin reductase (TvTrxR) is very different from human thioredoxin reductase (HsTrxR), TvTrxR has been identified as a target for trichomoniasis chemotherapy [12] and is the target of metronidazole and other nitroimidazole drugs [13].
Purine nucleoside phosphorylase (PNP) catalyzes the phosphorolysis of the N-glycosidic bonds of purine nucleosides (or deoxynucleosides) to give α-ribose-1-phosphate and the purine base and functions in the purine salvage pathway [14]. Purine salvage is essential for obligate parasitic protozoa, including T. vaginalis, and T. vaginalis purine nucleoside phosphorylase (TvPNP) has been identified as an attractive chemotherapeutic target [15].

Homology Modeling
In the absence of an experimentally-determined protein structure by crystallographic or nuclear magnetic resonance (NMR) methods, homology modeling can provide useful three-dimensional structures for proteins that are related to known protein structures. The premise is that the proteins have similar structures and binding and/or active sites of the proteins retain identical structures. Several computational methods for predicting protein structures based on homology models are currently available [16], and homology modeling has been shown to be a valuable tool for in silico screening of biomolecular targets [17]. In this work, we have used homology modeling to predict T. vaginalis protein structures for which there are no experimentally-determined structures.

Molecular Docking
Molecular docking is a well-accepted tool in drug discovery and complements X-ray crystallography and NMR spectroscopy in analyzing small molecule-protein interactions. In many cases, it has replaced high-throughput screening for initial investigations in lead discovery. Nevertheless, there are some limitations to the method that mostly arise from not accounting for local and global protein dynamics, as well as the inability to accurately predict ligand-protein covalent interactions and solvent accessibilities. The protein is typically modeled as a rigid structure without flexibility; solvation in the active or binding site is usually not included in the models, and free-energies of the ligand-protein complexes are generally ignored [18][19][20]. In spite of these limitations, molecular docking studies of natural product ligands with potential drug targets can serve to identify natural product drugs or drug leads to treat human infections [6]. In this work, we have carried out in silico screening of antiprotozoal natural products with several potential protein targets of Trichomonas vaginalis.

Homology Modeling
Homology models for each of the Trichomonas proteins that are not currently available from the Protein Data Bank (PDB) were constructed from crystal structure templates found in the Protein Data Bank using FASTA sequences downloaded from the National Center for Biotechnology Information's (NCBI) GenBank. Sequences with high sequence similarity in the PDB were identified with NCBI's BLAST utility using the BLOSUM80 scoring matrix (BLOcks SUbstitution Matrix). Sequences with high similarity to the reference sequences, as well as having good coverage for the active sites in the proteins, were identified using NCBI's BLAST utility with the BLOSUM80 scoring matrix and used as templates for single-reference homology modeling.
The protein sequences were first aligned to their respective template sequences using the BLOSUM62 substitution matrix and a protein backbone constructed and superposed to the reference structure using the protein alignment tool in Molecular Operating Environment, MOE 2014.0901. The homology modeling interface in MOE was used to generate a set of putative protein structures by aligning atomic coordinates of the amino acid sequence to those of the template sequence backbone and minimizing permutations of side chain orientations using the AMBER12:EHT force field [21][22][23] with reaction field solvation. The candidate structure with the lowest root-mean-square deviation of atomic positions (RMSD) deviation from the template backbone was selected and optimized using a constrained minimization.
Prior to docking, all solvent molecules and the co-crystallized ligands were removed from the structures. If co-factors were present, they were retained in each protein model (i.e., dihydroflavin adenine dinucleotide (FDA) in TvTrxR and 1,4-dihydronicotinamide adenine dinucleotide (NADH) in TvLDH). Molecular docking calculations for all compounds with each of the proteins were undertaken using Molegro Virtual Docker (Version 6.0.1, Molegro ApS, Aarhus, Denmark) [41], with a sphere (15 Å radius) large enough to accommodate the cavity centered on the binding sites of each protein structure in order to allow each ligand to search. If a co-crystallized inhibitor or substrate was present in the structure, then that site was chosen as the binding site. If no co-crystallized ligand was present, then suitably-sized (>25 Å 3 ) cavities were used as potential binding sites. Standard protonation states of the proteins based on neutral pH were used in the docking studies. Each protein was used as a rigid model structure; no relaxation of the protein was performed. Assignments of the charges on each protein were based on standard templates as part of the Molegro Virtual Docker program; no other charges were necessary to set. Our in-house library of antiprotozoal phytochemicals (obtained by searching the phytochemical literature and the Dictionary of Natural Products [42]), which was filtered for drug-like properties based on Lipinski's rule-of-five [43], was used for the molecular docking study. Overall, 952 antiprotozoal phytochemicals have been docked. This molecule set was comprised of 214 alkaloids, 369 terpenoids, 174 flavonoids and 195 polyphenolic compounds. Each ligand structure was built using Spartan'16 for Windows (Version 1.1.0, Wavefunction Inc., Irvine, CA, USA). For each ligand, a conformational search and geometry optimization was carried out using the MMFF (Merck Molecular Force Field) [44]. Flexible ligand models were used in the docking and subsequent optimization scheme. Variable orientations of each of the ligands were searched and ranked based on their re-rank score. For each docking simulation, the maximum number of iterations for the docking algorithm was set to 1500, with a maximum population size of 50 and 100 runs per ligand. The RMSD threshold for multiple poses was set to 1.00 Å. The generated poses from each ligand were sorted by the calculated re-rank score. In analyzing the docking scores, we have attempted to account for the recognized bias toward high molecular weight compounds [45][46][47][48][49][50], using the scheme: DS norm = 7.2 × E dock /MW 1 / 3 , where DS norm is the normalized docking score, E dock is the MolDock re-rank score, MW is the molecular weight and 7.2 is a scaling constant to bring the average DS norm values comparable to E dock [51]. As a test of docking accuracy and for docking energy comparison, co-crystallized ligands were re-docked into the protein structures (see Table 1). In addition, as a validation of the docking method, we have carried out docking of picomolar and nanomolar synthetic purine nucleoside phosphorylase inhibitors with human PNP and T. vaginalis PNP (see Table 2). The docking shows good docking scores and good docking pose orientations for these compounds serving to confirm the validity of the docking method.

Homology Modeling
Structures in the PDB with co-crystallized ligands in the active sites were chosen as templates for the homology models in order to retain essential binding site topology once the best candidate models were subjected to minimization with the AMBER force field. The Ramachandran plots for the homology models, along with the plots for the template structures are shown in Figures 1-3. The Ramachandran plot for modeled T. vaginalis thioredoxin reductase (TvTrxR) shows that the phi and psi angles cluster in the typical regions for helices and sheets ( Figure 1). As can be seen in the Ramachandran plots of the homology models, outlier residues are not located in the active sites of the target proteins, and the binding site residues lie within the allowed regions of the psi-phi angle islands. Thus, although there are three outliers in the homology model Ramachandran plot of TvTxR, there are none in the active site of this enzyme. As shown in the comparative sequences, the binding site residues are also principally conserved ( Figure 1). Likewise, the Ramachandran plots of the homology models of T. vaginalis cysteine proteases TvCPCAC1 ( Figure 2) and TvCP2 ( Figure 3) show the phi and psi angles to cluster in the typical regions for helices and sheets, particularly with the active sites, which showed no residual outliers. In addition, the active sites of the homology modeled protein structures are structurally very similar to the active sites of the proteins from which the models were based (Figures 4-6). In the case of the TvTrxR homology model, the RMSD between the crystallized ligand and that of the homology model is 0.703 Å. For the TvCPCAC1 homology model, the RMSD is 0.800 Å. Additionally, in the case of the TvCP2 homology model, the RMSD between the crystallized ligand and that of the homology model is 3.036 Å. Calculated RMSD values include all ligand atoms. Because the models in this study gave reliable backbone conformations, as well as residue interactions around the active sites, these homology models are deemed to be trustworthy, particularly in the regions where molecular docking takes place.
Sci. Pharm. 2017, 85, 5 5 the phi and psi angles to cluster in the typical regions for helices and sheets, particularly with the active sites, which showed no residual outliers. In addition, the active sites of the homology modeled protein structures are structurally very similar to the active sites of the proteins from which the models were based (Figures 4-6). In the case of the TvTrxR homology model, the RMSD between the crystallized ligand and that of the homology model is 0.703 Å. For the TvCPCAC1 homology model, the RMSD is 0.800 Å. Additionally, in the case of the TvCP2 homology model, the RMSD between the crystallized ligand and that of the homology model is 3.036 Å. Calculated RMSD values include all ligand atoms. Because the models in this study gave reliable backbone conformations, as well as residue interactions around the active sites, these homology models are deemed to be trustworthy, particularly in the regions where molecular docking takes place.

Molecular Docking Validation
In order to confirm the validity of the docking method using MolDock, those protein structures with co-crystallized ligands were re-docked to confirm the docking orientation. The docking energies and root mean squared deviations (RMSD, Å) are shown in Table 1. Those structures with relatively rigid co-crystallized ligands reproduced the ligand orientation very well. Thus, for example, TvMGL (PDB 1E5F) and TvPNP (PDB 1Z36) showed excellent re-docking properties with RMSD = 0.57 and 0.21 Å, respectively ( Figure 7). On the other hand, co-crystallized ligands with many rotatable bonds were not as successfully re-docked; protein structures with floppy cysteine protease inhibitors, such as HsCatK (PDB 1U9V) and HsCatL (PDB 3HWN), had RMSD values of 4.03 Å and 6.65 Å, respectively (see Figure 7). As an additional validation of the docking method, docking of synthetic picomolar and nanomolar synthetic purine nucleoside phosphorylase inhibitors with human PNP [52] and T. vaginalis PNP [53,54] were carried out. The docking of these synthetic inhibitors shows good docking scores, generally <−100 kJ/mol ( Table 2) and consistent docking pose orientations ( Figure 8) for these compounds compared with the co-crystallized ligands. Figure 6. Overlay of the protein structures of rabbit (Oryctolagus cuniculus) cathepsin K, PDB 2F7D [25] (red ribbon), and the homology model of Trichomonas vaginalis papain-like protein, TvCP2 (blue ribbon). The co-crystallized ligand is shown as a wireframe structure.

Molecular Docking Validation
In order to confirm the validity of the docking method using MolDock, those protein structures with co-crystallized ligands were re-docked to confirm the docking orientation. The docking energies and root mean squared deviations (RMSD, Å) are shown in Table 1. Those structures with relatively rigid co-crystallized ligands reproduced the ligand orientation very well. Thus, for example, TvMGL (PDB 1E5F) and TvPNP (PDB 1Z36) showed excellent re-docking properties with RMSD = 0.57 and 0.21 Å, respectively ( Figure 7). On the other hand, co-crystallized ligands with many rotatable bonds were not as successfully re-docked; protein structures with floppy cysteine protease inhibitors, such as HsCatK (PDB 1U9V) and HsCatL (PDB 3HWN), had RMSD values of 4.03 Å and 6.65 Å, respectively (see Figure 7). As an additional validation of the docking method, docking of synthetic picomolar and nanomolar synthetic purine nucleoside phosphorylase inhibitors with human PNP [52] and T. vaginalis PNP [53,54] were carried out. The docking of these synthetic inhibitors shows good docking scores, generally <−100 kJ/mol ( Table 2) and consistent docking pose orientations ( Figure 8) for these compounds compared with the co-crystallized ligands. Table 1. MolDock docking energies of co-crystallized ligands and root mean squared deviations between the co-crystallized ligand and the re-docked poses of the co-crystallized ligand with Trichomonas vaginalis and homologous human protein structures.

Molecular Docking of Phytochemicals
Of the 952 antiprotozoal phytochemicals examined in this molecular docking study, a total of 42 showed notable docking scores (<−125 kJ/mol) (Tables 3-5, Figure 9). The −125 kJ cut-off was chosen based on docking scores of co-crystallized ligands and synthetic inhibitors (see Tables 1 and 2). T. vaginalis cysteine proteases do not look to be promising targets for antiprotozoal phytochemicals. None of the phytochemical ligands in this study showed docking energies <−114 kJ/mol (Table 3). Furthermore, most of the ligands showed comparable or better docking to the human cathepsins than to T. vaginalis cysteine proteases. Likewise, T. vaginalis triosephosphate isomerase does not look to be a promising target for antiprotozoal phytochemicals. Only one ligand, the chalcone 2ʹ,4,4ʹ-trihydroxy-3ʹ,5ʹ-diprenylchalcone, showed promising docking to human triosephosphate isomerase with a normalized docking score of −128.4 kJ/mol (Table 4). There were no phytochemical ligands that showed docking scores to T. vaginalis lactate dehydrogenase more exothermic than −116 kJ/mol (Table 4); TvLDH does not look to be a promising target for antiprotozoal phytochemicals. Similarly, the lowest docking score for T. vaginalis thioredoxin reductase was −119.8 kJ/mol (Table 4), so TvTrxR cannot be regarded as a potential target for the phytochemicals examined.

Conclusions
Drug resistance of Trichomonas vaginalis is increasing, and trichomoniasis can be regarded as a re-emerging infectious disease. There is a need for new chemotherapeutic agents to treat trichomoniasis, and natural products are an attractive resource. Molecular docking of antiprotozoal phytochemical agents has revealed two potential protein drug targets of T. vaginalis, methionine gamma-lyase (TvMGL) and purine nucleoside phosphorylase (TvPNP). The best docking ligands were polyphenolic compounds, including aurones, chalcones, flavonoids and lignans. Several phytochemicals showed docking properties superior to co-crystallized ligands or synthetic enzyme inhibitors. This preliminary computational investigation predicts several phytochemicals as potential inhibitors of T. vaginalis. However, at this time, we have not provided experimental evidence, in vitro or in vivo, that these predictions will necessarily lead to effective treatments. Additional experimental validation of these predictions is necessary; experimental validation experiments are underway in our laboratories.