In silico screening of the phytochemicals present in Clitoria ternatea L. as the inhibitors of snake venom phospholipase A2 (PLA2)

Abstract Millions of people suffer from snake bite envenomation, and its management is a challenge, even today. Medicinal plants have attracted the researcher’s attention for their outstanding advantages in treating many diseases, including snake venom poisoning. Clitoria ternatea L, is a plant popularly known for its various pharmacological effects especially, anti-snake venom property. However, the molecular mechanism behind this is poorly understood. It is reported that snake venom PLA2 is an extensively studied toxic factor. This study is meant to screen the compound’s capability to act as inhibitors of the Daboia russelli snake venom PLA2 through molecular docking and dynamics studies. Our results show that among the 27 compounds taken for the study, only Kaempferol showed good interaction profile with the conserved catalytic active site residues, His48 and Asp49. The pharmacophore features of the compound also demonstrate its exact fitting at the binding pocket. Further RMSD, RMSF, Rg, and hydrogen bond analysis confirmed the stable binding of Kaempferol with PLA2 through molecular dynamic simulations for 100 ns. In addition, the MM/PBSA binding free energy calculation of the complex was also affirming the docking results. The binding free energy (BFE) of Kaempferolis better than the reference compound. ADME and Lipinski’s rule of five reveals its drug like properties. Communicated by Ramaswamy H. Sarma


