Phytochemical analysis, molecular docking and molecular dynamics simulations of selected phytoconstituents from four herbs as anti-dotes for snake bites

The Irula community constitutes a small tribal community living in different parts of India and their main occupation is snake catching. They have rich knowledge about medicinal plants and their uses against various ailments. Snake bites are the common acute medical emergency faced by rural people. Even though Anti Venom Serum (AVS) is used as a remedy for snake bites in hospitals and in primary health centres, the rural people still depend on antidotes from medicinal plants. To counter the emergencies due to snake bites, the plants/parts of the plants are used alone or in combinations with other plants to make anti-dotes. The survey reports thirty traditional medicinal plants which are used by Irula against snake bites. These plants must contain some compounds responsible for snake venom inhibition. Venom neutralization or anti-inflammatory properties of the compounds identified from these plants have not been reported. Our present study reports the phytochemical analysis of the extracts from four plants which are commonly found and some of these compounds are identified as anti-inflammatory agents. Screened compounds are chosen for finding the binding affinity with PLA2 targets using the commercial Schrodinger software. Molecular Dynamics simulations were also carried out for the most favourable compounds and stability was checked upto 50ns. Correspondence to: D. Velmurugan, UGC-BSR Faculty, CAS in Crystallography and Biophysics, University of Madras, Chennai-25. E-mail: shirai2011@gmail.com


Introduction
The medicinal plants are the reservoirs of phytochemicals for curing human ailments and they play an important role in healing. They constitute an effective source of both traditional and modern medicines. Phytochemicals are naturally occurring in the different parts of plants (leaves, flowers, vegetables, root, stem, and bark) that have definite physiological action on the human body [1]. Herbal medicine has been shown to have genuine utility of about 80% of rural population which depends on it as primary health care (WHO, 2005). But, with the advent of modern technology and transformation of culture, this traditional practice has been decreasing gradually [2]. Irula is the tribal community in different parts of India, their occupation is snake catching and they know well about the medicinal plants and their uses against various diseases, particularly against snake bites [3]. Since snake bites lead to high mortality, the remedies are of great importance. To counter this, most of the tribal remedies are a combination of medicinal plants [4]. The survey conducted with the herbal healers and Irula tribes of Tamilnadu, resulted in thirty traditional medicinal plants, which are used by Irula against snake bites [5]. These plants must contain some compounds responsible for snake venom inhibition. Venom neutralization or anti-inflammatory properties of the compounds identified from these plants have not been reported. In this connection, phytochemical analysis of four herbal plants was carried out and the screened compounds were subjected to computational studies to find the binding mechanism with the suitable macromolecular PLA 2 targets. Many phytoconstituents that are present in the four plants are identified as anti-inflammatory agents.

Inflammation
A major component of snake venom was known to be Phospholipase A 2 (PLA 2 ) [6]. PLA 2 catalyses the hydrolysis of the ester linkage at the sn-2 position of phospholipids, leading to the production of free arachidonic acid and lysophospholipids [1]. This further digested by cyclooxygenase leads to the biosynthesis of proinflammatory compounds that are known as eicosinoids. The eicosinoids are implicated the triggering of inflammation. The PLA 2 family is classified into two forms; cytosolic isoforms and secretary isoforms, involved in signal transduction pathway and inflammation pathway [8]. The catalytic mechanism of PLA 2 is common throughout the family with the conserved three-dimensional structure, but only at the sequence level, differences can be observed [9,10]. PLA 2 exists as monomer, dimer and trimeric form depending on the source [11]. The three dimensional structure of PLA 2 includes three helices (H1, H2, H3), two short helices (SH4, SH5) and with β sheets (β1, β2) in antiparallel direction [12,13]. The catalytic triad comprises His 48, Asp49, Gly30, Cys45 residues and it possesses calcium coordination for its catalytic activity [14]. The calcium binding loop and Asp49 are most important for calcium coordination and catalytic mechanism. The venom PLA 2 does not have calcium ion coordination in the native and complex form. Calcium ion coordination can acquire only when it reacts with aggregated substrates [15,16]. The compounds, which are inhibiting PLA 2 , can be potent anti-inflammatory agents. Molecular docking was carried out with the screened compounds at the active site of PLA 2 . To confirm the stability of docked complexes and also for the analysis of radius of gyration, RMSD, RMSF, hydrogen bonds, etc.,

