This article was downloaded by: [University of Sydney] On: 07 September 2014, At: 20:59 Publisher: Taylor & Francis Informa Ltd Registered in England and Wales Registered Number: 1072954 Registered office: Mortimer House, 37-41 Mortimer Street, London W1T 3JH, UK

Journal of Biomolecular Structure and Dynamics Publication details, including instructions for authors and subscription information: http://www.tandfonline.com/loi/tbsd20

Identification of potent inhibitors against snake venom metalloproteinase (SVMP) using molecular docking and molecular dynamics studies a

b

a

Sathishkumar Chinnasamy , Selvakkumar Chinnasamy , Selvaraman Nagamani & Karthikeyan Muthusamy a

a

Department of Bioinformatics, Alagappa University, Karaikudi – 630004, Tamil Nadu, India

b

Department of Microbiology, Faculty of Medicine, Misurata University, Libya Accepted author version posted online: 05 Sep 2014.

To cite this article: Sathishkumar Chinnasamy, Selvakkumar Chinnasamy, Selvaraman Nagamani & Karthikeyan Muthusamy (2014): Identification of potent inhibitors against snake venom metalloproteinase (SVMP) using molecular docking and molecular dynamics studies, Journal of Biomolecular Structure and Dynamics, DOI: 10.1080/07391102.2014.963146 To link to this article: http://dx.doi.org/10.1080/07391102.2014.963146

Disclaimer: This is a version of an unedited manuscript that has been accepted for publication. As a service to authors and researchers we are providing this version of the accepted manuscript (AM). Copyediting, typesetting, and review of the resulting proof will be undertaken on this manuscript before final publication of the Version of Record (VoR). During production and pre-press, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal relate to this version also.

PLEASE SCROLL DOWN FOR ARTICLE Taylor & Francis makes every effort to ensure the accuracy of all the information (the “Content”) contained in the publications on our platform. However, Taylor & Francis, our agents, and our licensors make no representations or warranties whatsoever as to the accuracy, completeness, or suitability for any purpose of the Content. Any opinions and views expressed in this publication are the opinions and views of the authors, and are not the views of or endorsed by Taylor & Francis. The accuracy of the Content should not be relied upon and should be independently verified with primary sources of information. Taylor and Francis shall not be liable for any losses, actions, claims, proceedings, demands, costs, expenses, damages, and other liabilities whatsoever or howsoever caused arising directly or indirectly in connection with, in relation to or arising out of the use of the Content. This article may be used for research, teaching, and private study purposes. Any substantial or systematic reproduction, redistribution, reselling, loan, sub-licensing, systematic supply, or distribution in any form to anyone is expressly forbidden. Terms & Conditions of access and use can be found at http:// www.tandfonline.com/page/terms-and-conditions

Publisher: Taylor & Francis Journal: Journal of Biomolecular Structure and Dynamics DOI: http://dx.doi.org/10.1080/07391102.2014.963146

Identification of potent inhibitors against snake venom metalloproteinase

Downloaded by [University of Sydney] at 20:59 07 September 2014

(SVMP) using molecular docking and molecular dynamics studies Sathishkumar

Chinnasamya,

Selvakkumar

Chinnasamyb,

Selvaraman

Nagamania,

Karthikeyan Muthusamya* a

Department of Bioinformatics, Alagappa University, Karaikudi – 630004, Tamil Nadu, India

b

Department of Microbiology, Faculty of Medicine, Misurata University, Libya

Running title:

Identification of potent inhibitors against SVMP

*Corresponding author: Dr. Karthikeyan Muthusamy Assistant Professor Phone (Off): +91-4565-223344 Fax (Off): +91-4565-225202, Email: [email protected] 1

