Probing the anticancer mechanism of prospective herbal drug Withaferin A on mammals: a case study on human and bovine proteasomes

Background The UPP (ubiquitin proteasome pathway) is the major proteolytic system in the cytosol and nucleus of all eukaryotic cells which regulates cellular events, including mitotis, differentiation, signal transduction, apoptosis, and inflammation. UPP controls activation of the transcriptional factor NF-κB (nuclear factor κB), which is a regulatory protein playing central role in a variety of cellular processes including immune and inflammatory responses, apoptosis, and cellular proliferation. Since the primary interaction of proteasomes occurs with endogenous proteins, the signalling action of transcription factor NF-κB can be blocked by inhibition of proteasomes. A great variety of natural and synthetic chemical compounds classified as peptide aldehydes, peptide boronates, nonpeptide inhibitors, peptide vinyl sulfones and epoxyketones are now widely used as research tools for probing their potential to inhibit proteolytic activities of different proteasomes and to investigate the underlying inhibition mechanisms. The present work reports a bio-computational study carried out with the aim of exploring the proteasome inhibition capability of WA (withaferin A), a steroidal lactone, by understanding the binding mode of WA as a ligand into the mammalian proteasomes (X-ray crystal structure of Bos taurus 20S proteasome and multiple template homology modelled structure of 20S proteasome of Homo sapiens) using molecular docking and molecular dynamics simulation studies. Results One possible mode of action which is proposed here for WA to act as a proteasome inhibitor is by suppression of the proteolytic activity which depends on the N-terminal threonine (Thr1) residue hydroxyl group. Docking studies carried out with herbal ligand WA into the structures of bovine and human proteasomes substantiate that WA has the ability to inhibit activity of mammalian 20S proteasomes by blocking the nucleophilic function of N-terminal Thr1. Results from molecular dynamics simulations in water show that the trajectories of both the native human 20S proteasome and the proteasome complexed with WA are stable over a considerably long time period of 4 ns suggesting the dynamic structural stability of human 20S proteasome/WA complex. Conclusions Inhibition of proteasomal activity are promising ways to retard or block degradation of specific proteins to correct diverse pathologies. Though quite a number of selective and efficient proteasomal inhibitors exist nowadays, their toxic side effects limit their potential in possible disease treatment. Thus there is an indispensable need for exploration of novel natural products as antitumor drug candidates. The present work supports the mammalian proteasomes inhibiting activity of WA along with elucidation of its possible mode of action. Since WA is a small herbal molecule, it is expected to provide one of the modest modes of inhibition along with added favours of ease in oral administration and decreased immunogenicity. The molecular docking results suggest that WA can inhibit the mammalian proteasomes irreversibly and with a high rate through acylation of the N-terminal Thr1 of the β-5 subunit.