Sample collections
The fresh and healthy plants and leaves were collected from Chennai and its surroundings. The samples were washed under running tap water, shade dried for a few days (tuber was chopped into fine pieces) and crushed into fine power. 10g of each sample was soaked in to 100ml methanol and continuous stirring was carried out for overnight. The extract was filtered and the solvent was allowed to evaporate at the room temperature. Finally the extracts were collected and GC-MS analysis of these extracts was undertaken at Bureu Veritas Consumaer Products Services (I) Pvt. Ltd., Chennai, Tamilnadu. Mass spectrogram of GC-MS was analyzed using the database of about 1,50,000 compounds available in NIST library. The spectrum of the known compounds stored in the library was compared with the unknown compounds in the spectrum for identification.

Molecular docking
Molecular docking studies are used to determine the interaction of two molecules and to find the best orientation of ligand, which would form a complex with overall minimum energy. When a ligand molecule is docked at the active site of macromolecular targets, conformational changes may occur in some of the amino acid residues at the active site. In order to find the best fit at the active site, flexibility has to be provided not only for the ligand but also for the active site residues. However, this involves lot of computational time for the analysis hence to start with, the active site is assumed rigid and flexibility is provided only for ligand. In this way, a filtration can be made to few ligands when the input is involved with many ligands. After the selection of the small number of ligands using rigid docking, flexibility will now be provided to the active site region also. This is called Induced Fit Docking (IFD).
Molecular mechanics force fields are used to estimate the binding affinity between receptor and ligand that have been docked. The various components are contributing to the binding free energy. It can be written as

∆G bind = ∆G solvent +∆G conf +∆G int +∆G rot +∆G tor +∆G vib
The components consist of solvent effects, conformational changes in the protein and ligand, free energy due to protein-ligand interactions, internal rotations, association energy of ligand and receptor to form a single complex and free energy due to changes in vibrational modes [17]. . Molecular docking analysis was performed by using the commercial Schrödinger software [18].. The best conformation was chosen based on the docking score and glide energy. The sum of energy such as lipophilic, hydrogen bonding, metal interaction, rotatable bond counts and salvation contribute to the docking score. The glide energy is binding free energy calculated based on the OPLS-AA force field. The best compound can be selected based on the least glide energy/docking score/both.

Target preparation
The crystal structures of human PLA 2 at the resolution of 2.8Å (1DB5) [19] and Daboia russelii viper PLA 2 at the resolution of 2.7Å (2B17) [20]. Proteins were downloaded from Protein Data Bank (PDB). As explained in the introduction, the inhibitors of human PLA 2 make calcium coordination whereas with viper PLA 2 they do not. The proteins were prepared by removing water molecules, stabilizing the charges, fixing the missing residues and side chains by the prime module. The force field OPLS 2005 [21]. (Optimized Potentials for Liquid Simulations) was used for minimizing the energy of these proteins.

Ligand preparation
The compounds that are identified by GC-MS analysis were retrieved from Pubchem, NCBI Database [22]. These are reported as anti-inflammatory agents. The structures of the compounds were minimized after adding hydrogens, and correcting the bond orders by using steepest descent method with 1000 cycles and conjugate gradient method with 5000 cycles during energy minimization.

Molecular dynamics simulation
Molecular Dynamics simulation was performed using AMBER [23]. to determine the conformational changes and binding free energy. Crystal structure of human PLA 2 and venom PLA 2 were systematically prepared by using tleap module of Amber 12 suit. Molecular mechanics force field ff99sb amber was used to parameterize amino acids of protein; antechamber module was used to obtain the charges of ligand and parameterized by using GAF force field for system preparation and simulation [24,25]. The protein was solvated in a cubic TIP3P box of three -point charged water molecule at the 10Å marginal radius. The system was neutralized with adding adequate number of Na + and Cl + ions appropriately. Initially, the solvent molecules were relaxed while all the solute atoms and calcium ion were harmonically restrained to their original position with a force constant of 100kcal/mol.Å 2 . The water molecules were relaxed by 5000 steps of steepest descent algorithm and 2000 cycles of conjugate gradient algorithm. Finally, the whole system was subjected to energy minimization for 2500 iterations by conjugate gradient without restraints.
Before the simulation, the system was equilibrated in three phases with 2fs time integration step. Berendsen temperature coupling method [26] was used to regulate the temperature of 300K inside the box with time constant 2ps. Electrostatic interactions were computed using the Particle Mesh Ewald method [27]. In the second phase, by using isotropic position scaling, the constant pressure was given to the system. In these two stages, all non-solvent atoms were restrained. Finally, equilibration was extended at constant temperature and pressure without restraints. Consequently, simulations were performed in explicit solvent environment using NPT ensemble with 1 fs integration time step. SHAKE algorithm [28]. was used to constrain the bonds involving hydrogen atoms. The initial velocities were assigned from a Maxwell distribution at a given temperature. Finally, the system was computed to MD simulation for 50ns each on an 8 NVIDIA GPU build cluster. Potential energy of each system was validated by equilibration and simulation processes. Trajectories were used for the analysis of RMSD, RMSF, Rg and MMGBSA calculations.

