Pharmacophore modeling, in silico screening, molecular docking and molecular dynamics approaches for potential alpha-delta bungarotoxin-4 inhibitors discovery.

BACKGROUND: The alpha-delta bungartoxin-4 (α-δ-Bgt-4) is a potent neurotoxin produced by highly venomous snake species, Bungarus caeruleus, mainly targeting neuronal acetylcholine receptors (nAchRs) and producing adverse biological malfunctions leading to respiratory paralysis and mortality.
OBJECTIVE: In this study, we predicted the three-dimensional structure of α-δ-Bgt-4 using homology modeling and investigated the conformational changes and the key residues responsible for nAchRs inhibiting activity.
MATERIALS AND METHODS: From the selected plants, which are traditionally used for snake bites, the active compounds are taken and performed molecular interaction studies and also used for modern techniques like pharmacophore modeling and mapping and absorption, distribution, metabolism, elimination and toxicity analysis which may increase the possibility of success.
RESULTS: Moreover, 100's of drug-like compounds were retrieved and analyzed through computational virtual screening and allowed for pharmacokinetic profiling, molecular docking and dynamics simulation.
CONCLUSION: Finally the top five drug-like compounds having competing level of inhibition toward α-δ-Bgt-4 toxin were suggested based on their interaction with α-δ-Bgt-4 toxin.

INTRODUCTION

Neuronal acetylcholine receptors (nAChRs) and their subunits including alpha-7, 8 and 9 are considered to be important receptor proteins in central nervous system which are responsible for several neuronal functions like signal transduction, modulatory effects on neurons in the nervous system, etc.[12] In the same manner muscular nicotinic acetylcholine receptor (muscle nAchRs) are also involved in nervous system function by stimulating Na+ and K+ ionic conductance. Both kinds of alpha 7 receptor of m-type and n-type AchR proteins are involved in many neuronal diseases like Parkinson's, Alzheimer's and schizophrenia.[345] The α7nAchRs are mainly present in brain, lymphocyte and spleen responsible for pre and postsynaptic excitation of nerve cells.[6] Apart from neuronal diseases and disorders, many potential neurotoxins act as antagonists and alarming the molecular functions of α7, α8 and α9 subunits of nAchRs.[7]

Alpha-delta bungarotoxin-4 (α-δ-Bgt-4) is a potent neurotoxin from elapidae family Indian snake called Bungarus caeruleus or Indian krait. From the recent years snake bites in India are increasing, the awareness and treatment strategies are comparatively slow and poor because of lack of antivenom, so the fatality rates in venomous snake bites are more in India.[8] The recent statistical study conducted across in India was reported the detailed snake bites and its average rate of deaths are 2, 50, 000/year.[910] The major snake bite deaths are caused by four highly venomous snake species, Daboia russelii, Echis cariatus, Naja naja and B. caeruleus are commonly called as “big four”.[11] Among these four snakes, Bungarus species causes death without showing local symptoms that are the main cause for death of the victim.[12] The venom of common krait contains most potent neurotoxins that have both presynaptic and postsynaptic neurotoxins and it stimulate muscular paralysis by affecting nerve ending situated near to the synaptic cleft of brain cells followed by respiratory paralysis, severe abdominal cramps, followed by death.[13] The krait bite is treated with antivenom treatment, and it shows several side effects like anaphylactic reactions that are viewed as risk to some of the victims.[14] The alternative way of treating the snake bite cases are using several plant based inhibitors compounds, which are used in ancient days and the people used folk medicine to treat the victims of poisonous snake, scorpions, etc., and it showed significant outcome against envenomation.[1516]

Many medicinally engrossed plants species were identified and used for several human ailments in earlier days. In each plant has 100's of bioactive compounds, and each one has their own biological and medicinal properties.[9] The two structures of selected bioactive compounds used to treat snake bites cases are given in Figure 1. The main purpose of this study is to endeavor present insights into the structural and functional role of α-δ-bgtx-4 and identification of potential α-δ-Bgt-4 inhibitors through in silico analysis, such as computational structure prediction, molecular dynamics (MD) simulation, pharmacophore mapping, pharmacokinetic and molecular docking analysis of α-δ-Bgt-4.

