1. INTRODUCTION
Glutamate cysteine ligase (GCL) is an enzyme that catalyzes the initial and rate-limiting step of glutathione (GSH) biosynthesis [1,2]. In the synthesis of the tripeptide GSH, the ATP-dependent mechanism proceeds via a γ-glutamylphosphate intermediate [2-4] with a subsequent nucleophilic attack by the α-amino group of l-cysteine to produce the dipeptide γ-glutamylcysteine [1,2]. Glutathione synthetase couples the resulting γ-glutamylcysteine to l-glycine to generate reduced GSH, an abundant cellular reducing agent [1]. GSH is the major low molecular weight cellular thiol and reduced GSH is present in most cell types at millimolar levels, whereas the oxidized form glutathione disulfide is less abundant. GSH plays roles in many cellular functions, such as maintenance of reduced protein thiols, detoxification of hydrogen peroxide and lipid peroxides, secondary metabolism and non-enzymatic scavenging of free radicals [5]. GSH also protects against apoptotic cell death following exposure to antineoplastic agents, radiation and receptorbased death signals [6-10].
GCL activity is modulated by free l-cysteine availability [11], transcriptional regulation [12] and post-translational modifications [13]. Because of its central role in GSH homeostasis, GCL is an attractive target for antitumor drug development. Upregulation of GCL catalytic subunit (GCLC) mRNA and GCL activity have been frequently observed in cells derived from human tumors resistant to chemotherapeutic agents [14-16]. Increased production of GSH protects against reactive oxygen and nitrogen species [5,17] and facilitates detoxification of electrophilic xenobiotics by the GSH S-transferases [18]. It has been reported that drug resistance in tumors can be overcome by the administration of the GCL inhibitor L-buthionine-S,R-sulfoximine (BSO) [19], which subsequently depletes GSH and sensitizes the tumor cells to radiation treatment and chemotherapy.
Our previous study revealed that GSH could be an important factor for the selective toxicity toward central nervous system-derived tumor cells [20], which suggests that targeting the key enzyme for GSH homeostasis, namely GCL, with highly selective inhibitors could be utilized for controlling the development and progression of cancer. GCL is a heterodimeric enzyme with a 73-kDa catalytic subunit and a 31-kDa modifier subunit, and GCLC contains the glutamate, cysteine and ATP binding sites and has catalytic activity even as a monomer [1,21]. Thus, structural analysis of GCLC with its possible inhibitory ligands could be of importance for successful antitumor drug development. Although a few GCLC models have been publicized [22,23], no human GCLC (hGCLC) model has been reported to the best of our knowledge. Molecular modeling has found widespread utility in the field of drug development [24-26], and in the present study we will report the homology modeling and structural analysis of hGCLC by a highly sophisticated software package, the Molecular Operating Environment (MOE) 2010.10 (Chemical Computing Group Inc., Montreal, Canada).
2. COMPUTATIONAL METHODS
2.1. Homology Modeling of hGCLC
Homology modeling of hGCLC was executed as previously reported [27]. In brief, the hGCLC (NCBI reference sequence: NP_001489.1) [28] sequences and the crystal structure coordinates of yeast GCLC (yGCLC; PDB code: 3LVV) [23] were loaded into the MOE. The primary structures of hGCLC and yGCLC were aligned, carefully checked to avoid deletions or insertions in conserved regions and corrected wherever necessary. A series of hGCLC models were independently constructed with the MOE using a Boltzmann-weighted randomized procedure [29] combined with specialized logic for the handling of sequence insertions and deletions [30]. The models with the best packing quality function were selected for full energy minimization and further inspection.
2.2. Assessment of the Modeled Structure
The qualities of the protein folds of the hGCLC homology model were evaluated with the MOE by calculating the effective atomic contact energies, comprising the desolvation free energies required to transfer atoms from water to the interior of the protein [31]. Briefly, the contact desolvation energies were calculated for 18 different atom types of the 20 common amino acids that were resolved based on the clustering pattern of their properties. The contact potentials for each atom type were measured within a contact range of 6 Å by explicitly accounting for neighboring interactions. The overall geometric and stereochemical qualities of the final modeled structure of hGCLC were examined using Ramachandran plots generated within the MOE [32,33].
2.3. Molecular Electrostatic Potential (MEP) Mapping
Electrostatic potential surfaces were calculated by solving the nonlinear Poisson-Boltzmann equation using a finite difference method as implemented in the MOE. The molecular electrostatic interactions form a crucial part of the non-covalent interaction energy between the molecules. The MEP on a molecular surface can be used to visually compare different molecules, analyze docking studies and identify sites that interact with ligands. For example, the surface MEP was utilized to relate a nucleotide mutation with the potential values [34]. In the present study, the MEP was colored in deep blue to indicate the most positive potential and in deep red to represent the most negative potential.
2.4. Binding Site Selection and Exploration
The binding site selection and exploration for hGCLC was executed as previously reported [27]. In brief, the Site Finder module of the MOE was used to identify possible substrate-binding pockets within the newly generated 3D structures of hGCLC. Hydrophobic or hydrophilic alpha spheres served as probes denoting zones of tight atom packing. These alpha spheres were utilized to define potential ligand-binding sites (LBSs) and as centroids for the creation of dummy ligand atoms [35,36]. The dummy atoms were matched to the corresponding alpha spheres during the characterization of the LBSs in hGCLC. This method generates bound conformations that approach crystallographic resolutions [37].
2.5. Alpha Sphere and Excluded Volume-Based Ligand-Protein Docking (ASE-Dock)
The docking and analysis of the ligand-protein interaction between GSH (or BSO) and hGCLC were also performed with ASE-Dock in the MOE [38]. In the ASEDock module, ligand atoms have alpha spheres within 1 Å. Based on this property, concave models are created and ligand atoms from a large number of conformations generated by superimposition with these points can be evaluated and scored by maximum overlap with alpha spheres and minimum overlap with the receptor atoms. The scoring function used by ASE-Dock is based on ligand-protein interaction energies and the score is expressed as a Utotal value. The ligand conformations were subjected to energy minimization using the MMF94S force field [39]. From the resulting 500,000 poses, the 200 poses with the lowest Utotal values were selected for further optimization with the MMF94S force field. During the refinement step, the ligand was free to move within the binding pocket.
3. RESULTS AND DISCUSSION
3.1. Homology Modeling of the hGCLC Structure
The sequence alignment of GLCL is shown in Figure 1. The alignment revealed that the critical active site residues Glu52 and Glu103 [40] were conserved in yGCLC and hGCLC. yGCLC (PDB code: 3LVV) was selected as a template (Figure 2(a)) for the present 3D structure modeling of hGCLC because of its good crystal structure resolution (2.2 Å) and its information was the latest (from 2010) [23]. The % sequence identity between yGCLC and hGCLC was 39.7%. For the construction of the hGCLC model, 100 independent models of the target protein were built using a Boltzmann-weighted randomized modeling procedure in the MOE that was adapted from reports by Levitt [29] and Fechteler et al. [30]. The intermediate models were evaluated by a residue packing quality function, which is sensitive to the degrees to which non-polar side-chain groups are buried and hydrogen bonding opportunities are satisfied. The hGCLC model with the best packing quality function and full energy minimization was utilized in the present study (Figure 2(b)). The relatively few β-sheets in the hGCLC model is probably due to the differences in the residues between yGCLC and hGCLC, such as Thr91-Arg94 and Glu137-Cys142 (in hGCLC) that could not form hydrogen bonds with adjacent structures to create β-sheets.
3.2. Analysis of the Contact Energies for the hGCLC Model
As reported by Zhang et al. [31], the effective atomic contact energies were calculated using the MOE for heavy atoms of standard amino acids within a contact range of 6 Å, assigning energy terms in kcal/mol for each contact pair. These energies were summed for each resi
Figure 1. Evaluation of the 2D structures of the hGCLC model. Homology-aligned sequences of yGCLC (PDB code: 3LVV; green) and hGCLC (magenta). Red line: α-helix; Blue line: turn; Yellow line: β-sheet.
(a)
(b)
Figure 2. The secondary structures of the GCLC models. Upper panel: yGCLC (PDB code: 3LVV). Lower panel: hGCLC.
due, and in general, a large negative value indicated that the residue was predominantly in contact with hydrophobic atoms and therefore expected to be in a buried protein environment. Conversely, residues with positive energy terms indicated contacts with predominantly hydrophilic atoms, and were expected to be in more solvent-exposed regions of the proteins. The contact energy profiles of the homology-modeled hGCLC (Figure 3(b)) were compared with those of the X-ray structure of yGCLC (Figure 3(a)). Although 5 out of 36 residues at the LBSs had differently charged contact energies, such as His94 (yGCLC; –3.4 kcal/mol) and Arg94 (hGCLC; 3.3 kcal/mol), most residues had quite similar contact ener gies, such as Phe271 (yGCLC; –19.8 kcal/mol) and Phe254 (hGCLC; –20.1 kcal/mol), and the trends of the variation in the contact energy in most parts of the hGCLC model
Figure 3. Contact energy profiles of the GCLC models. (a) yGCLC (PDB code: 3LVV); (b) The constructed hGCLC model. The positions of the amino acid residues are shown on the x-axis, while the contact energies are shown on the y-axis.
were in good agreement with those of the X-ray structure of yGCLC.
3.3. Evaluation of the Stereochemical Qualities of the hGCLC Model
The phi and psi backbone dihedral angles for hGCLC were scored using 2D probability distributions calculated on a high-resolution collection of X-ray structures containing approximately 100,000 residues from 500 protein structures [41]. Each probability distribution was estimated with 2-degree spacing for each of the phi and psi backbone dihedral angles with separate histograms for preproline, proline, glycine and general amino acids. The stereochemical qualities of the hGCLC model were assessed by Ramachandran plots (Figure 4). 86.2% of the residues were in the favored region, 12.4% were in the allowed region and only 1.4% were in the disfavored region. These results indicate that the phi and psi backbone dihedral angles in the hGCLC model are reasonably accurate.
3.4. Structural Comparisons of the GCLC Models
The GCLC sequences were realigned and reinspected with those of other species [42] for the superimposition
Figure 4. Ramachandran plots for the hGCLS model. Green: favored region; Light-brown: allowed region.
analysis in the MOE (Figure 5). Root mean square deviation (RMSD) values between the main chain atoms of the yGCLC (3LVV) vs hGCLC after main chain fit were 1.83 Å. RMSD values for each residue were also analyzed. The RMSD values for the residues located in the LBS were about 2 Å or less (Figure 6). A superimposition of the template yGCLC (green) and hGCLC (magenta) models revealed that the GCLC models exhibited significant 3D similarities (Figure 7(a)). They also presented similar structures at their LBSs (Figure 7(b)).
3.5. The MEP Maps for the GCLC Models at LBSs
The MEP maps can play a vital role on analyses and predictions of molecular interactive behaviors and properties. For instance, they can be used to compare two molecules visually, which helps in identifying sites that act attractively on ligands by matching opposite electrostatics. Electrostatic interactions are one of the main parts of the interaction energy between ligands and receptors, and govern the strength of non-bonded interactions and
Figure 5. Realigned GCLC sequences for the superimposition analysis. Residues in the LBS are enclosed in red rectangles. Upper sequence: yGCLC (PDB code: 3LVV); Lower sequence: hGCLC.
molecular reactivity. In the case of a ligand-receptor interaction at the catalytic site, the ligand experiences a unique environment in terms of the electrostatic, steric and hydrophobic properties. Variations in these properties near the catalytic site of receptors can contribute to their selectivity/specificity [43]. The MEP maps of the GCLC models are shown in Figure 8. The hGCLC model had a positive potential at the LBS (colored in blue, Figure 8(b)), which indicates that hGCLC possibly attracts negatively charged parts of ligands at the LBS. On the other hand, the yGCLC model had a negative potential at the LBS (colored in red, Figure 8(a)). Although the sequence identity of the LBS between hGCLC and yGCLC was 86%, the MEP maps were different possibly due to the structural differences at and adjacent to the LBSs. These results suggest that binding orientations of ligands at the LBS can be different between hGCLC and yGCLC.
Figure 6. RMSD values between the main chain atoms of the yGCLC (3LVV) vs hGCLC after main chain fit. The positions of the amino acid residues are shown on the x-axis, while the RMSD values are shown on the y-axis.
(a)
(b)
Figure 7. Structural comparison of the GCLC models. (a) A superimposition of the template yGCLC (green) and hGCLC (magenta) models; (b) A superimposition of the yGCLC (green) and hGCLC (magenta) models at the LBSs.
(a)
(b)
Figure 8. The MEP maps for the GCLC models. (a) The MEP map for the yGCLC model; (b) The MEP map for the hGCLC model. Arrow: LBS. Deep blue: Most positive potential. Deep red: Most negative potential. The LBSs are enclosed in a yellow circle.
3.6. Docking Simulations of GSH to GCLC
It has been reported that GCL is feedback regulated by the end product, GSH. GSH inhibits GCL by competing with L-glutamate [44], suggesting that the two binding sites are coincident. Further studies with GSH analogues such as ophthalmic acid, S-methylglutathione, and glutathione disulfide have demonstrated that the free thiol group of GSH is necessary for maximal inhibition [1,44]. However, the precise mode of GSH binding is unknown. The ASE-Dock was performed to evaluate the present docking simulation and showed that GSH bound at the LBS but had a slightly different binding orientation between the yGCLC (Figure 9(a)) and hGCLC (Figure 9(b)) models. The similarity of the bound location at the LBS between the docked GSH-yGCLC pose and the hGCLC model suggests that the present methods are capable of generating the GSH-hGCLC model similar to
(a)
(b)
(c)
Figure 9. Docking simulations of GSH and BSO-P to the GCLC models. (a) ASE-Dock of GSH for yGCLC; (b) ASEDock of GSH for hGCLC; (c) ASE-Dock of BSO-P for hGCLC. Blue, nitrogen; gray, carbon; magenta, phosphorus; red, oxygen and yellow; sulfur.
the reported near-native GCLC complex. The results of the slightly different binding orientations of GSH between GSH-yGCLC and GSH-hGCLC complexes reflected the different MEP maps at the LBSs in the GCLC models. This also suggests that the homology modeling of hGCLC and the docking simulations in the present study were performed reasonably well. Further, a clinically available specific GCL inhibitor BSO in the phosphorylated form (BSO-P) was docked to the hGCLC model (Figure 9(c)). The bound location of BSO-P in hGCLC was similar to that of GSH at the LBSs in the yGCLC and hGCLC models, which further indicates the accuracy of the present docking simulation.
4. CONCLUSION
Examination of the hGCLC structures provides considerable insight into the catalytic mechanism of the enzyme and suggests approaches by which GCL inhibitors with greater selectivity may be attainable. The analysis of the GSH-binding region in hGCLC revealed that a subtle difference of the GSH binding orientation can be found between species. Consequently, the location of the ligand possibly influences the physico-chemical properties of the LBS in the enzyme and has some effects on the binding of the inhibitor to the enzyme. Thus, detailed analysis of the ligand-protein interaction is of great significance in designing in silico hGCLC-inhibiotor models for successful development of antitumor drugs. The main objective in the present study was to create a hGCLC model. Analyses of the structural properties of the hGCLC and the docking simulations of the GSH-hGCLC pose suggest that the present methods are capable of generating the hGCLC model similar to the nearnative yGCLC. Consequently, it is proposed that the hGCLC in the present study will be suitable for further insilico structure-based de novo drug design. Furthermore, to the best of our knowledge, this is the first report of a hGCLC model.
5. ACKNOWLEDGEMENTS
This study was partially supported by a grant-in-aid from the Promotion and Mutual Aid Corporation for Private Schools of Japan.