Targeting Glutaminase by Natural Compounds: Structure-Based Virtual Screening and Molecular Dynamics Simulation Approach to Suppress Cancer Progression

Cancer cells change their glucose and glutamine (GLU) metabolism to obtain the energy required to continue growing. Glutaminase (GLS) plays a crucial role in promoting cell metabolism for cancer cell growth; targeting GLU metabolism by inhibiting GLS has attracted interest as a potential cancer management strategy. Herein, we employed a sequential screening of traditional Chinese medicine (TCM) database followed by drug-likeness and molecular dynamics simulations against the active site of GLS. We report 12 potent compounds after screening the TCM database against GLS, followed by a drug-likeness filter with Lipinski and Veber rule criteria. Among them, ZINC03978829 and ZINC32296657 were found to have higher binding energy (BE) values than the control compound 6-Diazo-5-Oxo-L-Norleucine, with BEs of −9.3 and −9.7 kcal/mol, respectively, compared to the BE of 6-Diazo-5-Oxo-L-Norleucine (−4.7 kcal/mol) with GLS. Molecular dynamics simulations were used to evaluate the results further, and a 100 ns MD simulation revealed that the hits form stable complexes with GLS and formed 2–5 hydrogen bond interactions. This study indicates that these hits might be employed as GLS inhibitors in the battle against cancer. However, more laboratory tests are a prerequisite to optimize them as GLS inhibitors.