Background
Ubiquitin is a small 76 amino acid protein conserved in all eukaryotic cells with a molecular weight of 8.6 kDa. When polyubiquitin is attached to target proteins, tagged proteins are selected for destruction by cytoplasmic organelles called proteasomes [1]. The mammalian 20S proteasome is characterized by a cylindrical shaped quaternary structure consisting of four heptameric stacked rings, α7β7β7α7, with 7 distinct α-type and 7 distinct β -type subunits, and with C2 symmetry similar to those of the yeast [2]. Within this multienzymatic proteasome system, a proteolytic multifunctional complex 26S proteasome is involved, which consists of a 19S regulatory particle and a 20S core particle [3,4]. The two outer α rings complex with the two 19S regulatory particles, forming a narrow channel through which only denatured proteins can pass [5]. The catalytic chamber is formed by the two inner rings, each of which contains three well-characterized proteolytic activities. In particular, three active subunits β1, β2, and β5 are responsible for the three major peptidase activities: the peptidylglutamil-hydrolase like, trypsin-like and chymotrypsin-like (ChT-L) activities respectively, as investigated by mutational and crystallographic studies [6,7]. The core particle is responsible for degradation of the proteins in a progressive manner, generating peptides of 3-25 amino acids in length [8].
The enzymatic activity of β-subunits is associated with the N-terminal threonine residues, which act as nucleophiles in the hydrolysis reaction catalyzing the cleavage of peptides through nucleophilic attack. Thus these proteasomes are classified as members of the Ntn (N-terminal nucleophilic) hydrolases group [1,9]. In the eukaryotic proteasome, out of the seven different β-subunit precursors which are processed during particle maturation by autolysis, only β1, β2 and β5 subunits can be activated by the autolytic process, with the release of the amino-terminal Thr1 functioning as the nucleophile [10][11][12].
UPP is the major proteolytic system in the cytosol and nucleus of all eukaryotic cells [13,14] which regulates cellular events including mitotis, differentiation, signal transduction, apoptosis, and inflammation [15]. UPP controls activation of the transcriptional factor NF-B (nuclear factor B), which is a regulatory protein playing central role in a variety of cellular processes, including immune and inflammatory responses, apoptosis, and cellular proliferation. Since the primary interaction of proteasomes occurs with endogenous proteins, the signalling action of transcription factor NF-B can be blocked by inhibition of proteasomes, thus inhibiting the completion of the cell cycle and mitotic proliferation of cancerous cells and ultimately leading to cell death. It has been suggested that proteasomal activity is essential for tumour cell proliferation and development of drug resistance. Therefore, the development of specific inhibitors of proteasome mediated degradation pathway is now of considerable interest in the drug discovery research for cancer therapy and prevention [16]. A great variety of natural and synthetic chemical compounds classified as peptide aldehydes, peptide boronates, nonpeptide inhibitors, peptide vinyl sulfones, and epoxyketones are now widely used as research tools for studying their ability to inhibit proteolytic activity of various proteasomes from diverse origins.
Currently two proteasome inhibitors, bortezomib and NPI-0052, have been stated in clinical trials [17,18]. However, undesirable side effects such as fatigue, nausea, vomiting, peripheral neuropathy, anaemia, diarrhoea, and constipation have also been reported for these drugs [19]. Therefore, there has been an intensive drive to develop new proteasome inhibitors especially those of natural origin having little or no side effects.
Naturally occurring inhibitors fall mainly in three groups being α'β'-epoxyketones, β-lactones and TMC-95s. WA, a principle constituent of the plant Withania somnifera, has received much attention in recent years owing to its various pharmacological properties like anti-inflammatory [20], antitumor [21], antibacterial [22], antioxidant [23], anticonvulsive [24,25] and immunosuppressive properties [26]. Most recently, it was shown to potentiate apoptosis of tumor cells by suppression of NF-B activation [27][28][29]. Targeting of UPP has been identified as one of the mechanisms of WA activity exerting two distinct pharmacological activities; antitumor and anti-inflammatory [30]. Since proteasomes are required for nuclear translocation of p65/NF-B which in turn results in activation of NF-B, it is worth substantial to consider proteasomes as the target of WA. WA belongs to a family of steroidal lactones having withanolide skeleton as their basic structure ( Figure  1A). It has been reported that C 1 and C 24 of WA are highly susceptible towards a nucleophilic attack [31]. As is evident from the structure of WA ( Figure 1B) that it contains a lactone ring enclosed ester group, two conjugated ketone bonds and a three membered epoxy ring, all of which are quite susceptible to a nucleophilic attack. It has been hypothesized that WA can be a potent proteasome inhibitor and the mode of its action can be irreversible covalent modification. There is an evidence rationalising the proteasome inhibitory action of WA in which WA is shown to inhibit chymotrypsin like activity of a purified rabbit 20S proteasome (IC 50 =4.5µM) and 26S proteasome in human prostrate cancer cultures (at 5-10µM) and in xenografts (4-8 mg/ Kg/day) [31]. Herein we report the ability of naturally occurring drug candidate WA and its mode of action for binding to mammalian 20S proteasomes as macromolecular receptors using computational approaches by analyzing the interactions between WA and the proteasomes.

