Exploration of the binding modes of l-asparaginase complexed with its amino acid substrates by molecular docking, dynamics and simulation

Acute lymphocytic leukemia (ALL) is an outrageous disease worldwide. l-Asparagine (l-Asn) and l-glutamine (l-Gln) deamination plays crucial role in ALL treatment. Role of Erwinaze® (l-asparaginase from Erwinia chrysanthemi) in regulation of l-Asn and l-Gln has been confirmed by the experimental studies. Therapeutic research against ALL remained elusive with the lack of structural information on Erwinaze® enzyme. In this present study, homology model of the Erwinaze® was developed using MODELLER and the same was validated by various quality indexing tools. For the apo state enzyme and ligand bound state complexes molecular dynamics (MD) simulation was performed. The trajectory analysis showed the confirmational changes of structures in the dynamic system. Ligand binding mechanisms were studied using different docking tools to interpret the various ligand-receptor interactions and binding free energies. MD simulation of docked complex with l-Gln ligand substrate showed the defined structural folding with stable conformation over the l-Asn complex in dynamic environment. This research reports give much more information on structural and functional aspects of Erwinaze® with its ligands which may be useful in designing of effective therapeutics for ALL. Electronic supplementary material The online version of this article (doi:10.1007/s13205-016-0422-x) contains supplementary material, which is available to authorized users.


Introduction
Leukemic cell's inability to synthesize L-Asn on their own has primed L-Asparaginase enzyme as a dynamic factor of chemotherapy in acute lymphoblastic leukemia (ALL) therapy (Ohnuma et al. 1970;Wriston and Yellin 1973). In blood vascular system, L-Asparaginase hydrolyses L-Asn into L-Aspartic acid and ammonia building the leukemic cells barren of essential exogenous L-Asn, leading to protein synthesis inhibition followed by apoptosis (Lubkowski et al. 1994;Neuman and McCoy 1956;Mashburn and Wriston 1964), hence making this potent against ALL.
Nowadays, salable L-Asparaginases are predominantly produced from E. coli and Erwinia chrysanthemi (Godfrin and Bertrand 2006). Apart from direct production from various bacterial and fungal strains, several attempts were also made for the production of recombinant enzyme by cloning and expression of L-Asparaginase encoding genes from Erwinia chrysanthemi (Kotzia and Labrou 2007) and Erwinia carotovora (Kotzia and Labrou 2005) and other microbes in E. coli. It is found that against the cells of reverted ALL patients, the commercial anti-cancer enzyme has in vitro resistance (Klumper et al. 1995) along with its high glutaminase activity and low substrate specificity, causing pancreatitis, liver dysfunction, coagulation anomalies causing intracranial thrombosis or hemorrhage, neurological seizures and leucopenia (Duval et al. 2002). L-Glutaminase activity of same enzyme leading to deamination of L-Gln to L-Glutamate in blood plasma induces some toxic effects in normal cells (Capizzi and Cheng 1981;Storti and Quaglino 1970). This facilitates the necessity for novel and healthy L-Asparaginases from innocuous microorganisms with elevated substrate affinity, amended stability, low glutaminase activity, adequate halflife and least KM value under physiological conditions to overcome the above said challenges encountered in the recent scenario. Though we have sufficient data on the production, optimization of bioprocess and purification of enzyme etc., no research has been done on the molecular aspects of the enzyme from Erwinia chrysanthemi. In the absence of crystal structure it is highly difficult to get the molecular information about the enzyme like its interactions with the substrates, enzyme stability etc. Therefore, the current study is aimed towards the in silico modeling of L-Asparaginase structure from Erwinia chrysanthemi, molecular interactions with the substrates through docking, and testing the stability of the enzyme and docked complexes under physiological conditions by molecular dynamic simulations methods.