Figure 1

2D structures of selected bioactive phytochemicals used for snake bites. (a) Aristolochic acid I; (b) Edunol; (c) Wedelolactone; (d) Ellagic Acid; (e) 4-nerolidylcatechol; (f) Cabenegrin A-I; (g) Salireposide; (h) Curcumin; (i) Melanins; (j) Cabenegrin A-II

MATERIALS AND METHODS

Molecular modeling and molecular dynamics simulation

In order to identify the structural and functional information of α-δ-Bgt-4, the three-dimensional (3D) structure is considered to be an important component. The experimental structure of α-δ-Bgt-4 is unavailable in structural databases. Hence, α-δ-Bgt-4 structure was predicted using an automated homology modeling method using Modeller 9 v11.[17] The predicted 3D model was validated with structure analysis and verification server (SAVS) and Mol probity servers by analyzing amino acids distribution in Φ and Ψ of Ramachandran plot.[1819] Energy minimization was performed to the predicted 3D using Steepest Descent and Conjugate Gradient algorithms and it was allowed for MD simulation using Standard Dynamics Cascade program of Accelrys Discovery Studio (ADS) 2.0 for 1 nanosecond (1 ns) and the final stabilized model was obtained. From the trajectory analysis tool, potential energy and root mean square deviation (RMSD) were calculated. The final simulated model was used for further computational studies.

Identification and selection of antivenomic plants and their compounds

Information on antivenomic compounds of various medicinal plant species was collected from various literature resources. From the collection of plants and their compounds used for snake bites were segregated out and used for further computational studies. There were 25 bioactive compounds identified from the literature that has the antivenomic properties against venomous snake (including king cobra, cobra, krait, etc.) bites.[9] The pharmacologically active plant compounds and their structural analogs were retrieved from chemical databases using drug-likeliness filters. Pharmacokinetic properties were analyzed using “absorption, distribution, metabolism, elimination and toxicity (ADMET) descriptors analysis” module of ADS 2.0. In addition to pharmacokinetic analysis, the compounds that share the common chemical features called pharmacophore were drawn from the list of antivenomic plant compounds.

Pharmacophore model generation and computational virtual screening

The compounds with active antivenomic properties [Table 1] were used as base structures for the generation of pharmacophore using “common feature pharmacophore generation” program of ADS 2.0.[20] There are 10 hypotheses were generated with the same kind of parameters such as hydrogen bond acceptor, aromatic features (AA) using HipHop program. The “Ligand pharmacophore mapping” protocol of ADS was used for mapping the best compounds with good pharmacophore features. The best pharmacophore models were chosen based on fit values as one of the basic standard for selecting the ligand compounds that have a comparatively higher structural likeness from the resultant pharmacophore mapping of the database compound. Computational virtual screening protocols were followed based on the generated pharmacophore models against the drug-like compound database, Minimaybridge. All retrieved compounds were validated by checking the drug likeliness and pharmacokinetics (Absorption, Distribution, Metabolism, Excretion, and Toxicity) filters using drug likeliness and ADMET Descriptors module of ADS.[21]

Table 1

Common feature generation using selected bioactive phytochemicals

Binding site analysis and molecular docking and consensus scoring