GC-MS analysis
The peaks of the compounds with retention time in mins identified by GC-MS analysis of methanolic extracts of four plants are shown in Figures 1-4, respectively. The name of the compound, retention time and the percentage of their presence indicated as peak area are listed in Tables 1-4, respectively, for the four plants extracts. GC-MS analysis of the samples confirm the presence 14 compounds namely, Amyrin, Lupeol, α -Sitosterol/Campesterol/stigmasterol, Octadecanoic acid, n-Hexadecanoic acid, 9-Octadecenoic Acid, 9,12-Octadecadienoic acid, 9,12,15-Octadecatrienoic acid, 2-methyl-3-phenyl-1Hindole, Vitamin E, isophytol, 9-Nanodecene, Bicyclo Heptane and 22,23-Dihydrospinaterone from these four extracts and , the rest of the compounds are present comparatively in lesser proportions (not more than 3% only). Induced fit Docking module was used to predict the    Table 1).
The compounds which are covering more than 3% of peak areas are tabulated in Table 1. Rest of the compounds are not shown in the Table.1. Amyrin and lupeol had the same molecular formula and they are known to be similar compounds.
The compounds which are covering more than 3% of peak areas are tabulated in Table 2. Rest of the compounds are not shown in the Table.2. Stigmasterol and 22, 23-Dihydrospinasterone had the same molecular formula and they are known to be similar compounds ( Figure 2).
The compounds which are covering more than 3% of peak areas are tabulated in (Table 3). Rest of the compounds are not shown in the (Table 3) (Figure 3).
The compounds which are covering more than 3% of peak areas are tabulated in Table 4. Rest of the compounds are not shown in the Table.4.

Molecular docking
Docking was performed for all the 14 compounds identified from Tables 1-4, to confirm the binding of these ligands at the active site of human and Venom PLA 2 . These compounds are listed below in Table 5. The available Human PLA 2 structure is the PLA 2 complex with Indole-6 (PDB: 1DB5). This has interactions with His47, Asp48, Gly31, and Gly29 and also with calcium ion. In case of venom PLA 2 (PDB: 2B17), the structure is complex with diclofenac which interacts with His48 and Asp49. From the docking studies with venom and human PLA 2 , most of the compounds were found to show favorable docking score and glide energy ( Table 5). Out of 14 compounds, Amyrin,              Gly30 mediated by hydrogen bond interactions at a distance of 3.3 Ǻ, 2.7 Ǻ and 3.2Ǻ, respectively (Figure 5a) and other residues Tyr 22, Trp31, Cys 45, Tyr52 and Lys69 interact by hydrophobic interactions (Figure  5b). But in human PLA 2 , Vitamin E maintains only one interaction with Asp48 and it coordinates with calcium ion at a distance of 2.7 Ǻ, 2.8 Ǻ respectively (Figure.6a). Hydrophobic interactions of vitaminE with His47, Gly29, In case of amyrin with venom PLA 2 , it maintains hydrogen interaction with Asp 49 at a distance of 3.0 Ǻ (Figure 7a) and hydrophobic interactions with His48, Gly30, Tyr 22, Trp31, Cys45, Tyr52 and Lys69 are shown in Figure 7b. In human PLA 2 , amyrin interacts with Asp48, Cys44 and it coordinates with calcium ion at a distance of 2.8 Ǻ, 3.2 Ǻ and 3.1 Ǻ respectively (Figure 8a). Hydrophobic interactions of amyrin with His47, Gly29, Gly31, Cys44 and Tyr 21are maintained (Figure 8b).