Methods
Comparative protein structure multiple template homology modelling of Human 20S proteasome The amino acid sequence of the β-5 subunit of h20S (human 20S proteasome) (GenBank: AAH57840.1) comprising of 263 amino acid residues was retrieved from NCBI. Since the crystal structure of h20S is not available, we used Modeller9v7 [32] to perform the homology modelling. The best possible templates were obtained using the build_profile python script. We built three dimensional structure of h20S using multiple template comparative homology modelling based on the following X-ray crystal structures obtained from PDB (Protein Data Bank) [33]: b20S (bovine 20S proteasome) at 2.75 Å resolution (PDB:1IRU), yeast 20S proteasome at 2.4 Å resolution (PDB:1RYP) and 20S proteasome from Archaeoglobus fulgidus (PDB:1J2Q). These template proteins were chosen based on a significant sequence similarity of h20S with these proteins in addition to their satisfactory crystallographic resolution. The alignment files of the targets and templates were then prepared using the align2D_mult python script (Sequence alignment module in MODEL-LER). These alignment files along with the x-ray crystal structures of the templates were used to generate the three dimensional structure models using the model_mult python script of MODELLER.

Energy minimization
Discovery Studio (Version 1.7, Accelrys Software Inc.) was then used to execute energy minimization and to perform stereochemical quality checks to arrive at the best possible three dimensional structure of the protein.
The force field applied was CHARMm and the energy minimization algorithm used was Conjugate Gradient with an RMS gradient of 0.1 using a maximum of 2000 steps. This resulted in model structures with considerably favourable potential energies. Furthermore, the variability among the models was used to evaluate the reliability of the modelling. The qualities of these models were analyzed by PROCHECKv3.4 [34].

Binding pocket analysis
Binding Site analysis module of Discovery Studio was used to identify the putative binding pockets and protein ligand binding sites in the energy minimized threedimensional structures of b20S and h20S.

Ligand Docking
The energy minimized crystal structure of Bos taurus b20S (PDB: 1IRU) and the modelled structure of h20S were used to carry out molecular dockings. For ease in dockings, the subunits which are not directly interacting with the β-5 domain of the protein were removed of b20s crystal structure, and only the domains K, L, M, X & Y were retained for further analysis. The ligand molecule Withaferin-A [PubChem:265237] was retrieved from NCBI-PubChem Compound database [35]. AutoDock 4.0 suite was used as molecular-docking tool in order to carry out the docking simulations [36]. AutoDock 4.0 was launched in a Cygwin interface in the Windows operating system. Docking logs were analyzed in the graphical user interface of ADT (Auto dock Tools) [37]. Water molecules were cleaned off from the protein crystal structure before docking. H-atoms were added to these target proteins for correct ionization and tautomeric states of amino acid residues and the non-polar hydrogens were then merged up. Kollman united atom charges and solvation parameters were assigned to the proteins. Gasteiger charge was assigned to the ligand and then nonpolar hydrogens were merged. Rigid roots were also assigned to the ligand and five bonds were made "active" or rotatable. The modified structures so obtained: x-ray crystal structure of bovine, modelled three dimensional structure of human 20S proteasome, and the structure of ligand WA accounting the flexibility of its bonds, were converted to PDBQT format in ADT, as required in AutoDock calculations. The Lamarckian Genetic Algorithm was used with a population size of 150 dockings. Five million energy evaluations were used in the docking experiments. All other parameters, e.g. crossover rate and mutation rate, were run with default settings. The grid size for specifying the search space was set at 40 × 30 × 30 centered on Thr1 of the β5 subunit with a default grid point spacing of 0.375 Å. Energy scoring function of AutoDock 4 is based upon the calculation of pairwise atomic terms including evaluations for different secondary interactions, dispersion/repulsion, hydrogen bonding, electrostatics, and desolvation [38]. Pre-calculated grid maps, which store grids of interaction energy based on the interaction of the ligand atom probes with receptor target, were obtained using AutoGrid. The user defined three dimensional grid must surround the region of interest in the macromolecule, and the ligand was limited to this search space during docking. The results are clustered into bins of similar conformations according to the cluster root mean square deviation (rmsd) and orientation.

Confirmation of the docking results
The docking results obtained using AutoDock were also confirmed using ParDOCK [39] , which is an all atom energy based monte carlo docking protocol. Docking using ParDOCK requires a reference complex (target protein bound to a reference ligand) and a candidate molecule along with specific mention of the centre of mass of the cavity on which the ligand is to be docked.