Binding site of α-δ-Bgt-4 was predicted using three different binding site prediction tools, which help us to confirm the definite binding pocket and the constituted amino acid residues. Initially Q-site finder was used to predict the binding site, then it was cross-checked by CASTp server, and then a final validation was done in “Binding Site Prediction tool” of ADS 2.0.[222324] The predicted binding site with the volume of 26.500 Å3 was used for molecular docking analysis. The compounds good in both drug-likeliness features and ADMET properties were chosen and docked to the binding site of α-δ-Bgt-4 using AutoDock 4.0 and ADS-Ligand Fit programs.[2526] The resulting docked complexes were selected based on the binding energy, inhibition constant, VdW and electrostatic energies. The “ADS-Consensus Score protocol” was used to categorize the docked complex on the basis of the measured scores through multiple scoring functions such as Ligand Score 1 and 2, piecewise linear potential (-PLP1) and -PLP2, Jain, potentials of mean force, Ludi, Dock score, etc., The top scoring ligands and their conformations were taken for further MD simulation studies of docked complexes.[2627]

Molecular dynamics simulation on toxin-inhibitor complexes

Molecular dynamics simulation is one of an extensive protocol for calculating physical properties of biological molecules, especially proteins and in addition to this, MD simulation will help us to predict the essential properties of biological interest using “ADS-simulation module”.[2829] All selected best scored docked complexes of α-δ-Bgt-4-inhibitors were allowed for MD simulation for the period of 1 ns. Initially, each docked complex was stabilized by CHARMM force field, and hydrogen were added to both protein and ligand structures then protein structure conformation was corrected using “ADS-Clean Protein tool”. All sets of docked complexes were further analyzed using MD simulation under periodic boundary conditions in all directions to simulate the entire molecular system. The final production steps were carried out for 1 ns for each complex system under constant volume and temperature ensemble and the atomic coordinates of each complexes were updated at 1 picosecond (ps) level. The resultant simulated docked complex conformational changes and the corresponding RMSD were analyzed using superimposition tool available in “ADS-Analyze Trajectory and Superimpose Protein tool”.

RESULTS AND DISCUSSION

Three-dimensional structure prediction of alpha-delta-bungarotoxin-4

The sequence of α-δ-Bgt-4 was downloaded from the Uniprot sequence database (www.uniprot.org) with the accession number D2N116. The 3D structure of α-δ-bgtx-4 was predicted by homology modeling method, and the known experimental template protein structure was selected based on the identity (sequence identity is 35.6%). Most of the homology modeling programs predict the three dimensional protein structures based on the sequence identity more than 30% and yield reasonably good protein models, Modeller program is also one of them.[1730] The crystal structure of irditoxin from colubrid snake (Boiga irregularis) was taken as a template structure with the resolution of 1.50 Å from PDB database (PDB ID: 1H7Z B chain).[31] The secondary structure prediction and sequence alignment were done for alpha-delta-bungarotoxin-4 protein with the template protein using Jpred (www.compbio.dundee.ac.uk/www-jpred/) and ClustalW (www.genome.jp/tools/clustalw/) respectively. The selected template and target protein sequences were aligned perfectly, and most of the amino acids were found conserved. Based on the alignment, the initial 3D protein model was predicted using Modeller 9 v11. The predicted 3D structure validated using SAVS server and it shows 100% of amino acids are present in protein are in the allowed region of Ramachandran plot. In addition the structure was also validated using Mol probity server and it shows 90.5% amino acids in favored region and 100% residues are in the allowed region of Ramachandran plot and there is no amino acid found in disallowed region. These results clearly indicate the quality of predicted protein structure.

The protein model refinement done using energy minimization and final model obtained with lowest energy was taken for further optimization and MD simulation. MD simulation was performed for 1 ns and the molecular features analysis and the conformational changes were studied by calculating the RMSD between the conformations. The final stabilized model was chosen after careful analysis of potential energy, total energy and kinetic energy, etc., The predicted structure and detailed potential energy and RMSD plot along with the simulated structure of α-δ-Bgt-4 are given in Figure 2. There was a considerable decline of potential energy and the RMSD between the conformations were found between 1.6 and 1.7 Å and throughout the simulation, stable conformations were found. Hence, the simulated protein was used as a target protein for further in silico inhibitory analysis with selected phytochemicals used for venomous animal bites especially snake bites.

Figure 2