Introduction
Snake bite envenomation is a critical public health emergency in tropical and subtropical regions causing high rate of mortality and morbidity (White et al.,2003, Chippaux, 1998.Every year around 5.4 million snake bites occur which results 1.8 to 2.7 million cases of envenomation. The consequent death rate is three times more than other disabilities each year (Guti� errez et al., 2015). Africa, Asia and Latin America are the major places of incidence. Women, children and farmers are the main victims in rural areas of the developing countries as per World Health Organization (WHO) (Chippaux, 2017).
Snake venoms are a mixture of diverse peptides and proteins. Many of these components are responsible for their patho physiological effects in which secreted phospholipase A 2 (sPLA 2 ), is one of the dominant families of protein that was most extensively studied and is a reported toxic factor (Zelanis & Tashima, 2014, Slagboom et al., 2017, Mala et al., 2019. PLA 2 is highly present in Russell's viper snake venom, one of the deadliest snakes in India, causing high mortality. It hydrolyses phospholipids at its sn-2 position and releases fatty acids (Funk, 2001). Group IA and group II PLA 2 has been identified so far (Faure et al., 2010;Schaloske & Dennis, 2006). It is reported that a minimum of seven isoforms are present in viper venom PLA 2 (Warrell, 1989). Monomeric, dimeric, or trimeric PLA 2 exists among snake venoms depending on the nature of source and concentration (Dennis et al., 2011).It is composed of 121 amino acids in which 14 cysteine residues with catalytically active diad residues His48 and Asp49. (Burke & Dennis, 2009).
The only effective treatment for snake envenomation in this era is monovalent and polyvalent antivenoms, which has many limitations such as quality, storage, maintenance and geographical differences. These factors have made it imperative to think of safer alternatives (Pore et al., 2015, Kasturiratne et al., 2008, Gutierrez et al., 2010, Gutierrez et al., 2014. The information regarding the plants and the therapeutic detoxifying agents that we gained from folklores is a fundamental blessing and foundation for scientific community. Due to the apparent effectiveness, less side effects, availability and low cost, plant-based medicines are receiving increased acceptance (Pereanez et al., 2014, Sivaramakrishnan et al., 2016. C. ternatea, mainly present in South East Asia, especially in Indian subcontinent, belonging to Fabaceae family, is widely accepted for its different medicinal importance (Pati & Patil, 2011), particularly anti-snake envenomation property. Each part of the plant has medicinal importance and is used for various diseases. The plant exhibits different properties such as anti-inflammatory, anti-stress, memory enhancer, anti-depressant, anti-microbial, anti-diabetic, blood platelet aggregation-inhibiting, local anaesthetic, analgesic, and antipyretic (Mukherjee et al., 2007). With varied pharmacological activities C. ternatea is well known for its anti-snake venom activity (Dolui et al., 2004;Jaiswal, 2010).Traditional wound healers use this plant for the treatment of snake bite envenomation either by oral administration or external application of different parts of the plant in different forms. But no data for scientific validation, particularly molecular basis for its protective/curative action against snake envenomation seem to be available.
Wide varieties of computational studies exist related to drug development (Aanouz et al., 2020, Gupta et al., 2020 in which docking is a highly accepted and established in silico structure-based method to predict the binding mode of two or more compounds, its interaction and hence the mechanism. A possible poses of bioactive compounds were able to generate within a short period of time by cost effective manner. In the present study, as the understanding of intermolecular interaction between the active ligands with the macromolecule are most significant, docking studies are performed to elucidate the binding affinity of the biomolecules present in C. ternatea with the PLA 2 toxic venom protein.
Usually, the potency of the drug like molecule in a target was determined based on dock score.
No study has been reported so far on the molecular docking of the compounds present in the plant C. ternatea with Daboia russelli venom PLA 2 receptor. Therefore, this study aims to investigate the binding mode and interaction of the active principles present in C. ternatea with the target PLA 2 from Viper Russel's venom through docking and molecular dynamics. Results of the study predict Kaempferol as a potent PLA 2 docking molecule that may act as an antidote.

Preparation of protein structure
The crystal structure of Daboia russelli viper venom PLA 2 in complex with (2-carbamonyl methyl-5-propyl-octahydro-inol-7-YL) acetic acid (PDB id: 1OXL) (Alam et al., 2016) was retrieved from the Protein Data Bank (PDB) (http://www.rcsb. org). It exhibits 1.80 A � resolution. The protein was prepared, for docking by removing the water and hetero atoms using Biovia Discovery Studio 2018 and then defined the sphere in the active site. The active binding site residues His 48 and Asp49 (Burke & Dennis, 2009) of PLA 2 were used to define the active site of receptor by the option in Define and Edit Binding site, and the radius of the SBD_Site_Sphereset to 8.000. The values of X, Y and Z axis of the sphere defined are 15.077902, 21.606804 and 45.214753 respectively.

Ligand preparation
The structures of the 27 compounds from C. ternatea are given in Table S1. Varidiflorene, Quercetin-3-rutinoside and Quercetin-3-glycoside were drawn in Discovery Studio and others were directly retrieved from the pubchem database (http://pubchem.ncbi.nlm.nih.gov). These compounds were optimized, energy minimized by CHARMm force field and prepared for molecular docking using Biovia Discovery Studio v18.1.100.18065 client version. The ligand (2-carbamonyl methyl-5-propyl-octahydro-inol-7-YL) acetic acid in the crystal complex 1OXL, was removed, prepared and then redocked with the receptor using same protocol. This ligand is considered as reference ligand.

Molecular docking
CDocker in BIOVIA Discovery studio 2018, a systemic computational simulation method that utilizes a molecular dynamics (MD) simulated-annealing-based algorithm was used to generate more accurate complexes. The phytocompounds from the plant selected for the study and the reference compound were docked with PLA 2 receptor and the docked poses were ranked according to their CDocker energy.

Pharmacophore modelling
The pharmacophore feature of Kaempferol and the reference compound was identified through Automatic pharmacophore generation technique in Biovia Discovery Studio. The interaction of the key residues in the receptor with the ligand was taken as the query for finding the receptor based pharmacophore. This procedure was repeated for the crystal ligand.

Molecular dynamics simulations
Based on docking results, we selected the lowest energy valued and best posed PLA 2 -Kaempferol and the PLA 2-reference compound complexes for performing molecular dynamic analysis studies. The MD simulation was carried out by GROMACS 2021 with CHARMm 27 force field. The SWISSPARAM server was used to generate ligand topology. Before minimization, the overall system charge was neutralized by adding ions. The complex is solvated in the TI3P model, using the volume occupancy in triclinic box with periodic boundary conditions for neutralizing the system. The overall charge of the complex was solvated with appropriate cations and anions along with a salt concentration 0.1 mol/L. The system was subjected for energy minimization during MD, which was (nsteps¼ 50000) conducted using steepest descent approach (1000 ps) and was terminated at a threshold 10 kcal/mol. The Particle Mesh Ewald (PME) method was employed for the calculation of electrostatic interaction and Van der Waals interaction was predicted based on linear constraint solver (LINCS) (Amiri et al., 2007) algorithm. The standard temperature was kept at 300 K by applying modified Berendsen thermostat and Parrinello-Rahman coupling method was used to set pressure. Finally, a 100 ns molecular dynamics simulation was carried out for both the complexes. The MD trajectory analysis was also conducted for exploring  Homopterocarpin À 5.45 His48, Cys45, Ile19, Ala18, Gly6, Phe5, Leu2 9 10 Reference ligand À 31.07 Ile19, Gly30, Gly32, Asp49, His48, Tyr52 10 The active site residues Asp49 and His48 indicated in bold. No. of interaction denotes all types of contacts such as hydrogen bond, electrostatic attractive force, hydrophobic and Pi-Pi interactions.

Binding free energy of PLA 2 -ligand complex
The binding free energy of PLA 2 -Kaempferol and PLA 2-reference compound complex was calculated using the tool gmx_MMPBSA which implements the MM-PBSA (Molecular Mechanic -Poisson Boltzmann Surface Area) approach. It is extensively used in the study of biomolecular interactions. MM/PBSA calculates the binding free energy, which is associated with several energy terms through given equations (Kumari et al., 2014).
Where, G complex , G ligand , and G protein are the total free energy of the complex, ligand and receptor respectively, where each term in the right hand side of the equation was given by, In turn Where DH represents the enthalpy of binding and TDS is the conformational entropy after ligand binding (Wang et al., 2019). In the present study, 500 frames of span were used when performing MM/PBSA calculation with igb5 (GB-OBC2) model and a salt concentration ¼0.15 M.

Lipinski rule of five & ADME/T prediction
The ADME properties of the compounds were calculated using Swiss ADME program. Lipinski rule of five (Lipinski et al., 2001), pharmacokinetic properties, the solubilities of the drug, and drug likeness were also studied.

Results & discussion
In the current study, 27 compounds (Figure 1) of the plant C. ternatea were collected from the literature.
A total of 27 compounds were collected for the study. 2D structures of 24 of them were retrieved from Pubchem database and the remaining three sketched through Discovery Studio.

Molecular docking
Molecular docking can be used to model the interaction between a molecule and a protein at the atomic level and hence characterize its behaviour in the active site. Due to lack of adequate information on the molecular mechanism of antivenom property of C. ternatea, docking studies were carried out to examine whether inhibitors for PLA 2 were present. The efficiency of the docking of compounds with the target was inspected with respect to following parameters: interacting residues, number of hydrogen bonds, hydrophobic interactions, dock score, and the residues of the target involved in bonding. The most important step in virtual screening is the proper docking of compounds to the active site of target. Due to accuracy and reliability, we have chosen Cdocker was chosen for the current study. All the compounds docked individually with the PLA 2 through CDocker program. Among these, 9 compounds were docked, of which Kaempferol showed highest dock score to the target protein followed by Delphinidin, Isoparvifuran, and Myricetin-3-glucoside in that order (Table 1) whose dock score ranges from À 33.98 to 5.45.The further insights into the binding interactions of these compounds with PLA 2 were also analysed and recorded. Although, each compound showed interactions with the active site residues, only Kaempferol indicated the highest dock score and number of interactions, compared to the reference ligand (2-carbamonyl methyl-5-propyl-octahydroinol-7-YL) acetic acid) (-31.07).
The detailed intermolecular interactions at the binding site of Kaempferol with protein are shown in Figure 2. We could identify two types of interactions: five hydrogen bond and five hydrophobic interactions. Hydrogen bond interactions provide significant contribution to the strength of the small molecule-protein complex (Bauer et al., 2019, Panigrahi, 2008.Detailed analysis of Kaempferol showed to have hydrogen bonds with the catalytic residues His48, Asp49 and some neighbouring aminoacids (Balasubramanya et al., 2005).
A hydrogen bond interaction was observed in HIS48 with Oxygen atoms located at the second position of the compound in a distance of 0.28 nm. Other significant interaction in the docked complex is formed by the involvement oxygen in Asp49 at the active site through conservative hydrogen bond with H in a distance of 0.24 nm. The Hydrogen bond distances obtained in the present study comes under the allowed range (Daniel & Margaux, 2018). Additionally, other non-covalent interactions such as ionic and hydrophobic also play an important role in protein-ligand interactions. A noncovalent electrostatic force of attraction in His48 was also observed. Five different types of Pi-Pi interactions such as: two T-shaped, Pi-lone pair, Pi-sulfur and Pi-Alkyl hydrophobic interactions with residues such as Phe5, His48, Asp49, Cys45, and Cys29 respectively. For reference ligand, there was observed one pi-pi-Tshaped, one pi-alkyl and one alkyl hydrophobic interactions with Phe5 and Leu2 respectively and one electrostatic interaction in His48.
Molecular docking studies showed that the aminoacids His48, and Asp49 located at the binding site of the enzyme played a crucial role in enzymatic activity of PLA 2 .The molecular interaction results of the investigated molecule showed promising with the reference ligand indicating that Kaempferol possessed favourable orientation within the binding site that is similar to that of the reference ligand. So Kaempferol was taken for further studies.

Pharmacophore modelling
The pharmacophore model was also generated for Kaempferol and reference ligand based on their interaction with the active site of the receptor. A total of five different chemical features were determined, among them 3Hydrogen bond donors (HBD), 1 hydrogen bond acceptor (HBA), 1 aromatic ring, and 1 hydrophobic features were identified in Kaempferol ( Figure 3A), while the reference ligand showed 1 hydrogen bond acceptor (HBA), 1 ring aromaticity, 1 ionizable and 1 hydrophobic feature ( Figure 3B). In Kaempferol-PLA 2 complex, hydrogen bond donors are the predominant pharmacophore feature formed. This model represents the analysis of complementary chemical features of the active site and their spatial relationships. These features are responsible for the insertion of the compound which corresponds to the residues in the active site. These are selected based on the interactions generated from the structure information of the target. In this study it was focused to explore Kaempferol's pharmacophore features when it is inside the active pocket of the target. The phenolic compound, Kaempferol was investigated for its structural importance as a potential inhibitor against snake venom PLA 2 . The benzene ring and the alcoholic hydroxyl group impart acidic nature for the compound and hence it possesses the probability to inhibit snake venom PLA 2 . The formyl group of the compound can easily interact with macromolecules such as enzymes having toxic nature present in the venom, via chemical bond or chelation (Alam et al., 2016).

Molecular dynamic simulation analysis
The MD simulation was carried out to determine the structural and conformational stability of the protein (apo), protein-reference compound complex and protein-investigated compound complex for 100 ns of dynamic simulation. The structural changes in the protein and complexes and their dynamic behaviours were analysed through RMSD, RMSF, Rg and Hbond.
RMSD represents the extent of deviation of group of atoms with respect to the initial reference structure (Schreiner et al., 2012). Thus, low RMSD would be correlated to significant stability and ligands with high RMSD suggests inadequate ligand accommodation with the studied active site across 100 ns MD simulation .
The binding modes of the best docked molecules for both the reference ligand-PLA 2 and Kaempferol-PLA 2 complexes along with native protein were further opted for molecular dynamics simulation study in an aqueous system for a simulation time of 100 ns. In order to monitor its conformational and structural changes, the RMSD of the backbone atoms of apo protein and above mentioned complexes were carried out. The RMSD estimates extend of root mean square deviation of group of atoms with respect to its initial position (Schreiner et al., 2012). Figure 4 represents the three MD runs exhibited a fluctuation between 0.15 nm to 0.2 nm from �6ns till the end. For native protein, the fluctuation pattern is stable from �6ns to �62.5 ns. From there to 90 ns, it shows a minute instability and finally exhibits the uniformity. In the case of reference ligand-PLA 2 complex, almost a plateau-like appearance was observed beyond 6 ns. A similar style of pattern was exhibited by the Kaempferol-PLA 2 complex. From the overall RMSD analysis of 100 ns trajectory, the average RMSD for apo protein, the reference ligand-PLA 2  and Kaempferol-PLA 2 complexes were found to be 0.2 nm, 0.15 nm and 0.175 nm respectively. Though, near the end of simulation, from 90 ns to the end, it is observed that the RMSD of the investigated complex was same as that of apo and reference complex.
Further investigations of ligand stability within the binding site of the protein was proceeded through monitoring the ligand RMSD as it provides valuable information regarding the orientation of ligands with respect to their binding pocket. RMSD of reference ligand and Kaempferol were plotted with respect to time ( Figure 5).For Kaempferol, two intermittent inconsistency patterns of RMSD were observed during the first 7 ns and 11ns -16ns in correlation with reference ligand. These two portions were the possible areas where Kaempferol possessed different conformations during 100 ns simulation. Beyond the 17 ns, the compound exhibits stable conformation till the end when compared with the reference ligand with as RMSD of 0.225 nm.
To confirm these results and also to get more insights regarding the stability of the complex binding site, per residue root mean square fluctuation (RMSF) profile was calculated for each complex relative to apo protein. It estimates the average duration of each residue from its reference position within the original structure (Benson & Daggett, 2012).
Higher the value of RMSF indicates the structure has more flexible regions like loops and turns, whereas the lower value represents, the secondary structure, especially sheets and helices. The fluctuations for each residue were analysed for protein when it is at native form, complexed with investigated compound and with reference compound for 100 ns simulation period. Most of the protein show almost similar pattern of fluctuations except a few in Kaempferol complex. Figure 6 indicates that all residues in the complexes and apo protein have an RMSD value lower than 0.275 nm, indicating perfectly acceptable stable ligand-protein complex system.
Greater fluctuation pattern was observed in case of the complex of Kaempferolupto�0.21 nm. The fluctuating residues for the complex are Lys16 at � 0.15 nm, Gly30 at �0.21 nm, Ala40 at �0.125 nm, Ser70 at � 0.16 nm, Tyr75 at �0.165 nm, Gly80 at �0.18 nm, Glu85 at � 0.165 nm, and Leu130 at 0.12 nm. It is observed that the active residues such as His48 and Asp49 didn't show any fluctuation even after accommodating the ligand molecules in the active site. Beyond 85 th position the RMSF curve of the Kaempferol complex move as similar as apo protein and reference ligand complex.
Hydrogen bonds play a significant role in ligand binding with receptor. The bonding pattern was further assessed by observing the fluctuations of hydrogen bonds in both complexes during   100 ns simulation. The total number of hydrogen bonds present in the complex is shown in Figure 7A and B. Around four hydrogen bonds were observed in Kaempferol complex, while in reference complex five were computed. This bonding parameters detected indicates that the investigated compound binds to the target as effectively and tightly as reference drug.
Radius of gyration (Rg) was also monitored to investigate global stability of the complex. It is the mass-weighted RMSD for a group of atoms relative to common mass center in the whole complex structure (Lickic 0 et al., 2005).The Rg value of the Kaempferol complex in the present study was approximately 1.475 nm, whereas the reference compound and apo protein showed around 1.425 nm (Figure 8). This average value of Rg suggest that the complex is stably folded throughout the simulation as that of apo and reference complex, which in turn increases the rigidity of the protein structure and hence the overall stability.

Binding free energy calculations
The strength of the protein-ligand complex can be indicated via binding free energy. In the present study, gmx_MMPBSA was used to calculate binding free energy, in order to elucidate the affinity between the protein and ligand. It was done for the opted MD trajectories after its RMSD attained equilibrium. It is another term used to represent the stability of the docked complex in terms of all the intermolecular interactions between the ligand and the target. Previous studies showed that in the protein-ligand complex, various energy components such as van der Waals, electrostatic, polar and non-polar energy, render its own role to the binding free energy (Arun et al., 2021).Here, the major contributing interaction to the binding free energy of both Keampferol-PLA 2 and reference compound-PLA 2 complex were tabulated in Table 2.  Red color represents the pairs of atoms within the range of 0.35nm, and black color stands as Hydrogen bonds present within the limit. Here, 4 hydrogen bonds retained after molecular dynamics of 100ns, whereas in reference compound there are five hydrogen bonds. In both complexes electrostatic energy contribute a major role in total binding free energy. Non-polar solvation energy favourable in both the Kaempferol-PLA 2 and reference ligand-PLA 2 complexes (-3.14 ± 0.19 kcal/mol and À 1.37 ± 1.13 kcal/mol respectively).The present dynamic simulation was conducted for 100 ns, in which 10000 frames were incorporated. The standard error of mean of BFE in an average of 500 frames for both Kaempferol and reference compound. The net binding energy for the complex predicted by gmx_MMPBSA was established to be À 13.28 ± 2.24 kcal/mol for Kaempferol and À 8.12 ± 8.15 kcal/mol for reference compound was calculated. By comparing the BFE of these two complexes, it is clear that the investigated compound has a very good BFE rather than the reference compound. Despite these factors, it might prove to be essential in terms of potential drug, once preceded with advanced studies whose graphical representation was depicted in Figure 9.
Natural plants are rich sources of phenolic compounds, whose multiple biological functions stand to play important roles in pharmaceutical industry. Previous in silico studies reveal that the CrepisideE from Elephantopus scaber (E.scaber) was able to inhibit Viper venom PLA 2 by blocking the active binding site His 48, Asp49 and other neighbouring residues and hence a potential PLA 2 inhibitor (Mala et al., 2019). Kaempferol, a polyphenolic compound is rich in fruits and vegetables exhibiting anti-inflammatory effect by blocking various targets such as NF-OEB (Zhang et al., 2017), COX-1, and COX-2 (Garc� ıa- Mediavilla et al. 2007;Lee & Kim, 2010). In vivo assays have already established its anti-inflammatory role in atherosclerosis by inhibiting various targets such as TNF-a, IL-1b, cytokines, leukocytes etc (Kong et al., 2013). The anti-cancer effects of Kaempferol was also studied by triggering apoptosis via (Chen & Chen, 2013, Yao et al., 2016 inhibiting Caspase-9. Triggering of apoptosis by inhibiting the telomerase and PI3K/Akt pathways in cervical cancer by Kaempferol was also studied earlier (Kashafi et al., 2017). In addition to that, pharmacokinetic properties of the compound reveal its potency and efficacy as a candidate drug against the target. It satisfies the parameters of non-violation of Lipinski's rule. The logP value 1.9 implies that it has effective cell membrane permeability. As the successful completion of preliminary in silico studies, Kaempferol can be considered for next higher level in vitro experiments and prove it a potential drug lead against Daobia russelli viper venom PLA 2 inhibitor.

Conclusion
In the quest to find a novel drug candidate for snake venom PLA 2 , screening of compounds from plants which are commonly used by the traditional practitioners prove to be vital. In this context, identification of a potential inhibitor of PLA 2 of Daboia russelli viper venom has been achieved by docking studies from natural plant product Kaempferol, on the crystallized structureofPLA 2 (1OXL).Results derived from the docking studies have been validated by dynamic simulation, followed by binding energy calculation for Kaempferol which showed a good proximity with the reference crystal ligand, to confirm its stability. Preliminary in silico investigations showed that indeed molecule like Kaempferol has properties which can further be exploited and tested for its drug candidacy against Daboia russelli viper venom PLA 2 . In vitro studies are required to elucidate its molecular mechanisms to confirm its inhibitory activity.