Abstract Snake venom metalloproteinase (SVMP) (Echis coloratus (Carpet viper) is a multifunctional enzyme that is involved in producing several of the symptoms that follow a snakebite, such as severe local hemorrhage, nervous system effects, tissue necrosis, etc. Because the three dimensional structure of SVMP is not known, models were constructed,

Downloaded by [University of Sydney] at 20:59 07 September 2014

and the best model was selected based on its stereo-chemical quality. The stability of the modeled protein was analyzed through molecular dynamics (MD) simulation studies. Structure-based virtual screening was performed, and fifteen potential molecules with the highest binding energies were selected. Further analysis was carried out with induced fit docking (IFD), Prime/MM–GBSA (ΔGBind calculations), quantum polarized ligand docking (QPLD) and density functional theory calculations (DFT). Further, the stability of the lead molecules in the SVMP active site was examined using MD simulation. The results showed that the selected lead molecules were highly stable in the active site of SVMP. Hence, these molecules could potentially be selective inhibitors of SVMP. These lead molecules can be experimentally validated, and their backbone structural scaffold could serve as building blocks in designing drug-like molecules for snake anti-venom. Key words: SVMP, Molecular dynamics simulation, induced-fit docking, Prime/MM– GBSA, QPLD and Density functional theory.

2

Introduction Snakebites are an important public health problem. Although it is difficult to be precise about the actual number of cases (Chippaux, 1998), every year approximately 35,000 to 50,000 people die from snakebites, and this is an especially serious concern in rural areas (Alam & Gomes, 2003; Bawaskar, 2004). Echiscoloratus and Naja n. Nigricollis are

Downloaded by [University of Sydney] at 20:59 07 September 2014

important envenoming snakes because their bites lead to clinical complications with systemic and local pathology. However, the venoms from the two species respond differently with respect to their pathophysiological functions (Chippaux et al., 1991). For example, an Echis coloratus (Carpet viper) bite leads to inflammation (such as swelling, blistering, and necrosis) and hemorrhages. Snake venoms are complex mixtures of proteins and peptides with a variety of biological activities. Snake envenomation is frequently tested with horse or sheep (Paul, 1993; Hasson et al., 2010). Snake venom metalloproteinase (SVMP) is an extracellular matrix (ECM) protein that belongs to the subfamily of zinc metalloproteinases (Bjarnason & Fox, 1994; Teixeira et al., 2005; Gutierrez et al., 2007; Pithayanukul et al., 2009; Lingott et al., 2012; Romero et al., 2012). Recent studies have suggested that SVMP is a mediator for edema, local tissue damage, inflammation and hemorrhage (Gutierrez & Rucavado, 2000; Lingott et al., 2012). SVMP is an enzymatic toxin because of the localized hemorrhage and damage to endothelial cells (Ownby et al., 1978). SVMP induces hemorrhaging, which subsequently leads to shock, tissue necrosis, hypotension, hypovolemia, inflammation and a reduced ability to regenerate muscle tissue (Gutierrez & Rucavado, 2000). SVMP is considered to be the key toxin involved in snake venom induced pathogenesis. Most of the time, SVMP induces hemorrhage by directly affecting the capillary blood vessels by clearing key bonds of the basement membrane components in a highly 3

selective fashion, and thus affecting the interaction between the basement membrane and the endothelium. It is evident that SVMP-targeted drug therapy would result in a beneficial outcome by reducing the mortality rate among patients (Yee, 2005; Gutierrez et al., 2007; Pithayanukul et al., 2009). In this study, a 3D structure of SVMP was constructed by homology modeling (SWISS-MODEL, Prime and Modeller) and was validated using PROCHECK, WhatIF analyses and the ProQ server. The stability of the modeled protein was

Downloaded by [University of Sydney] at 20:59 07 September 2014

analyzed using the GROMACS software (version 4.5) (Spoel et al., 2005; Hess et al., 2008; Ling et al., 2012; Karina et al., 2012). An average structure was generated from a 30 ns molecular dynamics (MD) simulation. Large collections of compounds were screened from different databases (Binding, Maybridge, Hitfinder and TOSLab Databases) against the selected target (SVMP). From these screening experiments, fifteen compounds were identified based on their potential interaction with SVMP. These fifteen compounds were then refined through various in silico techniques. Materials and methods Sequence alignment and Homology modeling The protein sequence of SVMP (Echis coloratus) (UNIPROT ID: E9JGE0) was retrieved from the UniProt database. To determine the sequences homologous to SVMP, the FASTA sequence of SVMP was submitted for BLAST (Basic Local Alignment Search Tool) analysis. The crystal structure of VAP2 from the Crotalus atrox venom (PDB ID: 2DW0) was selected as the template. The VAP2 has a 49% sequence identity with the SVMP query sequence. The homology modeling of SVMP was carried out by three different homology modeling programs (MODELLER9v7 (Sali & Blundell, 1993), SWISSMODEL (Nicolas & 4

Manuel, 1997) (http: //swissmodel. expasy.org/) and PRIME (Schrödinger, LLC, 2012)). The VAP2 from the Crotalus atrox venom (PDB ID: 2DW0) was selected as the template. The best geometry and energy profile was constructed using PRIME. Structure validation The modeled structure of SVMP was visualized using PyMOL (DeLano Scientific LLC, San Carlos, CA, USA, http://www.pymol.org). Superimposition of modeled structure

Downloaded by [University of Sydney] at 20:59 07 September 2014

was enclosed with crystal structure of VAP2 from the Crotalus atrox venom (PDB ID: 2DW0). The backbone of the modeled structure was calculated by analyzing phi (Φ) and psi (ψ) torsion angles using PROCHECK. (Sali & Blundell, 1993; Laskowski et al., 1993). Further, the modeled structure (SVMP) was evaluated by ProSA analysis. The energy criteria of the modeled structures were compared with large set of known protein structures. In addition, the quality of the modeled protein was checked by the WhatIF analyses (Padmavathi et al., 2008) and ProQ server. Molecular Dynamics Simulations The stability of the modeled protein and the protein ligand complexes were analyzed through MD simulation studies (Karina et al., 2012). The MD simulation was performed using GROMACS version 4.5 (Spoel et al., 2005; Hess et al., 2008; Ling et al., 2012; Karina et al., 2012) with the GROMOS96 53a6 force field. The NVT and NPT ensembles (constant number of atoms, volume, pressure, and temperature) were applied to the system for the MD simulation studies (Lei et al., 2013). Each system was placed in the center of a 10 Å×10 Å cubic box with periodic boundary conditions and solvated by SPC (simple point charge) water molecules. Na+ counter ions were added to satisfy electroneutrality, and periodic boundary conditions were applied to the system. A stable environment was maintained by 5

applying a 300 K Berendsen temperature, and the coupling constants were set for 2 ps (0.1 for temperature and1.0 for pressure). Long-range electrostatic interactions were calculated with the Particle-Mesh Ewald (PME) (Darden et al., 1993) algorithm with a cutoff value of 1.0 nm, and a cutoff of 1.4 nm was set for van der Waals interactions (Lei et al., 2013). The LINCS algorithm was applied for bond constraints (Hess et al., 1997). The energy of the system was minimized using the steepest descent algorithm for 1,000 steps (Payne et al.,

Downloaded by [University of Sydney] at 20:59 07 September 2014

1992; Kirubakaran et al., 2011). The solvent, ions and ligand were equilibrated for 100 ps, and the heavy atoms of the protein were controlled by a harmonic restraint with a force constant equal to 1000 kJ mol-1 nm-2 using the NPT and NVT ensembles (Ling et al., 2012). Totally, 30 ns of MD simulations were performed for all the complexes and the trajectories were stored every 1.2 ps. Active site predictions The active site of the SVMP was identified through SiteMap (version 2.3) (Anne Marie & Sid, 2008). Contour maps were generated to show the binding cavity of the target protein. Different physical descriptors such as size, degree of enclosure, degree of exposure, tightness, hydrophobicity, hydrophilicity and hydrogen-bonding possibilities were used to identify possible binding sites. Site points that are most likely to contribute to the proteinligand binding are linked together. Virtual screening Large collections of compound databases were screened against SVMP using the virtual screening workflow (Schrödinger, LLC, 2012) in Maestro. The virtual screening is performed using a collection of packages including LigPrep (ligand preparation), QikProp (ADME prediction and filter), HTVS (high throughput virtual screening) and other structure 6

tools. HTVS was carried out with the Binding database, the Maybridge database, Hit finder and the TOSLab collection of the compounds, which together contain a total of 460,572 compounds. Molecular docking The Glide docking program was performed against modeled SVMP with best inhibitory molecules. The selected inhibitory molecules were performed by Glide XP of

Downloaded by [University of Sydney] at 20:59 07 September 2014

Glide module. The shape and properties of the receptors were represented on a grid by several sets of fields that progressively provide more accurate scoring to the ligand poses. The grid was set to be 12 Å X 12 Å 12 X Å for the analysis of docking study. Geometric or hydrogen-bonding constraints were not introduced for substrate docking. Experiments were carried out using the default parameters. The standard precision (SP) level was used for the generation and scoring of each ligand. The top score conformations of the molecules were further re-docked using extra precision (XP) algorithm. Finally the top fifteen compounds were selected for further analysis. ADME Predicting Activity The QikProp is an absorption, distribution, metabolism, and excretion (ADME) prediction program. The QikProp program was used to obtain the ADME properties of the compounds. The ADME predicts both physically significant descriptors and pharmaceutically relevant properties of the molecules, either in individuals or in batches. All the molecules were neutralized before being used by QikProp program (QikProp, version 3.4, Schrödinger, LLC, 2012). The program was processed in normal mode, and predicted 44 properties for the molecules, consisting of principal descriptors and physiochemical properties, along with a detailed analysis of log p P (octanol/water), QP% and log HERG. It also evaluated the 7

acceptability of the compounds based on Lipinski’s rule of five which is essential for rational drug design. Induced fit docking (IFD) The selected compounds were subjected to the IFD procedure using the Schrödinger package (Schrödinger, LLC, 2012). The IFD study was performed using the following three steps: (a) the ligand was docked to a rigid receptor model, (b) van der Waals (vdW) scaling

Downloaded by [University of Sydney] at 20:59 07 September 2014

factor 0.5 was applied to the ligand’s non polar atoms and the protein, (c) the protein structure was energy minimized by keeping it close to the original crystal structure, while removing bad steric contacts (Friesner et al., 2004). Energy minimization was performed using the OPLS 2005 force-field with an implicit solvation model. Quantum mechanical polarized ligand docking (QPLD) In QPLD, selected ligands were docked with the Glide docking tool, and charges were calculated by Qsite. All the screened compounds were docked in the active site of SVMP by QPLD. The charges of each atom were calculated by Qsite, a regular standard precision (SP). Glide docking followed by extra precision (XP) refinement was performed by generating 15 poses per docked molecule. The best poses were submitted to the QM-ESP (Quantum mechanics- Electrostatic potential) charge calculation at the B3LYP/3-21G* level of theory within the protein environment defined by the OPLS-2005 force field. Finally, the poses were redocked for another Glide run performing the ESP (Electrostatic potential) atomic charge calculations and obtaining the XP scoring models (QSite, version 5.5, Schrödinger, LLC, 2012). Binding free energy calculations

8

Prime MM-GBSA was performed using the OPLS 2005 force field to predict the free energy of binding for a set of ligands and their receptor. The binding free energy, (ΔGbind) equation is calculated using the formula ΔGbind=ER:L-(ER+EL)+ ΔGsolv+ΔGSA (Tono et al., 2005), where ER:L is the energy complex, ER+EL is the sum of the energies of the ligand and the protein, ΔGSolv is the difference between the GBSA solvation energy (surface area energy) of the complex and ΔGSA is the sum of the corresponding energies for the ligand and

Downloaded by [University of Sydney] at 20:59 07 September 2014

the protein (Alam & Naik, 2009; Ramakrishnan et al.,2012). Density functional theory (DFT) calculations The best active compound conformations from the virtual screening process were subjected to the DFT calculations. The DFT calculations were used to study electronic molecular features such as molecular electrostatic map, electron density, and frontier molecular orbital density fields (i.e., HOMO, LUMO), which can explain molecular properties and biological activity. The DFT calculations were performed with Jaguar, version 7.8. DFT calculations were carried out based on the solvation state. The DFT calculations were analyzed through Becke’s three-parameter exchange potential and the Lee-Yang-Parr correlation functional (B3LYP) (Becke, 1993; Lee et al., 1998) using a 3-21G* basis set level (Binkley et al., 1980; Gordon et al., 1982; Pietro, 1982) and PBF solvation. In the present study, 3D-molecular electrostatic potentials (MESP) V(r) at a point r were derived using the following equation because of a molecular system with nuclear charges {ZA} located at {RA}, and with the electron density ρ(r):

ρ (r ' )d 3r ' ZA −∫ | r − r '| A=1 | r − R A | N

V (r ) = ∑

9

In this equation, N represents the total number of nuclei in the molecule and the two terms refer to the bare nuclear potential and the electronic contributions, respectively. We computed the Jaguar dipole moment, molecular electrostatic properties, lowest unoccupied molecular orbital (LUMO), including MESP, and the highest occupied molecular orbital (HOMO) energy. We calculated the electrostatic potentials using the van der Waals contact surface area of the molecule. The overall molecular size and the positive electrostatic potentials were

Downloaded by [University of Sydney] at 20:59 07 September 2014

indicated by color coded surface values. The most negative electrostatic potential was colored deepest red and the most positive electrostatic potential was colored deepest blue. The intermediate yellow, orange, green shades indicate intermediate ranges of reactivity. Results and Discussion Homology modeling The SVMP sequence was obtained from UniProt (Boeckmann et al., 2003) (UNIPROT ID: E9JGE0). Using the results of the BLAST analysis, the crystal structure of VAP2 from Crotalus atrox venom was selected as a suitable template (PDB ID: 2DW0, Chain-A) for SVMP. VAP2 has a 49% sequence identity with the SVMP query sequence (Figure 1). The SVMP structure was modeled using three different tools: Swiss model, Modeller 9v7 and Prime. The predicted models were checked for psi and phi torsion angles using the Ramachandran plot (The torsion angle values for the SWISS MODEL, MODELLER & PRIME models were 81.4, 85.7 and 86.5, respectively). The best model (PRIME) was selected based on the stereochemical properties. The modeled protein is depicted in Figure 2. The protein contains six β sheets and five α-helices, and the β sheets are joined with α-helices. Structure validation

10

The PROCHECK program was used to validate the modeled structure. The SVMP showed 86.5% of the residues to be in the most favored region, 10.7% in the allowed regions and 0.8% in the disallowed region. The overall quality factor Errat is 77.477. In addition to PROCHECK, we checked the quality of the modeled protein based on WhatIF analyses and the ProQ server. All the contacts between amino acids were analyzed using the WhatIF program. The Z score of -1.58 (Close to -2) indicated that the modeled structure was of good

Downloaded by [University of Sydney] at 20:59 07 September 2014

quality. Further, we analyzed our modeled protein in the ProQ (Protein Quality Predictor) server, which gave a LG score of 3.466 and a MaxSub score of 0.369. The template protein in ProQ gave a LG score of 4.557 and a MaxSub score of 0.202. Based on the PROCHECK, ERRAT and WhatIF results, we confirmed the quality of the predicted 3D structure. The interaction energy of each residue was analyzed using a ProSA II energy plot. The result of the ProSA analysis shows a maximum number of residues with negative interaction energy and very few residues with positive interaction energy. The overall interaction energy of the model was 9.05 kcal/mol, which is quite similar to that of the template 1E19 (-9.67 kcal/mol). The structure was subsequently checked with a VERIFY-3D graph. The quality of the VERIFY-3D score was above zero in the graph and corresponds to the satisfactory side chain environment. The back bone RMSD value was 0.4 Å. In terms of geometry and energy profiles, the final model was well validated. The overall homology model suggested that the model was good enough to be a preparatory point for the next phase of docking studies. Molecular dynamics The MD simulation studies were performed to analyze the stability of the protein using GROMACS version 4.5 (Hess et al., 2008). The root mean square deviation (RMSD) increased at the beginning of the simulation. The modeled SVMP protein had a low RMSD value (0.2 to 0.4 Ǻ) for the backbone, which indicated stability and a stable dynamic behavior 11

of the structure. The MD simulation study showed that the energy of the molecule was constant throughout the simulation time period. The root mean square fluctuation (RMSF) values were between 0 and 0.7 Ǻ. Very few atoms of the RMSF were at the C and N terminal loop regions, indicating the accuracy and stability of SVMP;. The simulation studies also indicated that the backbone RMSD increased between 15000 ps and 20000 ps, and it was constant for the rest of the duration of the simulation. Thus, the modeled protein was quite

Downloaded by [University of Sydney] at 20:59 07 September 2014

stable (Figure 3). Active site prediction The active site of the SVMP protein was predicted using the SiteMap program in Schrödinger, LLC (New York, NY, USA). Interestingly, we found three amino acids among the active site residues that were present in the template protein structure (Gly105, Glu139 and His148). The close structural similarity was observed by the superposition of modeled and template (PDB ID: 2DW0). Both structures were contains five β sheets and four α-helices, and the β sheets are joined with α-helices. A conserved methionine (Met162) is present in the downstream of the consensus sequence (HEXXHXXGXXHD) as observed in template structure (Met357). The three histidines (His138, His142 and His148) that function as ligands of the catalytic zinc atom (for template His333, His337 and His343) and a glutamate (Glu139) residue which also functions as the general base in our model protein (for template Glu334) (figure 1). Calcium ion was identified opposite to the active site and close to the crossover point of the N- and C-terminal segments of M-domain (Ca2+ -binding site). The Mdomain is followed by the D-domain, which can be divided into “shoulder” (Ds) and “arm” (Da) segments which was also reported in the template structure.

12

Virtual screening and molecular docking The active site of SVMP was used in the virtual screening work flow in Maestro, against the Binding, Maybridge, Hit finder and TOSLab Databases. The ligands were prepared at pH 7.0 ± 2.0 using the large penalties of high energy ionization, and the Epik state was removed from the tautomer states. The protein was maintained by scaling the van der Waals radius by 1.0 Ǻ; the partial atomic charge was less than 0.25 Ǻ with default

Downloaded by [University of Sydney] at 20:59 07 September 2014

constraint parameters. We used the Glide flexible docking algorithm using the OPLS-2005 force field in each grid. Approximately 81,592 compounds were screened from the Binding, May Bridge, TOS Lab Collection and Hit finder databases using Glide HTVS. Further docking analyses were carried out using Glide XP mode in Schrödinger. Finally, fifteen compounds were identified from the Glide XP results. The Glide XP results are shown in the Table 1. Among these fifteen molecules, the molecule of BD 17344 has high binding affinity with the Glide score (-12.05). The molecules of BD 16837 and BD 7951 were interacted with more hydrogen bond interactions. The interacted amino acid residues were (BD 16837) ZN1, ASP101, ASP102, HIS148, ALA164, LEU166 and (BD 7951) ZN1, ASP102, THR103, ALA134, HIS138, GLU139. The ZN1 interaction was observed in all the 15 ligand molecules. In Figure 4, we have detailed the fifteen docked protein ligand complexes with the hydrogen bond interactions and the superimposition of the molecules in the binding pocket. ADME Predicting Activity The pharmacokinetic properties of the best 15 compounds (Table 2) are predicted using QikProp (QikProp, version 3.4, Schrödinger, LLC, 2012) from Schrödinger, LLC. The predicted significant ADME properties such as the QikProp predicted the log IC50 value for blockage of K+ channels (QPlogHERG), predicted the octanol/water partition coefficient 13

(QPlogPo/w), skin permeability (QPlogKp), aqueous solubility (QPlogS), and the brain/blood partition coefficient (QPlogBB); the number of violations of Lipinski’s rule of five (Rule Of Five) are listed in Table 2. The number of stars indicates the deviations from 95% of the known drugs. The percentage of human oral absorption is based on the number of metabolites, the number of rotatable bonds, logP, solubility, and cell permeability. All these predicted ADME properties were acceptable range.

Downloaded by [University of Sydney] at 20:59 07 September 2014

Induced fit docking Molecular docking simulation helps us understand the reasonable binding modes of the protein and ligand. The IFD is a mixed protocol for molecular docking and molecular dynamics. The receptor was made flexible, which helps refine the active site. The whole protein was used as the input for grid generation. Fifteen compounds from the end result of Glide XP docking were subjected to IFD. All compounds were perfectly bound to the active site of the protein with the docking score ranging from -13.73 to -7.73. Then, these compounds were subjected to further computational analysis. Quantum mechanical polarized ligand docking All the protein ligand complexes were subjected to the QPLD analysis. This docking method can provide a more accurate treatment for the electrostatic interactions and helps improve the docking accuracy. It mainly helps in assigning the proper states to the receptor and ligand molecules. The docking accuracy improved in QPLD because it used quantum mechanical-based calculations for the ligand molecule and the new charges applied to the ligand. The quantum mechanical calculations were carried out using Qsite (Silva et al., 2010). The QPLD results were also well correlated with the Glide XP results. All these results

14

revealed that the selected fifteen compounds bind well with the active site of SVMP, even though the scores were calculated at the quantum level. Binding free energy calculations The Prime/MM–GBSA method based on the docking complex and used to calculate the binding-free energy (ΔGbind) of the ligands. PRIME MM/GBSA solvation energy (ΔGbind) was ranging from -29.77 to -79.02 kcal/mol and favorable binding G-score was -10.03 to -

Downloaded by [University of Sydney] at 20:59 07 September 2014

12.03 (Table 3). Among these fifteen molecules, increased ΔGbind scores were observed for the molecules of BD 904 and BD 16837 (-79.02kcal/mol and -74.20kcal/mol). Further, the binding free energy was calculated from the QPLD output. Here, we found the similar binding free energy which is correlating with the Glide XP results (Table 3). Density Functional Theory (DFT) Calculations The electrostatic features of SVMP inhibitors have been investigated using DFT calculations. All the screened compounds were optimized at the B3LYP/3-21G* level. The reactivity of the molecules was analyzed using various parameters such as HOMO, LUMO and MESP. The fragile nature of all the compounds was explained by smaller values of both HOMO and LUMO. All these lead molecules (BD17344, BD905, BD 904, BD 25279, BD16837, BD13364, BD5228, BD17338, BD7951, BD26458, BD20708, BD16010, BD4991, BD5001, and BD15993) were very reactive because electrons can be transferred rapidly. Moreover, the stability of all the lead molecules was proved by the small HOMO and LUMO gap (Table 4 and Figure 5). Molecular dynamics simulation of protein ligand complexes The selected fifteen compounds were subjected to a MD simulation study to analyze the stability of ligands in the active site of SVMP. The fifteen different docked complexes 15

(SVMP_BD17344, SVMP_BD905, SVMP_BD904, SVMP_BD25279, SVMP_BD16837, SVMP_BD13364, SVMP_BD5228, SVMP_BD17338, SVMP_BD7951, SVMP_BD26458, SVMP_BD20708,

SVMP_BD16010,

SVMP_BD4991,

SVMP_BD5001

and

SVMP_BD15993) were subjected to 30 ns of MD simulation. The conformational changes of the protein-ligand complexes were analyzed through the MD simulation. The RMSD of the back bone atoms were obtained over the 30 ns of MD simulation trajectories. Figure 6 shows

Downloaded by [University of Sydney] at 20:59 07 September 2014

the RMSD of the trajectory for the complex with respect to the initial structure. The RMSD of the protein was always maintained at 0.5 Å after 2.5 ns, except for BD 904, BD 5001 and BD 25279. The protein RMSD of BD 904 and BD 5001fluctuated at 0.6 Å from 2.5 ns to 22.5 ns, but it increased abruptly to 1.5 Å at 22.5 ns and finally stabilized at 0.8 Å after 25 ns. The RMSDs of the ligands in the active site of SVMP protein are shown in Figure 7. Although the RMSD of the BD4991 ligand fluctuated at 0.5 Å from 2.5 ns to 12.5 ns, it increased abruptly to 0.6 Å at 12.5 ns and finally stabilized at 0.45 Å after 25 ns. The number of hydrogen bonds was calculated for the entire 30 ns simulation time. All lead compounds have at least three hydrogen bond interactions (Figure 8) for the MD simulation performed for 30 ns. Conclusion Our main objective of this work is to model the SVMP protein and to identify potent lead compounds against SVMP. The 3D structure of the SVMP was modeled by different software and the best model was selected based on the stereochemical properties. The stability of the modeled protein was checked with MD simulations. The best lead compounds for SVMP were screened using different databases viz. Binding, Maybridge, Hitfinder and TOSLab Databases. Fifteen potent hits were identified through the Glide score and Glide energy. The ADME properties of these fifteen compounds are also in an acceptable range. 16

We performed a different set of docking procedures viz. IFD and QPLD to confirm the docking accuracy and the binding site of SVMP. All the lead molecules were also optimized through DFT calculations. The fragile nature of all the compounds was explained by smaller values of both the HOMO and the LUMO. The binding free energies were also well correlated with the docking result. Finally, MD simulation studies were performed for these fifteen lead molecules to check the protein-ligand complexes (SVMP_BD17344,

Downloaded by [University of Sydney] at 20:59 07 September 2014

SVMP_BD905, SVMP_BD904, SVMP_BD25279, SVMP_BD16837, SVMP_BD13364, SVMP_BD5228, SVMP_BD17338, SVMP_BD7951, SVMP_BD26458, SVMP_BD20708, SVMP_BD16010, SVMP_BD4991, SVMP_BD5001 and SVMP_BD15993) at 30 ns. All the compounds are quite stable in the active site of SVMP. From all these results, it is concluded that these lead compounds might show promising activity against SVMP. References Afroz Alam, M., & Naik, P. (2009). Molecular modelling evaluation of the cytotoxic activity of podophyllotoxin analogues. Journal of Computer-Aided Molecular Design. 23, 209–225. Alam, M. I., & Gomes, A. (2003). Snake venom neutralization by Indian medicinal plants (Vitex negundo and Emblica officinalis) root extracts. Journal of Journal of Ethnopharmacology, 86, 75–80. Anne Marie, J., & Sid, T. (2008). Driving forces for ligand migration in the leucine transporter. Chemical Biology & Drug Design, 72(4), 265–272. Baoping Ling, Min Sun, Siwei Bi, Zhihong Jing, & Yongjun Liu. (2012). Molecular dynamics simulations of the coenzyme induced conformational changes of Mycobacterium tuberculosis l-alanine dehydrogenase. Journal of Molecular Graphics and Modelling, 35, 1–10. Bawaskar, H. S. (2004). Snake venoms and antivenoms critical supply issues. Journal of the Association of Physicians of India, 52, 11–13. Becke, A. D. (1993). A new mixing of Hartree-Fock and local density-funtional theories. Journal of Chemical Physics, 98, 1372–1377. 17

Becke, A. D. (1993). Density-functional thermochemistry. III. The role of exact exchange. Journal of Chemical Physics, 98, 5648-5652. Binkley, J. S., Pople, J. A., & Hehre, W. J. (1980). Self–consistent molecular orbital methods. 21 Small split-valence basis sets for first-row elements. Journal of the American Chemical Society, 102, 939-947. Bjarnason, J. B., & Fox, J. W. (1994). Hemorrhagic metalloproteinases from snake venoms. Pharmacology & Therapeutics, 62, 325–372. Boeckmann, B., Bairoch, A., Apweiler, R., Blatter, M. C., Estreicher, A., Gasteiger, E., Martin, M. J., Michoud, K., Donovan. C., & Phan, I. (2003). The SWISS-PROT Downloaded by [University of Sydney] at 20:59 07 September 2014

protein knowledgebase and its supplement TrEMBL in (2003). Nucleic Acids Research, 31, 365–370. Catarina de Fatima Pereira Teixeira, Cristina Maria Fernandes, Juliana Pavan Zuliani, & Silvia Fernanda Zamuner. (2005). Inflammatory effects of snake venom metalloproteinases. Memorias do Instituto Oswaldo Cruz, 100, 181-184. Carlos Henrique Tomich de Paula da Silva, Vinicius Barreto da Silva, Resende, J., Rodrigues, P. F., Bononi, F. C., Benevenuto, C. G., & Taft, C. A. (2010). Computeraided drug design and ADMET predictions for identification and evaluation of novel potential farnesyltransferase inhibitors in cancer therapy. Journal of Molecular Graphics and Modelling, 28, 513–523. Chippaux, J. P., (1998). Snake-bites: appraisal of the global situation, Bulletin of the World Health Organization, 76(5), 515–524. Chippaux, J. P., Williams, V., & White, J. (1991). Snake venom variability: methods of study, results and interpretation. Toxicon, 29(11), 1279–1303. Darden, T., York, D., & Pedersen, L. (1993). Particle mesh Ewald: An Nlog (N) method for Ewald sums in large systems. Journal of Chemical Physics, 98, 10089–10092. Dongsheng Lei, Xing Zhang, Shengbo Jiang, Zhaodi Cai, Matthew Rames, Lei Zhang, Gang Ren, & Shengli Zhang. (2013). Structural features of cholesteryl ester transfer protein: A molecular dynamics simulation study. Proteins: Structure, Function, and Bioinformatics, 81, 415–425. Fabiaan Villalta Romero, Anna Gortat, Andrees Herrera, Rebeca Arguedas, Javier Quesada, Robson Lopes de Melo, Juan Calvete, Mavis Montero, Renato Murillo, Alexandra Rucavado, Josee Maria Gutieerrez, & Enrique Peerez-Payaa. (2012). Identification of New Snake Venom Metalloproteinase Inhibitors Using Compound Screening and 18

Rational Peptide Design. ACS Medicinal Chemistry Letters (ACS Publications), 3, 540−543. Friesner, R. A., Banks, J. L., Murphy, R. B., Halgren, T. A., Klicic, J. J., Mainz, D. T., Repasky, M. P., Knoll, E. H., Shelley, M., Perry, J. K., Shaw, D. E., Francis, P., & Shenkin, P. S. (2004). Glide: a new approach for rapid, accurate docking and scoring. 1. Method and assessment of docking accuracy. Journal of Medicinal Chemistry, 47, 1739–1749. Cheuk Yee, L. H. O. (2005). Purification of Metalloprotease Inhibitors That Neutralize Snake Venom Toxins in California Ground Squirrel Blood. Explorations: An undergraduate Downloaded by [University of Sydney] at 20:59 07 September 2014

Research Journal. 95–102. Glide, 5.0. Schrödinger, LLC; New York, NY, USA. (2012) Gordon, M. S., Binkley, J. S., Pople, J. S., Pietro, W. J., & Hehre, W. J. (1982). Selfconsistent molecular orbital methods. 22. Small splitvalence basis sets for secondrow elements. Journal of the American Chemical Society, 104, 2797-2803. Gutierrez, J. M., & Rucavado, A. (2000). Snake venom metalloproteinases: Their role in the pathogenesis of local tissue damage. Biochimie, 82, 841–850. Hasson, S. S., Jabri, A. A. A., Sallam, T. A., Balushi, M. S. A., & Mothana, R. A. A. (2010). Antisnake Venom Activity of Hibiscus aethiopicus L. against Echis ocellatus and Naja n. nigricollis. Journal of Toxicology, 2010, 1- 8. Hess, B., Bekker, H., Berendsen, H. J. C., & Fraaije, J. G. E. M. (1997). LINCS A linear constraint solver for molecular simulations. Journal of Computational Chemistry, 18(12), 1463–1472. Hess, B., Kutzner, C., van der Spoel, D., & Lindahl, E. (2008). GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation. Journal of Chemical Theory and Computation, 4,\ 435– 447. Jorgensen, W. L., & Duffy, E. M. (2000). Prediction of drug solubility from Monte Carlo simulations. Bioorganic & Medicinal Chemistry Letters, 10(11), 1155–1158. Jose Maria Gutierrez, Bruno Lomonte, Guillermo Leon, Alexandra Rucavado, Fernando Chaves & Yamileth Angulo. (2007). Trends in Snakebite Envenomation Therapy: Scientific, Technological and Public Health Considerations. Current Pharmaceutical Design, 13(28), 2935-2950. Karina, S., Matos, Elaine, F. F., da Cunha, Arlan da Silva Goncalves, Alan Wilter, Kamil Kuca, Tanos, C. C., Franca, & Teodorico, C. Ramalho. (2012). First principles 19

calculations of thermodynamics and kinetic parameters and molecular dynamics simulations of acetylcholinesterase reactivators: can mouse data provide new insights into humans. Journal of Biomolecular Structure and Dynamics, 30(5), 546–558. Kirubakaran, P., Kothapalli, R., Dhanachandra Singh, K. h., Nagamani, S., Arjunan, S., & Karthikeyan, M. (2011). In silico studies on marine actinomycetes as potential inhibitors for Glioblastoma multiforme. Bioinformation, 6(3), 100-106. Laskowski, R. A., MacArthur, M. W., Moss, D. S., & Thornton, J. M. (1993). PROCHECK: a program to check the stereochemical quality of protein structures. Journal of Applied Crystallography, 26(2), 283–291. Downloaded by [University of Sydney] at 20:59 07 September 2014

Lee, C., Yang, W., & paar, R. G. (1998). Development of the Colle-Salvetti correlationenergy formula in to a functional of the electron density. Journal of Physics Condensed Matter, 37, 785-789. Lipinski, C. A., Lombardo, F., Dominy, B. W., & Feeney, P. J. (1997). Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings. Advanced Drug Delivery Reviews, 23, 3–25. Nicolas, G., & Manuel, C. P. (1997). SWISS-MODEL and the Swiss Viewer an environment for comparative protein modeling. Electrophoresis, 18(15), 2714–2723. Ownby, C. L., Bjarnason, J., & Tu, A. T. (1978). Hemorrhagic toxins from rattlesnake (Crotalus atrox) venom. American Journal of Pathology, 93, 201-218. Padmavathi, G. V., Natraj Sekhar, P., Kavi Kishor, P. B., Vishal Kumar, A., & Saralakumari, D. (2008). Homology modeling and docking studies of human G-Protein coupled receptor involved in taste perception. International Journal of Integrative Biology, 2(1), 1-12 Paul, V. K., (1993). Animal and insect bites. In Medical Emergencies in Children, New Delhi, India, 2nd edition. Payne, M. C. et al., (1992). Rev Reviews of Modern Physics, 64: 1045 Pietro, W. J., Francl, M. M., Hehre. W. J., Defrees, D. J., Pople, J. A., & Binklley, J. S. (1982). Self-consistant molecular orbital methods. 24. Supplemented small splitvalence basis sets for second-row elements. Journal of the American Chemical Society, 104, 5039-5048. Pimolpan Pithayanukul, Jiraporn Leanpolchareanchai, & Patchreenart Saparpakorn. (2009). Molecular Docking Studies and Anti Snake Venom Metalloproteinase Activity of Thai Mango Seed Kernel Extract. Molecules, 14, 3198-3213. 20

Prime, 2.0. Schrödinger, LLC; New York, NY, USA, 2012. QikProp, version 3.4, Schrödinger, LLC, New York, NY, 2012. QSite, version 5.5, Schrödinger, LLC, 2012. QM-Polarized Ligand Docking Protocol; Glide, version 5.5; Jaguar, version 7.6; QSite, version 5.5; Schrödinger, LLC: New York, 2012. Ramakrishnan, C., Subramanian, V., Balamurugan, K. & Velmurugan, D. (2012). Molecular dynamics simulations of retinoblastoma protein. Journal of Biomolecular Structure and Dynamics, 2012, 1–16. Sali, A., & Blundell, T. L. (1993). Comparative protein modelling bysatisfaction of spatial Downloaded by [University of Sydney] at 20:59 07 September 2014

restraints. Journal of Molecular Biology, 234(3), 779–815. Schrödinger Suite 2012 Virtual Screening Workflow; Glide version. Schrödinger, LLC, New York, NY, 2012; LigPrep version. Schrödinger, LLC, New York, NY, 2012; QikProp version. Tono, C., Xu, G., Toki, T., Takahashi, Y., Sasaki, S., Terui, K., & Ito, E. (2005). JAK2 Val617Phe activating tyrosine kinase mutation in juvenile myelomonocytic leukemia, Leukemia, 19, 1843–1844. Torsten Lingott, Irmgard Merfort, & Thomas Steinbrecher. (2012). Free Energy Calculations on Snake Venom Metalloproteinase BaP1. Chemical Biology & Drug Design, 79, 990–1000. Van Der Spoel, D., Lindahl, E., Hess, B., Groenhof, G., Mark, A. E., & Berendsen, H. J. (2005). GROMACS: fast, flexible, and free. Journal of Computational Chemistry, 26, 1701–1718.

21

Downloaded by [University of Sydney] at 20:59 07 September 2014

Figguree 1.. Seequencce aaliggnm mentt beetweenn SV VM MP pprooteinn annd VAP22 froom Crrotaaluss attroxx veenom (P PDB B ID D: 2D DW00).

22 2

Downloaded by [University of Sydney] at 20:59 07 September 2014

Figure 2. Superimposition of modeled structure with the crystal structure of VAP2 from Crotalus atrox venom

23

Downloaded by [University of Sydney] at 20:59 07 September 2014

Figure 3. RMSD of the backbone atoms of the modeled protein (SVMP) over a period of 30ns.

24

Downloaded by [University of Sydney] at 20:59 07 September 2014

Figure 4. The fifteen docked protein ligand complexes with hydrogen bond interactions

25

Downloaded by [University of Sydney] at 20:59 07 September 2014

Figure 5. Three-dimensional-molecular electrostatic potentials (MESP) superimposed on to a surface of constant electron density (0.01 e ⁄ au3) showing the most positive potential region (deepest blue color) and the most negative potential region (deepest red color).

26

Downloaded by [University of Sydney] at 20:59 07 September 2014

Figure 6. The backbone RMSD of the protein the whole simulation line.

27

Downloaded by [University of Sydney] at 20:59 07 September 2014

Figure 7. RMSD of the selected fifteen ligands in the active site of SVMP.

28

Downloaded by [University of Sydney] at 20:59 07 September 2014

Figure 8. The total number of the Hydrogen bonds produced during the simulation time.

29

Table 1. List of compounds analyzed by docking simulations studies and their corresponding

Downloaded by [University of Sydney] at 20:59 07 September 2014

Glide scores, Glide energies and interacting residue S. No

Database ID

Glide score

Glide energy (kcal/mol)

1 2

BD 17344 BD 905

-12.05 -11.76

-56.43 -74.83

3

BD 904

-11.74

-74.47

4 5

BD 25279 BD 16837

-11.68 -11.14

-52.20 -53.47

6 7 8 9

BD 13364 BD 5228 BD 17338 BD 7951

-10.90 -11.04 -10.61 -11.09

-59.09 -46.17 -53.16 -49.77

10 11

BD 26458 BD 20708

-10.57 -10.42

-43.73 -57.46

12 13 14 15

BD 16010 BD 4991 BD 5001 BD 15993

-10.34 -10.03 -10.18 -10.07

-52.11 -36.66 -37.92 -51.04

30

Number H-bond interactions of SVMP of Hydrogen bonds 4 ZN1, ASP102, HIS138, GLU139 5 ZN1, ASP101, ASP102, THR103, HIS138 5 ZN1, ASP101, ASP102, THR103, HIS138 4 ZN1, GLY105, HIS138, ALA164 6 ZN1, ASP101, ASP102, HIS148, ALA164, LEU166 4 ZN1, ASP101, ASP102, LEU166 3 ZN1, GLY105, LEU166 3 ZN1, ASP102, HIS138 6 ZN1, ASP102, THR103, ALA134, HIS138, GLU139 3 ZN1, GLY105, HIS138 5 ZN1, GLY105, HIE125, ALA164, LEU166 4 ZN1, HIS138, CYX160, SER163 3 ZN1, GLU139, LEU166 3 ZN1, GLU139, LEU166 4 ZN1, HIS138, CYX160, SER163

Downloaded by [University of Sydney] at 20:59 07 September 2014

Table 2. ADME properties of the fifteen lead molecules predicted by the QikProp. S.L Compound No IDa

Mol MWb

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

471.55 569.70 530.66 434.32 562.45 542.70 331.41 439.47 314.36 337.39 525.61 372.38 270.32 290.31 322.32

BD 17344 BD 905 BD 904 BD 25279 BD 16837 BD 13364 BD 5228 BD 17338 BD 7951 BD 26458 BD 20708 BD 16010 BD 4991 BD 5001 BD 15993

Role of Fivec 1 2 2 1 1 2 0 0 1 0 2 0 0 0 0

Qplog HERGd

Qplog Po/we

Qplog Kpf

QplogSg

Qplog BBh

0.813 -4.500 -4.020 0.887 -2.657 0.136 -3.579 -2.505 -2.942 -3.824 -4.818 -4.467 -1.773 -2.189 -3.631

-1.514 2.414 2.140 5.001 0.392 0.310 -1.288 -0.740 -1.296 2.625 5.597 3.774 -1.823 -1.434 2.755

-5.929 -4.279 -4.210 -3.202 -6.396 -5.830 -7.608 -5.475 -6.038 -3.756 -2.952 -2.552 -6.424 -5.756 -3.174

-1.213 -3.724 -2.639 -4.140 -5.022 -3.918 -1.088 -2.466 -1.848 -4.109 -6.494 -4.611 -1.486 -1.364 -3.631

-1.487 -2.414 -2.173 -2.316 -2.226 -2.100 -0.841 -1.233 -1.660 -1.984 -2.363 -1.435 -1.189 -0.960 -1.555

a

Binding database compound ID

b

Molecular weight of the molecule (acceptable range: 130.0 - 725.0)

c

Number of permissible violations of Lipinski’s rule of five (acceptable range: maximum is

4) d

Predicted IC50 value for blockage of HERG K+ channels (concern below –5.0).

e

Predicted octanol/water partition co-efficient log p (acceptable range: –2.0 to 6.5).

f

Predicted skin permeability, log Kp (acceptable range: -8.0 to -1.0)

g

Predicted aqueous solubility; S in mol/L (acceptable range: –6.5 to 0.5).

h

Predicted brain/blood partition coefficient (acceptable range: –3.0 to 1.2).

31

Table 3. Docking score and binding energy fifteen lead molecules which bind to SVMP protein. Compound IDa

S.

Downloaded by [University of Sydney] at 20:59 07 September 2014

NO

XP Docking scoreb

Glide energyc(kc al/mol)

MM-GBSA ∆GBindd (kcal/mol)

IFD XP Docking scoree (kcal/mol)

QPLD XP Docking scoref (kcal/mol)

MM-GBSA ∆GBindg (kcal/mol)

1

BD 17344

-12.03

-56.43

-39.28

-11.07

-12.40

-49.36

2

BD 905

-11.74

-74.83

-72.48

-8.48

-12.12

-73.90

3

BD 904

-11.72

-74.47

-79.02

-11.56

-12.00

-79.13

4

BD 25279

-11.15

-52.20

-56.36

-11.47

-12.46

-59.80

5

BD 16837

-11.14

-53.47

-74.20

-12.38

-12.30

-74.79

6

BD 13364

-10.83

-59.09

-47.63

-11.17

-10.79

-67.36

7

BD 5228

-10.77

-46.17

-54.28

-10.37

-11.43

-55.78

8

BD 17338

-10.60

-53.16

-29.77

-11.40

-11.90

-28.89

9

BD 7951

-10.74

-49.77

-53.52

-11.51

-11.55

-53.72

10

BD 26458

-10.57

-43.73

-40.78

-13.73

-11.60

-33.91

11

BD 20708

-10.43

-57.46

-66.59

-10.65

-11.01

-64.60

12

BD 16010

-10.30

-52.11

-56.17

-9.71

-10.95

-54.80

13

BD 4991

-10.03

-36.66

-50.35

-7.73

-10.85

-45.34

14

BD 5001

-10.18

-37.92

-45.24

-11.51

-11.26

-43.58

15

BD 15993

-10.03

-51.04

-44.14

-9.13

-10.65

-48.42

a

Binding database compound ID

b

Glide XP docking score

c

Glide docking energy

d

Binding free energy calculated by Prime MM-GBSA

e

Induced Fit docking score

f

Quantum Polarized Ligand Docking score

g

Binding free energy calculated by Prime MM-GBSA

32

Table 4. Calculated by HOMO and LUMO gaps of selected molecules

Downloaded by [University of Sydney] at 20:59 07 September 2014

S. No

Compound IDa

HOMO (eV)b

LUMO (eV)c

Solvation energy

gap

(kcal/mol)e

1

BD 17344

-0.194

-0.007

-0.187

-0.09

2

BD 905

-0.206

-0.015

-0.191

-0.19

3

BD 904

-0.234

-0.011

-0.223

-0.18

4

BD 25279

-0.201

-0.037

-0.164

-0.36

5

BD 16837

-0.190

-0.068

-0.122

-0.19

6

BD 13364

-0.196

-0.003

-0.193

-0.22

7

BD 5228

-0.221

-0.011

-0.210

-0.22

8

BD 17338

-0.200

-0.036

-0.164

-0.10

9

BD 7951

-0.197

-0.012

-0.185

-0.16

10

BD 26458

-0.197

-0.060

-0.137

-0.15

11

BD 20708

-0.211

-0.021

-0.189

-0.13

12

BD 16010

-0.196

-0.050

-0.146

-0.13

13

BD 4991

-0.122

-0.011

-0.111

-0.11

14

BD 5001

-0.212

-0.008

-0.204

-0.13

15

BD 15993

-0.199

-0.046

-0.153

-0.14

a

Binding database compound ID

b

HOMOLUMOd

Highest Occupied Molecular Orbitals

c

Lowest Occupied Molecular Oribtals

33