(a) The predicted and simulated structure of α-δ-bgtx-4. (b) Potential energy variation in molecular dynamics simulation. (c) RMSD analysis of simulated structure of α-δ-Bgt-4

Inhibitors selection, validation and common feature pharmacophore generation

As a result of rigorous literature search, small molecules with antivenomic properties were chosen for molecular interaction studies against of α-δ-Bgt-4. Several plant species were used for many highly venomous snake bites, especially for the “big four” were taken for this study. Then bioactive compounds were identified and retrieved from selected plant species and used for in silico analysis. The three dimensional structure of selected compounds were obtained from the PubChem database (www.ncbi.nlm.nih.gov/pc compound). The compounds are aristolochic acid from Aristolochia indica, 4-nerolidylcatechol from Pothomorphe peltata, wedelolactone from Eclipta alba, salireposide from Symplocos racemosa, cabenegrin A-1 from cabeca de negra, melanins from Lycium chinense, edunol from Harpalyce brasiliana were used for the treatment of snake bites and they were considered for this study to explore more number of inhibitors from pharmacophore based compound search.[93233] Pharmacophore models were generated with the selected bioactive molecules to explore the more number of desired drug-like molecules from the MiniMayBridge database of ADS 2.0. The best ten pharmacophore models were generated using HipHop program and their molecular overlays along with best-fit values are shown in the Figure 3 and Table 1.

Figure 3

(a) The generated pharmacophore model. Molecular overlay of best four Pharmacophore model of (b) CID 6325610. (c) CID 5281813. (d) CID14982. (e) CID44263865

The generated pharmacophore contains features like hydrogen bond acceptors, hydrogen bond donors, charge interactions, hydrophobic areas, aromatic rings with good fit values for top three hypotheses (3, 2.076, 2.036, 1.998) and rank (22.519, 20.30, 18.86, 18.24). This fitness value and rank indicates that how selected molecules shared common chemical features with each other. The top four models were used for further screening of drug-like compound search from MiniMaybridge database.[34] As a result of database screening, 916 compounds were retrieved and used for further pharmacokinetic studies.

Absorption, distribution, metabolism, elimination and toxicity descriptors analysis of selected small molecules

The ADMET properties of compounds play a crucial role in the drug discovery process as these are largely responsible for around 60% failures of drugs during various clinical phases.[35] In silico pharmacokinetic studies using ADMET descriptors analysis was carried to find the suitable inhibitors by analyzing their ADMET properties. From this analysis, among selected 916 ligands, 732 compounds showed poor results in either one or more ADMET parameters. As a result, only 184 compounds were showed better probability scores and the detailed result of this analysis plot is given in Figure 4. The compounds passed at all levels of ADMET descriptors analysis were subjected to molecular docking analysis to explore the mode of binding and other interactions with target protein, α-δ-Bgt-4 were analyzed using AutoDock 4.0 and ADS-Ligand Fit programs [Table 2].

Figure 4

Compounds passed at all levels of ADMET Descriptor analysis

Table 2

Ligand Fit scoring functions and details of interactions

Binding site prediction and molecular interaction studies on alpha-delta bungartoxin-4 with selected inhibitors

From the identified list of binding pockets of alpha delta bungarotoxin-4 using the Q-site finder and confirmed with ADS-binding site prediction tool, the largest one with key residues involved in inhibition was selected for docking studies. All compounds passed in ADMET descriptors analysis were used for further molecular interaction studies. The AutoDock 4.0 and ADS-Ligand Fit programs were used to execute the molecular docking based genetic algorithm (AutoDock) and shape-based complementarity and Monte Carlo methods (Ligand Fit).[26] The final docked poses were evaluated based on dock scores and hydrogen bonding interactions with the binding site residues. The resulting docked complexes were further validated their strength of binding and conformational changes in active site regions using MD simulation. From the result of docking, the selected ten compounds showed different kind of interactions with target protein and most of them involved in H-bonding with binding site amino acids and their binding energy, VdW energies are considerably better and the detailed interactions are depicted in Figure 5 and the energy terms and amino acids involved in interactions are given in Table 3. Among the selected ten compounds the first five compounds namely aristolochic acid-I, edunol, wedelolacetone, ellagic acid and 4-nerolidyleatechol were shown good interactions with α-δ-bgtx-4.