MD simulation
Based on the docking studies, amyrin and octadecanoic acid bound complexes were chosen for MD study since the phytoconstituents amyrin and octadecanic acid are from the three plants, namely, Corollacarpus epigaeuss, Leucas aspera Spreng and Tinospora cordifolia. They showed minimum glide energy and have better binding affinity in terms of hydrogen bond and hydrophobic interactions with both venom and human PLA 2 ( Table 5). MD Simulation was performed to observe the molecular perturbation associated with the stability loss of the venom, human PLA 2 complexes and also the co-crystal complex.    The dynamic behavior of both complexes were calculated using RMSD, RMSF and R g. RMSD for all the Cα atoms were calculated from the initial structure to study the convergence of the protein system (Figure 14a). In a case of Venom PLA 2, cocrystal complex and amyrin bound complex showed similar trend of deviation till the end of simulation. At 15ns, the cocrystal complex remains much deviated than amyrin bound complex, resulting in final backbone RMSD of 2.0 to 2.3Ǻ during the simulation. Amyrin bound complex retained less deviation till the end of simulation. Octadecanoic acid bound complex exhibited higher fluctuation from 20ns to the end of simulation resulting in a backbone RMSD of 2.5 to 2.9 Ǻ. RMSD of human PLA2 is shown in Figure 14b. Cocrystal complex and octadecanoic acid bound complex deviate much at 42ns and 35ns, respectively, resulting in backbone RMSD of 2.2Ǻ. Amyrin bound complex retained less deviation till the end of simulation and it maintains the stability.

S.No
To observe the dynamic behavior of residues of different complexes of PLA 2 , atomic positional fluctuations of backbone residues of each PLA 2 were computed. In venom PLA 2, the fluctuation score exposed the occurrence of higher degree of flexibility in backbone residues present in the octadecanoic acid bound complex. Residues in amyrin bound complex depicted more fluctuations in the loop region (residues 35-40), (65-70), otherwise it maintains less fluctuations. Amino acid residues in cocrystal maintain less fluctuation compared with amyrin and octadecanoic acid bound complex residues and it maintains the stability (Figure 15a). Fluctuation score of amino acid residues in amyrin bound complex showed less degree of flexibility proved by more conformational restrictions in human PLA 2. In the co-crystal complex, residues deviate more than the residues in the amyrin and octadecanoic acid bound complex (Figure 15b). Active site residues of Cys45, His48 and Asp49 fluctuate very less for all complexes proving that this region is very stable in both venom and human PLA 2.
The radius of gyration (Rg) is the mass-weighted root mean square distance of group of atoms from their common centre of mass. Hence it provides an observation into global dimension of protein. In case of venom PLA 2 , major fluctuations are observed in co crystal complex and octadecanoic acid bound complex and it becomes less stable ( Figure  16a). In human PLA 2 , all the three complexes maintain the similar way of fluctuations till the end of simulation (Figure 16b).    The residue wise decomposition energy was calculated using the residues involved in hydrogen bond and hydrophobic interactions that can contribute to the delta free energy. In venom PLA 2 complex, Leu2, Cys29, Cys45 and His48 residues contribute the minimum range of delta free energy (Figure 17a). The major contribution of free energy comes from Leu2, Phe5, Ala18, Tyr21, and Gly22 for all human PLA 2 complexes shown in Figure 17b. In cocrystal complex, Asp48 showed maximum energy of about 13 kcal/Mol which is not favorable. In amyrin and octadecanoic acid bound complex, Cys44 and His47 maintain the least energy than the cocrystal complex.
The binding free energy is the sum of bonded and non-bonded interactions in the molecule. Here the binding energy is calculated for all venom and human PLA 2 complexes ( From the insilico analysis, it is proved that the best compounds amyrin and octadecanoic acid from the herbal plants behave like a cocrystal ligand at the active site of venom and human PLA 2 structures.

Conclusion
Medicinal plants are used for screening and discovering the phytoconstituents which are helpful for finding novel medicine and clues for modern medicine. In reviewing the above study, it can be concluded that the venom neutralization is possible with the phytoconstituents that are present in herbal plants or they can act as anti-inflammatory agents/antidotes for snake bites. Computational studies also clearly prove that the compounds have great medicinal and pharmacological action towards curing inflammation. Further experimental studies are needed to confirm the above finding.