Preparation of receptor and ligands
Ligand molecules L-Asn (C 4 H 8 N 2 O 3 ) and L-Gln (C 5 H 10 N 2 O 3 ) whose molecular masses are 132.12 g/mol and 146.14 g/mol were retrieved from Zinc database with ID numbers 1532556 and 1532526, respectively. Then they were subjected to energy minimization using the MMFF (Merck Molecular Force Field) (Halgren 1996) of VLi-feMDS v 4.3 that works based on MM3 force fields until reaching global minima.
Homology modeling approach was used to investigate the tertiary structure of the enzyme. Hypothetical configuration of enzyme was obtained by MODELLER v 9.13 tool using amino acid sequence provided by Drugbank (http://www.drugbank.ca/drugs/DB08886). Further the modeled enzyme was validated using Ramachandran plot analysis by Rampage (Lovell et al. 2003), followed by determination of QMEAN 6 score (Benkert et al. 2008), DFire energy value (Zhou and Zhou 2002) using SWISS-MODEL server and ERRAT 2.0 (Colovos and Yeates 1993) tools to verify the stereochemical quality of the model by analyzing the phi (ø) and psi (w) torsion angles, estimation of local quality of the modeled enzyme, assessment of non-bonded atomic interactions and for appraising the growth of crystallographic model construction and refinement, respectively.

Molecular docking
iGEMDOCK v 2.1 is a graphical atmosphere for identifying the pharmacological interactions and virtual screening that are beneficial for pinpointing the lead compounds and understanding the mechanism of ligand binding against a therapeutic target. iGEMDOCK, a flexible docking tool can be used for the docking, virtual screening and postscreening analysis. The post analysis tools works by using K means and hierarchal clustering methods (Hsu et al. 2011).
Interactive interface was provided initially for preparing target protein's binding site and compound library screening in GEMDOCK. The complete modeled structure of Erwinaze Ò was uploaded in the ''Prepare Binding Site'' and the ''By current file'' choice was kept because the uncut surface will be checked for binding, instead of specifying a single cavity. Then the in-house docking tool GEMDOCK docks the compounds from library into receptor binding site. Protein-compound interaction profiles were generated and analyzed by post analysis tools. iGEMDOCK finally ranks and visualizes the compound based on energy based scoring function and pharmacological interactions (Kaladhar 2012). During docking process, GA parameters were set as population size of 300, generations of 80 and solutions number as 10. Stable docking (slow) was performed with the docking scoring function as GEMDOCK scoring function and 1.00 was set to ligand hydrophobic and electrostatic preferences. Low energy profile indicates the stable system and it represents the likely binding interaction.
To verify the results obtained by iGEMDOCK, again molecular docking was performed by Patch Dock server by uploading the structures to web server (Schneidman-Duhovny et al. 2005) that works based on shape complementarity principles and again the outcomes were refined with FireDock server (Andrusier et al. 2007;Mashiach et al. 2008) that reshuffles the interface side chains and amends the molecule's relative orientation. Analysis of ligand binding interactions and docking viability was established on Fire Dock scores and visualizations with Pymol. Docking parameters were set as 30 as rotation angle, 30 as number of placements and 10 as ligand wise results. Intel Core i7 7470 CPU @ 3.40 GHz of DELL origin, with 8 GB DRR RAM under Windows 8.1 OS was used to perform all the docking runs.

Molecular dynamics and simulation
MD simulations were executed for the homology modeled Erwinaze Ò enzyme, Erwinaze Ò -L-Asn (complex 1) and Erwinaze Ò -L-Gln (complex 2) docked complexes developed from molecular docking to ratify the stability for anticancer enzyme in apo state and bound state with the substrates in dynamic system. Generating the L-Asn and L-Gln topologies using PRODRG server was the early step in MD simulations (SchuÈ ttelkopf and Van Aalten 2004). After defining ligand topologies, MD simulation for modeled enzyme and docked complex were carried using GRO-MACS 4.6.5 program package under Ubuntu 14.04 operating system. Steepest algorithm using the GROMOS 96 43a1 (van Gunsteren et al. 1996;Scott et al. 1999) was used for energy minimization, dismissing when the maximum force was found lesser than 1000 kJ mol -1 nm -1 . The modeled enzyme and the docked complexes were solvated in the system of a cubic box with a size of 1.0 nm and at least 2.0 nm between any two periodic images of a protein to provide an aqueous environment. With the addition of one chlorine ion the system was neutralized and periodic boundary conditions were applied in all directions. The cubic interpolation order in particle mesh Ewald (PME) simulation method was 4 and the grid spacing for FFT (Fourier spacing) was 0.16. In the neighbor searching method, the short range neighbor list cutoff of 1 Å was taken commonly for electrostatic interactions and Vander Walls interactions. The LINCS (Hess et al. 1997) and SETTLE algorithms (Miyamoto and Kollman 1992) were applied to constrain all bond lengths and geometry of water molecules, respectively. For 100 ps duration, the two equilibration phases NVT ensemble with a constant temperature of 300 K, coupling constant of 0.1 ps and NPT ensemble with constant pressure of 1 bar, coupling constant of 2 ps were applied for all the molecules. Modified Berendsen thermostat coupling scheme algorithm was employed for both ensembles of equilibration. Once the system equilibration with constant temperature and pressure was done, to carry out the structural and energy analyses30 ns production MD run was performed. Run trajectories were obtained and with GROMACS utilities analysis was carried out. Using g_energy, g_rms, g_rmsf and g_gyrate quality assurance of all the molecules was performed. g_hbond was used to cross check the docking results in terms of hydrogen bonding pattern between the receptor and ligand substrates. Using Xmgrace tool the entire trajectory results were analyzed (Turner 2005).

Results and discussion
Molecular modeling of L-asparaginase from Erwinia chrysanthemi For predicting protein three-dimensional structure based on user provided alignment of a sequence to be modeled with known related structures, homology or comparative modeling tool MODELLER can be used. BLAST search was performed to identify the most similar sequences whose crystal structures were experimentally determined. Search results suggested 1HFJ_A, 1O7J_A, 1ZCF_A and 2JK0_A as potential templates of bacterial origin having 99, 98, 79 and 79 % of sequence identity and SuperPose (http:// wishart.biology.ualberta.ca/superpose/) gave the backbone RMSD values of 0.31, 0.37, 0.83 and 0.43 nm, respectively, from the query. Calculation of a model containing all non-hydrogen atoms was automatically done by MODELLER using all the four templates to cover the complete query sequence. Best homology model was chosen based on the low discrete optimized protein energy (DOPE) score and RMSD values amongst a total of ten models predicted by MODELLER (Š ali and Blundell 1993;Fiser et al. 2000).
The hypothetical model of Erwinaze Ò enzyme ( Fig. 1) was validated using Rampage geometric evaluations to get Ramachandran plot. The plot has 97.2, 2.8 and 0.0 %of residues in favored, allowed and disallowed regions, respectively, and this strongly supports the geometric fitness of the modeled enzyme. Qmean6 score of 0.671, DFire energy of -455.15 and 84.326 % overall quality factors from ERRAT 2.0 indicates the good resolution structure. Secondary structure of the protein was also analyzed using PredictProtein tool (http://www. predictprotein.org). Protein was rich in loop region with a percentage of 46.01 residues and the remaining 29.14 and 24.85 % ones as helix and strands with two N-glycosylation sites in it. Along with these it also has eight sites for protein kinase C phosphorylation, three sites for Casein kinase II phosphorylation and one Tyrosine kinase phosphorylation. The final structure was further used to study the molecular interactions with ligand molecule.

Docking with L-Asn and L-Gln
Two ligand substrates namely L-Asn and L-Gln were docked into the catalytic site of Erwinaze Ò . Dock runs of ligands on the enzyme were done using iGEMDOCK v 2.1 and PatchDock server. When docking runs were carried out by iGEMDOCK, Erwinaze Ò has shown very high binding affinity with L-Asn and L-Gln (-52.6373 kcal/mol and -57.2457 kcal/mol, respectively). Results obtained by iGEMDOCK were again validated by PatchDock Fig. 1 Structure of Erwinaze Ò predicted by MODELLER webserver followed by refinement of obtained results with FireDock server.
Outcomes strongly supported the previous results with a very good binding efficiency between receptor and ligands with the calculated global energy values of -20.41 kcal/mol and -25.32 kcal/mol for L-Asn and L-Gln, respectively, by FireDock server. L-Asn made inter molecular hydrogen bonds with THR-15, THR-26, THR-27 and MET-121 of Erwinaze Ò and on the other side L-Gln generated bonds with THR-165, THR-167 and ASN-180 of same enzyme (Table 1). In fact, the enzyme has more affinity towards L-Asn; as it has formed the four inter molecular hydrogen bonds with the ligand atoms with the distances of 2.6, 2.5, 3.4 and 2.7 Å which are lesser compared to the bond lengths resulted with L-Gln.
Summary of the binding energies by all the docking methods were described in Table 2. The mode of binding and cavity of ligands with receptor was visualized in PyMol molecular graphics viewer (Fig. 2). Cross checking of the modeled enzyme binding sites in docking was done with active site prediction tool (Ngan et al. 2012;Brenke et al. 2009). Pooled results exhibited matching of most of the residues, and thus, the binding sites used in docking approach were well defined (Supplementary material). As the exertion of biological activity of any drug is based on the binding between protein and ligand, present investigation results reveal that the nearly equal affinity of Erwinaze Ò towards both the ligands strongly supports Asparaginase and Glutaminase activities of the enzyme (Keating et al. 1993;Derst et al. 2000). Numerous adverse effects of L-Asparaginase usage have also been stated, distant from severe immunogenic reactions (Kwon et al. 2009).

MD simulations: protein in water
To check the structural behavior of modeled enzyme and docked complexes in the dynamic system MD simulations were performed and the trajectory analysis was done using GROMACS utilities. The total energy found from g_energy results for Erwinaze Ò modeled enzyme was found to   be -9,30,919 kJ/mol demonstrating that the model was energetically stable (Fig. 3). The structural convergence includes terms like root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg) and hydrogen bond (H-bond) analysis. In the g_rms based results with backbone atoms (Fig. 4), apo state enzyme is stable right from the launch of MD run with a drift at 10 ns followed by a stable confirmation till the finale of entire run. In case of docked complexes, though the complex 1 is unwavering till 20 ns it exhibited an increasing path later and extended the RMSD value to virtually 0.4 nm. RMSD of complex 1 with L-Asn till 10 ns specified the ligand binding in apo state and after 10 a significant conformational change with deviation was observed which might be the plausible cause to obtain L-Asn bound state. Its high RMSD value shows the instability of the complex which is undesirous. On the other hand, complex 2 has shown a perpetual RMSD value around 0.3 nm showing its high stability in the 300°K atmosphere supporting the     (Capizzi and Cheng 1981;Storti and Quaglino 1970) of L-Glutaminase nature of enzyme. This leaves the plasma pool L-Asn free for the uptake by cancerous cells on one hand because of unstable nature of complex 1, and causing the toxicity in normal cells due to its stable Glutaminase activity. Both the complexes are displaying their undesirable results against the actual desired L-Asparaginase role of enzyme in cancer therapy.
In the analysis of g_rmsf results, many oscillations in the residues with Ca atoms were observed in case of first complex with a fall and amendment in initial peak due to the binding of L-Asn in that region. This substrate binding also influenced the entire enzyme with several numbers of residual fluctuations with a near RMSF value of 0.3 nm (Fig. 5). The same plot also disclosed the presence of a second high peak in the vicinity of bounded residues establishing H-bonds with ligand leading to a very high RMSF value of around 0.6 nm in complex 2, leaving most of the other residues stable. The comparison of RMSF outcomes exhibited minor variations in ligand binding sites (Table 3) and their effect on complex formation.
Radius of gyration (Rg) was done to analyze the structural compactness of enzyme and was calculated using g_gyrate tool of GROMACS. Though the modeled enzyme in its free state has some structural compactness, bounded complexes were shown their fluctuating nature from Rg plots. Complex 1 has displayed the drifts during entire MD run in the array of 2.0 and 2.08 nm. In the plot of complex 2, it revealed its compactness with a little drift at 15 ns time point with the Rg values oscillating between 2.0 and 2.06 nm (Fig. 6).
To authenticate the hydrogen bonding pattern docking results, H-bond analysis was also executed along with the intact MD trajectories (Fig. 7). Surprisingly, the desired stable H-bonding pattern for complex 1 described missing linkage of hydrogen bonds between receptor and ligand. During 7.5 ns point the H-bonds were increased to 5, but after 4.03 ns h-bonds were missed and again formed from 6 ns. The same was again repeated at around 9 ns and 11 ns time points. Finally, after 12.6 ns there were no h-bonds, showing the disappearance of H-bond linkage between enzyme and substrate (Fig. 7a). The loss of intermolecular hydrogen bands may induce a spatial conformational change in the tertiary structure of the enzyme i.e., unstable structure of docked complex. This becomes a major shortcoming for the enzyme in therapeutic perspective, supporting possibility of its lower effect on cancer cells. On the other hand, it has shown the constant H-bond pattern with L-Gln throughout MD run with an average of 7 H-bonds throughout entire trajectory (Fig. 7b). H-bond pattern of both complexes did not draw a parallel with the docking result and this gives scope for the further investigation for efficient and firm drug in ALL therapy. Entire computational analysis describes the necessity of further intensive investigation on other sources of L-Asparaginase enzyme which is more effective on cancerous cells with fewer side effects.

Conclusion
L-Asparaginase is one of the most attractive enzymes in cancer research. Investigational studies on Erwinaze Ò sustained it as one of the therapeutic agent for ALL. In the present computational study, hypothetical model of Erwinaze Ò enzyme was developed and MD simulations for the same in apo state confirmed its stability. Ligand binding studies on Erwinaze Ò with L-Asn and L-Gln using molecular docking explained the binding mechanism of Fig. 7 Inter-hydrogen bonding for docked complexes a Erwinaze Ò ? L-Asn, b Erwinaze Ò ? L-Gln docked complexes. Key residues in binding were identified as Thr 15, Thr 26, Thr 27and Met 121 in L-Asn binding followed by Tyr 165, Thr 167 and Asn 180 in L-Gln binding based on the docking results. MD simulations of complex 1 for 30 ns run confirmed the instability and unlikeliness of enzyme towards L-Asn and the kinship with L-Gln in the dynamic system with stability. Quite a few conformational changes in enzyme were observed with Erwinaze Ò structure due to the L-Asn ligand binding compared to L-Gln binding. This in silico studies favor added extensive structural research on L-Asparaginase towards scheming of potential inhibitors that aim the L-Asparaginase for effective treatment of ALL.