Figure 5

Molecular interaction of α-δ-Bgt-4 with selected inhbitors (a) aristolochic acid-I. (b) Edunol. (c) Wedelolacetone. (d) Ellagic acid. (e) 4-nerolidyleatechol

Table 3

Molecular docking and their energy calculation of α-δ-Bgt-4 and inhibitor complex

In a similar manner, drug-like database compounds were also allowed for molecular interactions and MD simulation studies. All selected 184 compounds with good pharmacokinetics profiles were taken for molecular docking studies. Among the selected compounds, the top eight molecules were chosen based on the various parameters including Ligand Fit scoring functions, binding free energy, H-bonding interactions with binding site amino acids, etc., From these studies, the drug-like database compounds such as BTB06769, BTB0384, BTB02340, BTB04932 and BTB5112 involved in active interaction with α-δ-Bgt-4 were selected based on scoring functions, energy calculations [Table 4] and all docked complexes were further analyzed with MD simulation and detailed structure of docked complex and their corresponding potential energy plots are illustrated in Figure 6.

Table 4

Ligand fit scoring functions of α-δ-Bgt-4 and selected inhibitor complex

Figure 6

Molecular interaction of α-δ-Bgt-4with selected inhbitors (a) BTB 06769. (b) BTB00384. (c) BTB02340. (d) BTB04932. (e) BTB05112. (f) BTB13358. (g) BTB13328. (h) BTB00722

Molecular dynamics simulation of docked complexes

In order to analyze the complex system of protein-ligand interaction and the stability of the interactions, MD simulation was performed over a period of 1000 ps with the temperature of 300 K for all selected top eight drug-like molecules with α-δ-Bgt-4. The initial and final potential energies were noted and analyzed for each bound complex. The complex structures showed stable conformation throughout the simulation process. Each docked complex was further analyzed using “Analyze Trajectory” tool and RMSD analysis. During the simulation studies, the compound BTB 00384 was shown better result and found that there wasn’t much significant variation in the conformation, and overall RMSD found ~ 1.9 Å. Hence, from the selected top eight compounds, BTB00384 is considered to be a good inhibitor for α-δ-Bgt-4 and then followed by BTB00722, BTB02340, BTB04932, and BTB06769.

CONCLUSION

In this work, we utilized various approaches, including computational molecular modeling of alpha delta bungarotoxin-4 protein, pharmacophore modeling and mapping of toxin inhibitors, molecular interaction and MD simulation studies which help to explore the suitable inhibitors from bioactive plant compounds for α-δ-Bgt-4. The predicted binding site amino acids of α-δ-Bgt-4 are Lys40, Glu58, Cys61 and Gly76 were shown to have a very prominent role in binding to the nAChR receptor protein. Among the selected bioactive phytochemicals aristolochic acid-I (CID13292117) from A. indica, cabenegrin A-I (CID13292115) and cabenegrin A-II (CID13292117) from Cabeca de negra and 4-nerolidylcatechol from P. peltata showed better interaction with α-δ-Bgt-4 toxin. Advanced techniques such as common feature pharmacophore models were developed based on compounds with good antivenomic effects such as aristolochic acid-I, cabenegrin A1 and A2 and 4-nerolidylcatechol. As a result of in silico studies, the following five compounds are selected namely BTB06769, BTB00384, BTB02340, BTB04932 and BTB007722, which are shown to have better interactions with the α-δ-Bgt-4 toxin. Subsequently, through MD simulation the stability of the docked complexes were examined. Hence, this identified compounds will help us to design novel and potential inhibitors, which are considered to be an alternative as well as a good antidote for B. caeruleus toxin, α-δ-Bgt-4.