Introduction
Cancer is the leading cause of mortality globally [1]. In the United States, approximately 2 million new cancer cases and 0.6 million cancer fatalities are expected in 2022, with 350 deaths a day from lung cancer, the biggest cause of cancer mortality [2]. Cellular metabolic reprogramming is a well-known feature of cancer to sustain the high demand for energy required for growth [3]. Cancer cells change their glucose and glutamine (GLU) metabolism to get the energy they need to keep growing. Enhanced GLU metabolism is a characteristic of cancer and is recognized as a major metabolic alteration in cancer cells [4]. Many cancers rely on GLU catabolism to generate metabolic intermediates that support bioenergetics and biosynthetic needs [5]. Glutaminase (GLS) is a mitochondrial enzyme that starts this process by converting GLU to glutamate, which is then used in a number of activities that help tumor cells grow and survive, including energy production, amino acid synthesis, and glutathione formation. GLS activity corresponds with poor clinical prognosis in patients with breast, colorectal, lung, liver, and brain tumors, and rapid-growing tumor cells have enhanced GLS mRNA levels as well as elevated GLS expression [4,[6][7][8][9]. GLS Molecules 2022, 27, 5042 2 of 10 inhibition by siRNA has antiproliferative effects on cancer cells [10]. There is currently no FDA-approved GLS inhibitor on the market, and CB-839 (telaglenastat) is the only GLS inhibitor tested for clinical trials [11].
Natural compounds with a wide chemical diversity have been studied intensively for their anticancer potential. Efforts have been undertaken in recent decades to identify new natural compounds from bacteria, plants, and other living species. Natural products are thought to have accounted for about 25% of all newly authorized anticancer medications between 1981 and 2019 [12,13]. Irinotecan, vincristine, etoposide, and paclitaxel are plant-derived anticancer medicines; actinomycin D and mitomycin C are bacteria-derived anticancer therapeutics; and bleomycin is a marine-derived anticancer medication [14].
Drug development and discovery is a difficult, slow, and expensive process, but it has been speeding up because of the development of computational tools and approaches. The use of computer-assisted drug design (CADD) tools in early investigations by prominent pharmaceutical corporations and research organizations have aided in the acceleration of the drug discovery process while reducing expenses and failures in the final stage [15]. The use of rational drug design as a component of CADD gives helpful insights into the binding affinity and molecular interaction between target protein and ligand. In addition, supercomputing facilities, advanced programs, algorithms, and tools have aided lead discovery in pharmaceutical research [16]. This work aimed to identify new possible compounds from a traditional Chinese medicine (TCM) database using in silico methodologies that might be employed as GLS inhibitors to combat cancer.

Protein Preparation
The 3D structure of GLS was obtained from the Protein Data Bank (PDB ID: 4O7D). It is a homo 4-mer structure, and we selected a monomer unit of GLS and minimized it for further screening.

Compound Library Preparation
The natural compound library comprising about 60 thousand compounds was downloaded from the TCM database. To prepare ligands (energy minimization and force field), the 'Prepare Ligands module' in Discovery Studio [17] 2020 was used. The control compound 6-Diazo-5-Oxo-L-Norleucine (DON; CID: 9087) was retrieved from the PubChem database.

Virtual Screening (VS) and Molecular Docking
VS plays an important role in the identification and development of new drugs [18]. Molecular docking-based VS is based on interacting receptors and small compounds. AutoDock Vina 1.1.2 and AutoDock 4.2.5.1 [19] were used for VS and in-depth docking studies. Top leads were retrieved based on the binding energy (BE) toward GLS.

Physiochemical and Drug-Likeness Properties Prediction
To estimate the physiochemical and drug-likeness properties, the top 20 compounds were selected from the VS results. The general physicochemical and drug-likeness properties of these compounds were analyzed using datawarrior tools [20], and these compounds were then filtered using Lipinski and Veber rules.

Molecular Dynamics Simulation
GROMACS 5.1.2 [21] was used to perform MD simulations on GLS-DON (control), GLS-ZINC32296657, and GLS-ZINC03978829 at 300 K, with the GROMOS96 43a1 forcefield [22]. The PRODRG server was used to produce the compound's topology and forcefield parameters [23]. Charges were manually corrected in the topology file, new compound atoms were added to the complex topology files, and all of the compounds' attributes were included in the system topology. GLS-DON (control), GLS-ZINC32296657, and GLS-ZINC03978829 were immersed in a 'cubic box' of water molecules with an initial diameter of 8 nm using the 'gmx editconf' module for boundary conditions and the 'gmx solvate' module for solvation. The charges on GLS-DON and screened complexes were neutralized by adding Na + and Cl − ions to preserve neutrality and a physiological concentration using the gmx genion module (0.15 M). All graphical representations of the 3D models were created using PyMOL and VMD [24].

Results and Discussion
GLU metabolism has long been a focus of cancer research due to its role in stressed settings. Because of GLS's critical involvement in promoting cell metabolism for tumor cell proliferation, targeting GLU metabolism via GLS inhibition has received a lot of interest as a possible therapeutic approach for cancer management [4]. This paper screened the TCM database of natural compounds against GLS. We chose the top 20 compounds to further investigate their physicochemical and drug-likeness properties ( Table 1). These (20 compounds) were further filtered using Lipinski and Veber rules to narrow down the screening results, and therefore only 13 compounds, including control, passed the prerequisites ( Table 2). By analyzing the binding poses of these 13 compounds with the GLS active site, it was determined that ZINC03978829 and ZINC32296657 have a strong affinity for GLS. ZINC03978829 was noted to interact with Lys245, Tyr249, Ser286, Lys289, Asn335, Val484, Leu505, and Gly509 residues of GLS ( Figure 1), while Tyr249, Ser286, Lys289, Asn319, Asn335, Glu381, Ser384, Arg387, Asn388, Tyr414, and Gly509 residues of GLS were found to interact with ZINC32296657 ( Figure 2). The BEs of ZINC03978829 and ZINC32296657 with GLS were noted to be −9.3 and −9.7 kcal/mol, respectively (Table 3).
Molecules 2022, 27, 5042 4 of 10 By analyzing the binding poses of these 13 compounds with the GLS active site, it was determined that ZINC03978829 and ZINC32296657 have a strong affinity for GLS. ZINC03978829 was noted to interact with Lys245, Tyr249, Ser286, Lys289, Asn335, Val484, Leu505, and Gly509 residues of GLS ( Figure 1), while Tyr249, Ser286, Lys289, Asn319, Asn335, Glu381, Ser384, Arg387, Asn388, Tyr414, and Gly509 residues of GLS were found to interact with ZINC32296657 ( Figure 2). The BEs of ZINC03978829 and ZINC32296657 with GLS were noted to be −9.3 and −9.7 kcal/mol, respectively (Table 3).   DON, a glutamine analog, is a competitive inhibitor of GLS [25], and was used as the reference compound in the present study. DON was observed to interact with Tyr249, Gln285, Ser286, Asn335, Glu381, Asn388, Tyr414, Tyr466, Gly483, and Val484 residues of GLS ( Figure 3). Tyr249, Ser286, and Asn335 were the common interacting residues of GLS with ZINC03978829, ZINC32296657, and the DON (Figures 1-3). The BE of DON with GLS was noted to be −4.7 kcal/mol (Table 1).   Figure 1. Visualization of ZINC03978829 (stick representation) in the binding pocket of GLS (a), interacting residues of GLS with ZINC03978829 (b,c).  Figure 1. Visualization of ZINC03978829 (stick representation) in the binding pocket of GLS (a), interacting residues of GLS with ZINC03978829 (b,c).   DON, a glutamine analog, is a competitive inhibitor of GLS [25], and was used as the reference compound in the present study. DON was observed to interact with Tyr249, Gln285, Ser286, Asn335, Glu381, Asn388, Tyr414, Tyr466, Gly483, and Val484 residues of GLS ( Figure 3). Tyr249, Ser286, and Asn335 were the common interacting residues of GLS with ZINC03978829, ZINC32296657, and the DON (Figures 1-3). The BE of DON with GLS was noted to be −4.7 kcal/mol (Table 1). DON, a glutamine analog, is a competitive inhibitor of GLS [25], and was used as the reference compound in the present study. DON was observed to interact with Tyr249, Gln285, Ser286, Asn335, Glu381, Asn388, Tyr414, Tyr466, Gly483, and Val484 residues of GLS ( Figure 3). Tyr249, Ser286, and Asn335 were the common interacting residues of GLS with ZINC03978829, ZINC32296657, and the DON (Figures 1-3). The BE of DON with GLS was noted to be −4.7 kcal/mol (Table 1). A docked complex with a higher (negative) BE is often viewed as an effective proteinligand interaction. As an outcome, the rate of dissociation of specific ligands from their binding targets would be slower, and their half-lives would be longer [26]. Interestingly, hits A docked complex with a higher (negative) BE is often viewed as an effective proteinligand interaction. As an outcome, the rate of dissociation of specific ligands from their binding targets would be slower, and their half-lives would be longer [26]. Interestingly, hits (ZINC03978829, ZINC32296657) exhibit higher BE with GLS than the DON (reference compound), indicating that these hits can be used as GLS inhibitors to fight cancer.
Bioinformatics plays an increasingly indispensable role in the integration of biological data as the diversity and complexity of data types have grown. For mining biologically active substances, there are now two types of data that are useful and accessible. The first is experimental biological activity data, such as anti-tuberculosis, anti-cancer, and anti-HIV biological activity data, which may be found in the PubChem database [27]. The second is curated data from multiple TCM databases regarding TCM plants and their derived compounds. Using bioinformatics to mine candidate molecules with diverse actions, the two forms of data present a new possibility [28]. TCM databases using the bioinformatics approaches have been used to identify inhibitors against many therapeutically important biological targets [17,[29][30][31].
The root mean square deviation (RMSD) is a protein stability metric; the lower the deviations, the more stable the protein structure [32]. GLS-DON (control), GLS-ZINC32296657, and GLS-ZINC03978829 had RMSD average values of 0.25, 0.33, and 0.31 nm, respectively. The RMSD plot revealed that GLU-ZINC03978829 binding increased GLS stability and resulted in fewer structural aberrations from its normal conformation. The bound structure of the GLS-ZINC32296657 complex is highly deviated (Figure 4a); it showed that the catalytic pocket of GLS enzyme was not forming a strong interaction with ZINC32296657, therefore it showed high deviation. In addition, the ligand RMSD also showed that DON and ZINC03978829 bind better than ZINC32296657 and are more stable (Figure 4b). Initially, ZINC32296657 formed a strong interaction with GLS until 22 ns, and after that the interaction started to lessen due to unstable bonding with the enzyme.
Molecules 2022, 27, 5042 6 of 10 (ZINC03978829, ZINC32296657) exhibit higher BE with GLS than the DON (reference compound), indicating that these hits can be used as GLS inhibitors to fight cancer. Bioinformatics plays an increasingly indispensable role in the integration of biological data as the diversity and complexity of data types have grown. For mining biologically active substances, there are now two types of data that are useful and accessible. The first is experimental biological activity data, such as anti-tuberculosis, anti-cancer, and anti-HIV biological activity data, which may be found in the PubChem database [27]. The second is curated data from multiple TCM databases regarding TCM plants and their derived compounds. Using bioinformatics to mine candidate molecules with diverse actions, the two forms of data present a new possibility [28]. TCM databases using the bioinformatics approaches have been used to identify inhibitors against many therapeutically important biological targets [17,[29][30][31].
The root mean square deviation (RMSD) is a protein stability metric; the lower the deviations, the more stable the protein structure [32]. GLS-DON (control), GLS-ZINC32296657, and GLS-ZINC03978829 had RMSD average values of 0.25, 0.33, and 0.31 nm, respectively. The RMSD plot revealed that GLU-ZINC03978829 binding increased GLS stability and resulted in fewer structural aberrations from its normal conformation. The bound structure of the GLS-ZINC32296657 complex is highly deviated (Figure 4a); it showed that the catalytic pocket of GLS enzyme was not forming a strong interaction with ZINC32296657, therefore it showed high deviation. In addition, the ligand RMSD also showed that DON and ZINC03978829 bind better than ZINC32296657 and are more stable (Figure 4b). Initially, ZINC32296657 formed a strong interaction with GLS until 22 ns, and after that the interaction started to lessen due to unstable bonding with the enzyme. The GLS-ZINC32296657 and GLS-ZINC03978829 backbones displayed consistent fluctuations in the GLS catalytic pocket, most probably due to differing orientations, with the biggest fluctuation regions found between 300-325 and 370-400 residues (Figure 4c). The vibrations around the equilibrium are not random; rather, they are determined by the flexibility of the local structure. The average fluctuation of all residues throughout the simulation, as well as the root mean square fluctuation (RMSF) of GLS during binding of DON (control), ZINC32296657, and ZINC03978829, were plotted as a function of GLS residue numbers. The RMSF plot revealed that GLS has residual variants in various protein domain areas. Due to their close interaction with the GLS enzyme, DON and ZINC03978829 have been shown The GLS-ZINC32296657 and GLS-ZINC03978829 backbones displayed consistent fluctuations in the GLS catalytic pocket, most probably due to differing orientations, with the biggest fluctuation regions found between 300-325 and 370-400 residues (Figure 4c). The vibrations around the equilibrium are not random; rather, they are determined by the flexibility of the local structure. The average fluctuation of all residues throughout the simulation, as well as the root mean square fluctuation (RMSF) of GLS during binding of DON (control), ZINC32296657, and ZINC03978829, were plotted as a function of GLS residue numbers. The RMSF plot revealed that GLS has residual variants in various protein domain areas. Due to their close interaction with the GLS enzyme, DON and ZINC03978829 have been shown to reduce the residual fluctuations of the enzyme. Further, the volume and density plot revealed no significant differences between the hits and control compound (Figure 4d).
A protein's solvent-accessible surface area (SASA) is the area of its surface involved in the interaction with its solvent molecules. Average SASA values for GLS-DON (control), GLS-ZINC32296657, and GLS-ZINC03978829 complexes were recorded during the 100 ns MD trajectory. The SASA values for the GLS-DON (control), GLS-ZINC32296657, and GLS-ZINC03978829 complexes were 120.02, 134.01, and 125.30 nm 2 , respectively (Figure 5a). SASA analysis showed that upon binding of DON and ZINC03978829, surface exposure was reduced more than ZINC32296657.
to reduce the residual fluctuations of the enzyme. Further, the volume and density plot revealed no significant differences between the hits and control compound (Figure 4d).
A protein's solvent-accessible surface area (SASA) is the area of its surface involved in the interaction with its solvent molecules. Average SASA values for GLS-DON (control), GLS-ZINC32296657, and GLS-ZINC03978829 complexes were recorded during the 100 ns MD trajectory. The SASA values for the GLS-DON (control), GLS-ZINC32296657, and GLS-ZINC03978829 complexes were 120.02, 134.01, and 125.30 nm 2 , respectively (Figure 5a). SASA analysis showed that upon binding of DON and ZINC03978829, surface exposure was reduced more than ZINC32296657.
The hydrogen bond is vital in the stability of the ligand-target complex [33]. Between protein and ligand, hydrogen bonds were formed within 0.35 nm. The stability of docked complexes was tested using 100 ns simulations of GLS-DON (control), GLS-ZINC32296657, and GLS-ZINC03978829 in a solvent environment. GLS-DON and GLS-ZINC03978829 form 2-5 hydrogen bonds with the GLS catalytic pocket, whereas GLS-ZINC32296657 forms 1-2 hydrogen bonds with the GLS catalytic pocket. GLS-ZINC32296657 formed fewer Hbonds due to instability in the pocket of the GLS enzyme (Figure 5b). The energy landscape of the complexes GLS-ZINC32296657, GLS-ZINC03978829, and GLS-DON has been plotted. The blue color denotes the location with the least amount of energy. The GLS-ZINC03978829 complex has two distinct global energy minima basins (in blue), whereas the GLS-DON and GLS-ZINC32296657 complexes have comparable energy minima states. The blue spots represent greater stability, whereas more blue areas represent shifts in the protein structure followed by the thermodynamically more favorable state ( Figure 6). The hydrogen bond is vital in the stability of the ligand-target complex [33]. Between protein and ligand, hydrogen bonds were formed within 0.35 nm. The stability of docked complexes was tested using 100 ns simulations of GLS-DON (control), GLS-ZINC32296657, and GLS-ZINC03978829 in a solvent environment. GLS-DON and GLS-ZINC03978829 form 2-5 hydrogen bonds with the GLS catalytic pocket, whereas GLS-ZINC32296657 forms 1-2 hydrogen bonds with the GLS catalytic pocket. GLS-ZINC32296657 formed fewer H-bonds due to instability in the pocket of the GLS enzyme (Figure 5b).
The energy landscape of the complexes GLS-ZINC32296657, GLS-ZINC03978829, and GLS-DON has been plotted. The blue color denotes the location with the least amount of energy. The GLS-ZINC03978829 complex has two distinct global energy minima basins (in blue), whereas the GLS-DON and GLS-ZINC32296657 complexes have comparable energy minima states. The blue spots represent greater stability, whereas more blue areas represent shifts in the protein structure followed by the thermodynamically more favorable state ( Figure 6).

Conclusions
Considering GLS as a viable therapeutic target because of its role in cancer development, targeting it with the selected hits is an interesting approach. The in silico technique utilized in this work might be effective in identifying possible hits from natural compounds as potent GLS inhibitors. MD simulations indicated that the hit compounds (ZINC03978829 and ZINC32296657) strongly bind to GLS and are stable. This study implies that these hits from the TCM database could be used as GLS inhibitors to fight cancer. It is noteworthy that the BE estimates and MD simulations only illustrate the interaction efficiency and stability of hits with GLS; further experimental studies are a prerequisite to optimize these hits (ZINC03978829 and ZINC32296657) as GLS inhibitors.

Conclusions
Considering GLS as a viable therapeutic target because of its role in cancer development, targeting it with the selected hits is an interesting approach. The in silico technique utilized in this work might be effective in identifying possible hits from natural compounds as potent GLS inhibitors. MD simulations indicated that the hit compounds (ZINC03978829 and ZINC32296657) strongly bind to GLS and are stable. This study implies that these hits from the TCM database could be used as GLS inhibitors to fight cancer. It is noteworthy that the BE estimates and MD simulations only illustrate the interaction efficiency and stability of hits with GLS; further experimental studies are a prerequisite to optimize these hits (ZINC03978829 and ZINC32296657) as GLS inhibitors.