Molecular Dynamics simulations of human proteasome in water
The AMBER v.10 package [40] was used to prepare the protein and the ligand files as well as for the MD (Molecular Dynamics) simulations. The binding complex of h20S/WA obtained using ParDOCK and the free protein simulated in this study were neutralized by adding appropriate number of chloride counterions and were solvated in a octahedron box of TIP4P water with a 10 Å distance between the protein surface and the box boundary [41]. The partial atomic charges for the ligand were obtained using "antechamber" [42] module of Amber. The energy minimization and MD simulations of h20S and its complex with WA were carried out with the aid of the PMEMD module of the AMBER 10 program. First of all, the simulated binding complex was effected with a 2500 step minimization using the steepest descent algorithm followed by a 1000 step minimization using conjugate gradient to remove bad steric contacts. Topology and parameter files for the protein were generated using "ff03" and for the drug using "gaff" based on the atom types of the force field model developed by Cornell et al [43]. Then the system was equilibrated beginning with the protein atom restrained simulations having 150 ps equilibration dynamics of the solvent molecules at 300 K and a harmonic potential with a 10 kcal/mol restraint force. Next step involved the equilibration of the solute molecules with a fixed configuration of the solvent molecules in which the system was slowly heated from T = 10 to 300 K in 58 small intervals of 2.5 ps each for a total period of 145 ps. The entire system was then equilibrated at 300 K for 100 ps before a sufficiently long MD simulation (4 ns) at room temperature. The MD simulations were performed with a periodic boundary condition in the NPT ensemble at T=298.15 K with Berendsen temperature coupling [44] and constant pressure P=1 atm with isotropic molecule-based scaling . The SHAKE algorithm [45] was applied to fix all covalent bonds containing hydrogen atoms. We used a time step of 2 fs and a nonbond-interaction cut-off radius of 10 A°. The Particle Mesh Ewald (PME) method [46] was used to treat longrange electrostatic interactions. The coordinates of the trajectory was sampled every 1 ps for analysis of the energy stabilization and RMSD values of the protein as well as that of the complex. MD simulations were performed on a 320 processors SUN Microsystems clusters at Supercomputing Facility (SCFBio) at Indian Institute of Technology Delhi.

Results and discussion
Docking of WA into b20S proteasome One possible mode of action which is proposed here for WA to act as a proteasome inhibitor is by suppression of the proteolytic activity which depends on the N-terminal threonine (Thr1) residue hydroxyl group, which is responsible for catalyzing the cleavage of peptides through nucleophilic attack. Using binding pocket analysis, S1 pocket of β-5 subunit was obtained as one of the putative binding site. As evident from the docking of WA into b20S (Figure 2A), WA is trapped inside this protein pocket. Figure 2B shows the ligand occupying the S1 cavity of the receptor being represented as a mesh surface. As AutoDock reports the best docking solution for each GA run and also performs a cluster analysis in which the total number of clusters and the rank of each docking mode (cluster rank) is reported, in 7 out of 10 docked conformations obtained by the clustering analysis at 2.0 Å, the carbonyl group of the lactone ring is found closest to the hydroxyl group of Thr1 ( Figure 3A). The various properties of the docked conformation are shown in Table  1. The binding energies of the conformations of this cluster range from -6.37 to -6.11 Kcal/mol. The highest binding energy of -7.62 Kcal/mol was obtained for a conformation in which the epoxy group of the ligand is close to the nucleophilic hydroxyl group of the protein but with only a 10% clustering frequency ( Figure  3B).

Homology modelling of human 20S proteasome
The three dimensional structure of human proteasome was determined by comparative homology modelling with satisfaction of spacial constraints using multiple known X-ray crystal structures as templates. The chosen templates showed significant similarities to the h20S with e-values equal to 0. Five modelled three dimensional structures of human proteasome were obtained using Modeller out of which the model having the least DOPE score ( Table 2) was chosen for the purpose of studying ligand and protein interactions. The quality and reliability of the model was ensured by assessing the backbone and side-chain conformations, bond lengths, angles, and residue contacts of the model through ProCheck, magnitudes of which are well within the criteria established for reliable structures (data not shown). The model was almost as good quality as those of the reference templates as evident from the results obtained using Ramachandran plot analysis (data not shown) for comparison of stereochemical and energetic properties of the models with those of the templates. The first 59 amino acids of the protein were then removed off from the structure as these are a part of a propeptide which is absent in the mature form. This model was energy minimized with an energy lowering of around 19,000 Kcal/mol and this energy minimized structure was used further for docking analysis.

Docking of WA into modelled h20S proteasome
Binding energy of -6.92 Kcal/mol was obtained from docking of WA into homology modelled h20S. The various properties listed in Table 1 provide sufficient results in order to support the ongoing mechanism of inhibition of h20S. Docked withefrin A positions itself into the S1 pocket of the receptor as shown in Figure 4. Moreover the ligand occupies the same conformation as required to facilitate the nucleophilic attack, positioning   its lactone ring in vicinity of the hydroxyl nucleophile ( Figure 5), with a clustering frequency of 30%. It has been shown using kinetic analysis and X-ray diffraction studies that the ester bond of specific proteasome inhibitor lactacystin covalently modifies the Nterminal threonine of the β5-subunit, which is critical for proteasome inhibition [1,[47][48][49]. Thus it is quite probable that lactacystin-like reaction occurs with WA also as it contains an internal ester bond in a δ-valero lactone ring ( Figure 1B). A similar kind of cleavage of the lactone ring by serine protease has been reported in which the 3-benzyl-2-oxetanone, a β-lactone has been found to be a slowly hydrolyzed substrate of α-chymotrypsin [50]. Other herbal ligands like polyphenols [51] especially green tea polyphenols like EGCG and its analogs, genistein etc. [52,53], which have an ester bond susceptible to nucleophilic attack by Thr1 have also been reported to possess proteasome inhibition activity. Our results obtained from docking of WA into bovine and human 20S proteasome structures substantiate the proposed inhibition mechanism.

MD simulations in water
The h20S/WA protein-drug binding complex with the binding energy of -6.91 kcal/mol obtained using Par-DOCK ( Figure 6) was used for carrying out MD simulations. After MD simulations, we calculated RMSDs between Cα trajectory of h20S and Cα of its modelled structure recorded every 1 ps. The RMSDs for the trajectory of h20S complexed with WA were also calculated using its initial model as a reference structure. The results in Figure 7A show that the RMSDs of the trajectory of the complex were always less than 2 Å for the entire simulation suggesting the stability of our simulation system. The trajectories were not greatly different from the modelled structure, with only minor movements of the Cα of the protein observed. The adherence of the total energy trajectories to more or less constant values for both the complex and the protein were seen    during the entire simulation length ( Figure 7B), with the energy values of the complex much lowered than that of the native protein indicating thermodynamic stability of the complex. The simulation length used in this study was long enough to allow rearrangement of side chains of the native as well as the drug complexed protein to find their most stable binding mode. Thus the present MD simulations along with the molecular docking experiments made clear the dynamic structural stability of h20S in complex with the drug WA, together with the inhibitory mechanism.

Conclusions
Since proteasomes play an essential role in the turnover of cellular proteins, modulation or inhibition of proteasomal activity are thus promising ways to retard or block degradation of specific proteins in order to correct diverse pathologies. Though nowadays there exist quite a number of selective and efficient proteasomal inhibitors, the toxic side effects of these compounds strongly limit their potential in possible disease treatment. Thus there is an indispensable need for exploration of novel natural products as anti-cancer drug candidates. The study conducted here makes use of molecular docking and molecular dynamics simulation approaches, which include the search in space for the energetically most favorable conformation of a protein-ligand complex and the scoring of the resulting geometries with respect to binding energy, to analyze the proteasome inhibitory potential of WA and to investigate the underlying inhibitory mechanism. We have obtained significant results delineating the mammalian proteasomes' inhibitory activity of WA alongwith elucidation of its possible mode of action. Since WA is a small herbal molecule, it is expected to provide one of the modest modes of inhibition alogwith added favors of ease in oral administration and decreased immunogenicity. Conclusively it is strongly suggested here that WA is a potent proteasome inhibitor and should be looked forward for further clinical investigations as a possible proteasome inhibitory drug candidate.