1. INTRODUCTION
Snakebite envenomation is a major global health concern, causing over 100,000 deaths and 400,000 disabilities annually, predominantly in tropical and subtropical regions [1–4]. Rural and impoverished communities, particularly agricultural workers, are the most vulnerable [1,3,5]. Despite its severe impact, snakebite remains a neglected tropical disease, with limited data, restricted treatment access, and insufficient public health prioritization [1–5].
Although research on snake venom and the development of antivenom has increased [6], antivenom currently faces significant limitations, including high costs, limited access, and variability in effectiveness due to the diversity of venoms [7–11]. Additional challenges include adverse reactions and the need for intravenous administration, which hinders timely treatment in rural areas [11,12]. Comparative studies on venom-induced pathologies have also highlighted how different venoms can elicit varying systemic and local physiological effects, further complicating effective treatments [13]. These challenges highlight the urgent need for improved, broadly effective, and accessible therapeutic strategies for snakebite envenomation.
Phospholipase A2 (PLA2) enzymes are recognized as key contributors to snake venom toxicity, exerting a wide range of pathological effects in envenomation victims [14–16]. PLA2s are implicated in neurotoxicity, myotoxicity, hemostatic disturbances, and inflammation, and their toxic effects arise from both their enzymatic activity and specific protein–protein interactions [14–16]. Structural studies have revealed that variations in amino acid composition (AAC) and the presence of distinct toxicity domains contribute to the diverse toxicological profiles observed in different snake species [15,17,18]. The central role of PLA2s in venom-induced pathology underscores their importance as targets for therapeutic intervention in snakebite envenomation.
Medicinal plants have been increasingly explored as alternative sources of snakebite therapeutics, particularly in regions where access to conventional antivenoms is limited or unavailable [19–21]. Ethnobotanical surveys and systematic reviews have documented the use of hundreds of plant species across diverse geographic areas, including Central America, Sri Lanka, India, Tanzania, and Ghana, to treat snakebite envenomation [20–24]. Scientific investigations have demonstrated that certain plant-derived compounds can inhibit key toxic enzymes and neutralize the deleterious effects of snake venom, supporting their potential as complementary or alternative therapies to antivenom [25–27]. However, while pre-clinical studies and traditional knowledge highlight the promise of medicinal plants, further pharmacological validation and clinical research are necessary to ensure their efficacy and safety in snakebite management.
The use of herbal remedies as alternatives or adjuncts to conventional anti-snake venom therapy has several limitations. One major challenge is the lack of standardization and quality control of herbal preparations, which can result in variable efficacy and safety profiles [28]. Additionally, many herbal antivenoms have not undergone rigorous in vitro and in vivo testing, and there is a scarcity of well-designed clinical trials to validate their therapeutic potential and determine the optimal dosing regimens [29,30]. The degradation of active compounds during extraction and formulation, as well as poor solubility and bioavailability, further hinders the effectiveness of herbal products [28]. These limitations highlight the need for more comprehensive pharmacological studies and the development of standardized protocols to ensure the reliability and safety of herbal antivenom therapies.
Cyanthillium cinereum has emerged as a promising candidate for anti-snake venom therapy. Recent in vitro studies have demonstrated that aqueous root extracts of C. cinereum possess significant venom-neutralizing activity against the venoms of Naja naja and Daboia russelii, two medically important snake species [31]. These findings suggest that C. cinereum could serve as a formidable therapeutic agent for snakebite envenomation, warranting further investigation and development.
Computational methods, especially AI-powered techniques, have dramatically enhanced the systematic and efficient exploration of plant-derived compounds for anti-snake venom applications. Tools such as drug-target interaction (DTI) prediction, molecular docking, and molecular dynamics (MD) simulations enable the rational identification of potential inhibitors by modeling their interactions with venom proteins at the molecular level [15,32]. For instance, molecular docking studies have demonstrated that phytochemicals such as 5-methylcoumarin-4-β-glucoside and lupeol from Vernonia glaberrima can bind with high affinity to the active site of PLA2, a key venom enzyme, thereby supporting their potential as venom inhibitors [32]. Beyond traditional simulation, artificial intelligence now plays a pivotal role, including genetic algorithms and artificial neural networks (AIS models), which further facilitate the classification and selection of medicinal plants with venom-neutralizing properties by analyzing relevant features and amino acid sequences [15]. Additionally, machine learning approaches such as Random Forest (RF), XGBoost (XGB), and stacked autoencoder (SAE-DNN) have been successfully employed to predict peptide-protein interactions involving venom-derived peptides, showing their potential for accelerating antivenom discovery [33]. Taken together, the integration of AI-enhanced DTI prediction, molecular docking, MD simulations, and predictive models for pharmacokinetics/toxicity provides a robust and streamlined framework to accelerate the development of effective, plant-based anti-snake venom agents.
This study provides an overview of the state of the art in plant-based approaches for snakebite management, highlighting the potential of computational strategies to accelerate the discovery of PLA2 inhibitors. Building on this context, we employed machine learning-based DTI prediction, molecular docking, and MD simulations to evaluate C. cinereum compounds, aiming to identify natural inhibitors that may complement existing antivenom therapies.
2. MATERIALS AND METHODS
2.1. DTI prediction
2.1.1. Data acquisition
DTI data were acquired from the Snakebite Envenoming Medicines Database [34], which contains products and investigational drug candidates for various snake venoms. The whole database was downloaded, and the data were filtered by choosing those that were categorized as drugs in the product column. Approximately 102 drugs, which consist of 44 unique compounds and 42 unique protein targets, were collected.
The full dataset was generated by first shuffling the pairs of unknown interactions of compounds and proteins to obtain the negative interactions. Five negative interactions were produced for each compound by randomly selecting five non-interacting proteins. The total number of compound-protein pairs in the training data was 322, with 102 positive interactions and 220 negative interactions.
For the test data, we pooled all the distinct compounds of sawi langit (C. cinereum), which contained 303 compounds, and protein targets from the dataset, yielding 12.978 compound-protein pairings.
2.1.2 Feature extraction
Feature extraction was performed using two different descriptors for compounds and proteins. For compounds, we used PubChem fingerprint [35] and Morgan fingerprint [36], whereas for proteins, AAC and ESM 2 embeddings [37] were used.
The PubChem fingerprint is a substructure-based fingerprint that represents the chemical substructure of a compound, with each substructure assigned to a specific location in the array. PubChem fingerprints are frequently used in virtual screening due to their interpretability and fast computation, as well as covering a broad and diverse range of biologically relevant [38,39]. For each substructure that occurs in the compound, the position of the corresponding substructure in the fingerprint vector is 1; otherwise, the position of the substructure is 0. PubChem fingerprint extraction yielded 881 features that explained the distinct substructure of a compound. The Morgan fingerprint is a topology-based fingerprint that is formed by analyzing the number of molecular fragments that emerge from a specified path or radius of a molecule. Each path or radius is then encrypted with a hash. The bit value in a topology-based fingerprint array is 1 if there is a molecular fragment at a particular path length or radius; otherwise, it is 0 [39]. Morgan fingerprint extraction resulted in a 2,048 array for compound features. The illustration of the fingerprint feature can be seen in Figure 1 [40]. Morgan fingerprints are regarded as state-of-the-art chemical descriptors in DTI tasks because they are fast to compute, capture critical substructure information, and consistently serve as practical features in machine learning models [41,42]. These two fingerprints provide complementary structural information: PubChem emphasizes biologically relevant substructures, while Morgan captures molecular topology and variations in connectivity. We evaluated them separately to determine which representation better supports compound–target interaction (CTI) prediction in this context.
 | Figure 1. (a) Illustration of the substructure-based fingerprint. (b) Illustration of the topology-based fingerprint. [Click here to view] |
AAC represents the protein sequence by calculating the normalized frequencies of the 20 common amino acids. The AAC extraction result was an array with a length of 20 for every native amino acid [43].
Evolutionary Scale Modeling 2 (ESM2) is a transformer-based language model designed for protein sequence analysis. ESM2 encodes crucial details about protein features, including evolutionary, structural, and functional traits, by creating embeddings directly from raw sequences. ESM2 embedding resulted in 1,280 features that represented protein sequences. ESM2 embeddings are considered state-of-the-art protein features in DTI and protein–ligand modeling, frequently outperforming or replacing traditional sequence encodings due to their superior ability to represent protein structure [37].
In this study, we evaluated AAC and ESM2 embedding separately to explore different levels of protein representation. AAC provides a straightforward composition-based feature set, which is computationally efficient and effective for capturing general sequence information. In contrast, ESM2 utilizes deep learning to extract rich contextual and structural insights from raw sequences. By assessing each representation independently, we aim to determine which type of protein descriptor yields better performance in predicting DTIs.
2.1.3 Prediction model
To construct robust DTI prediction models, we employed three complementary machine learning algorithms: RF, XGB, and a deep neural network pretrained with a SAE-DNN. The rationale behind combining these models was to leverage their distinct strengths: tree-based ensemble learning for stability and interpretability, and deep learning for capturing complex, nonlinear relationships.
We used the SAE-DNN model to enhance feature representation in high-dimensional input spaces. Unlike conventional DNNs initialized with random weights, SAE pretraining provides optimized initial weights and biases, which enhance both convergence and generalization [44]. Since the prediction task was a binary classification, the output layer of the DNN used a sigmoid activation function. At the same time, we need to tune hyperparameters to achieve optimal performance.
For ensemble-based learning, we chose RF and XGB due to their proven robustness in handling high-dimensional bioinformatics datasets [45,46]. RF employs bagging to reduce variance and improve stability, whereas XGB applies boosting to minimize errors and sequentially enhance predictive accuracy.
We used two types of molecular descriptors for all models: compound fingerprints (PubChem and Morgan) and protein descriptors (AAC and ESM2 embeddings). This multimodal representation ensured that both chemical and protein features were captured, thereby improving the reliability of predictions.
We optimized hyperparameters to maximize predictive performance. For SAE-DNN, we employed Bayesian optimization due to its high-dimensional hyperparameter space, where exhaustive search would be computationally infeasible [47]. In contrast, we optimized the hyperparameters of RF and XGB using GridSearchCV [48], as their hyperparameter sets were smaller and more interpretable, allowing for an efficient exhaustive search. Table 1 presents the final hyperparameters for each model.
Table 1. Hyperparameter search space.
| SAE-DNN hyperparameters | |
|---|---|
| Hyperparameter | Values |
| HL0 node | 100–2000 |
| HLi node | 0.5 × (HL i-1)–0.75 × (HLi-1) |
| Hidden layer | 1–6 |
| Learning rate | 0.01–0.1 |
| Dropout rate | 0.2–0.7 |
| L2 kernel regularizer | 0,0.001,0.01 |
| RF hyperparameters | |
| n_estimators | 100–500 |
| Max_samples | 0.1–1, None |
| Max_depth | 3–6, none |
| XGB hyperparameters | |
| eta | 0.01–0.1 |
| colsample_bytree | 0.1–1 |
| n_estimators | 100–500 |
2.1.4 Model evaluation and prediction
The training dataset consisted of 102 positive interactions and 220 negative interactions, resulting in a moderately imbalanced class distribution (~1:2 ratio). To address this imbalance without introducing synthetic noise, no resampling method such as SMOTE is used, as such a process may generate artificial patterns and risk biasing the model, especially for high-dimensional datasets [49,50]. To balance the number of positive and negative interactions, the model was trained using stratified K-fold cross-validation with several folds, K = 5. This approach helps to reduce class imbalance and provides a more reliable estimate for model performance, especially in data sets where the number of negative samples is greater than that of positive samples. Accuracy, recall, precision, F-score, receiver operating characteristic curve (ROC), and precision-recall curve (PRC) were the evaluation metrics. The accuracy quantifies the predictive power of the test findings. The precision metric quantifies the proportion of accurate predictions compared to a positive class. Recall gauges how well the model predicts a favorable outcome. The performance of the minority class was measured using the F-score. The trade-off between specificity (1-FPR) and sensitivity (TPR) was displayed in the ROC curve [51]. To assess the effectiveness of the model in predicting the class under unbalanced conditions, the PRC curve compares precision and recall [52].
For each prediction model, the best model was chosen based on its metric score, particularly in terms of ROC and PRC. ROC and PRC are two commonly used metrics for model comparison in machine learning due to their ability to capture performance across classification thresholds, especially under class imbalance [53,54]. The prediction of C. cinereum compounds was then conducted using the best model of each algorithm, and the prediction results were compared to obtain similar prediction results.
2.2 Molecular docking
2.2.1. Preparation of 3D structures of target proteins (PLA2)
Three-dimensional (3D) structures of PLA2 from three different organisms (Crotalus durissus terrificus, Daboia russelii, and Bothrops asper) were retrieved from the Protein Data Bank (PDB) [55]. The structural quality of each model was assessed using PROCHECK [56] to evaluate the stereochemical properties. Refinement was conducted using GalaxyRefine2 [57], and the structural quality was compared before and after refinement by analyzing residues in the most favored regions in the Ramachandran plot. The model with the best quality was selected for molecular docking simulation. UCSF Chimera [58] was used to remove water molecules and final protein structures.
2.2.2 Preparation of 3D structures of ligands
The 3D structures of the compounds from C. cinereum, along with antivenom drugs Varespladib and Batimastat, were obtained from the PubChem Compound Database [59] in .sdf format. These structures were converted into .pdbqt format using OpenBabel [60] to ensure compatibility with the molecular docking simulations. During this preparation, OpenBabel was used to assign the dominant protonation states of the ligands at a physiological pH of 7.4. The prepared ligands were then used as input data for docking analysis conducted using AutoDock Vina [61,62].
2.2.3 Molecular docking simulation
Molecular docking simulations were performed for the top ten pairs of compound-protein from DTI prediction using AutoDock Vina within UCSF Chimera. The refined 3D structures of both the target proteins and ligands were imported into the software. A docking grid with predefined coordinates (X, Y, Z) and dimensions of 45 × 45 × 45 Å was established to cover the whole protein. The docking procedure was executed using the appropriate commands within the software. After docking, binding affinities (ΔG values) and root mean square deviation (RMSD) values were calculated. Binding poses, molecular interactions, and docking scores were visualized using Ligplot+ [63] and PyMOL [64]. These tools facilitated a detailed examination of bond types and key residues involved in protein–ligand interactions, providing insights into the stability and effectiveness of the docking results.
2.3 MD simulation
MD simulations were conducted to assess the stability of the top-performing ligands with PLA2 isoforms identified from molecular docking studies using GROningen MAchine for Chemical Simulations (GROMACS) version 2024.2 [65,66]. The simulations were performed on an NVIDIA DGX A100 system equipped with a 2g.20 GB slice of an NVIDIA A100 GPU, 16 CPU cores, and 64 GB of system memory. The protein topology was generated using the Chemistry at Harvard Macromolecular Mechanics 36 (CHARMM36) force field [67,68]. In contrast, ligand topology files were obtained from the SwissParam server with an MMFF-based approach, which provides parameters compatible with the CHARMM36 all-atom force field [69].
The protein–ligand complexes were placed inside a cubic simulation box with a 10 Å buffer between the molecular structure and box boundaries and solvated using the CHARMM-modified TIP3P water model. To maintain system neutrality, Na+ or Cl− ions were added as counter ions, and energy minimization was conducted using the steepest descent algorithm [70].
The system underwent equilibration in two phases: an NVT ensemble for 100 ps, followed by an NPT ensemble for 200 ps, using the V-rescale thermostat. Position restraints were applied to both protein and ligand to preserve their initial conformation. Isotropic pressure coupling was applied using the Parrinello-Rahman barostat at 1 bar, with restraints still applied. A 100 ns Production MD simulation was then executed under constant pressure (~1 atm) and temperature (~300 K). All bonds were constrained using LINCS, and a 2 fs integration timestep was applied. Trajectory snapshots were saved every 10 ps, and the “trjconv” module in GROMACS was employed to eliminate coordinate artifacts and correct periodic boundary conditions.
Structural parameters such as RMSD, root-mean-square fluctuation (RMSF), radius of gyration (Rg), and hydrogen bond formation (HBond) were extracted from the simulation trajectory. Visualization and graphical analyses were performed using the Matplotlib library in Python. The binding free energy of the top hits and reference compounds was calculated by the MMGBSA method using the gmx_MMPBSA tool, integrated with GROMACS and employing the AMBER force field [71].
3. RESULTS AND DISCUSSION
3.1. DTI prediction
3.1.1. SAE-DNN results
SAE-DNN was trained using the optimal parameter: HL0 node = 256, learning_rate = 0.0001, Hidden Layer = 3, dropout rate = 0.5, HLi node = 0.66. The average of SAE-DNN results using five cross-validation can be seen in Figure 2a.
 | Figure 2. (a) SAE-DNN average for five CV, (b) ROC, and (c) PRC curve of Morgan-ESM2 SAE-DNN model. [Click here to view] |
Among these results, the model trained using Morgan-ESM2 produced the best metric score compared to others, with a good score of accuracy (0.7925) and recall (0.72215), and reasonably satisfactory results in precision (0.64147) and F-measure (0.67837). A good score in accuracy indicates that the model is quite reliable in its overall class predictions, but owing to data imbalance, recall, precision, and F-measure metrics are more suitable for interpreting the model’s performance. A good recall score and satisfactory results in precision and F-measure indicate that the model has a good performance in its positive prediction. However, it needs to improve the precision and F-measure. The ROC and PRC curves of the model are shown in Figure 2b and c.
The Morgan-ESM2 model produced a good AUC-ROC of 0.83, suggesting that the model achieves a high true positive rate while maintaining a low false positive rate, making it effective for overall classification. However, due to data imbalance, the ROC curve might overestimate the model’s effectiveness. The AUC-PR of 0.64 implied that while the model is reasonably effective at identifying true positives, there is still room for improvement in reducing false positives. Fold 5, which produced the highest results among other folds, seems more predictive, resulting in 0.78 AUC-PR.
3.1.2 RF results
RF model trained using optimal parameters: n_estimators = 100, max_depth = None, max_samples = None. The average of the RF results using five cross-validations can be seen in Figure 3a.
 | Figure 3. (a) RF average score for five CV, (b) ROC, and (c) PRC curve of Morgan-AAC trained RF. [Click here to view] |
RF produced the best results when trained using Morgan-AAC features, with a score of 0.82642 in accuracy, 0.62366 in recall, 0.76364 in precision, and 0.67353 in F-measure. In contrast to the SAE-DNN model, the RF model produced a higher precision than recall, indicating that it is better at avoiding false positives (non-interacting pairs) than at finding all positive interactions. Both models produce a similar F-measure. The ROC and PRC curves of the model can be seen in Figure 3b and c.
The Morgan-AAC model of RF produced better results in terms of mean AUC of ROC (0.87) and PR curve (0.78) than the Morgan-ESM2 SAE-DNN model. The Morgan-AAC model is capable of distinguishing positive and negative classes, marked by its high AUROC, and has a good performance in identifying true interacting pairs while balancing false positives, as shown by its high AUPR.
3.1.3 XGB results
XGB was trained using the optimal parameters: colsample_bytree = 0.2, eta = 0.3, n_estimators = 100. The average of XGB results using five cross-validations can be seen in Figure 4a.
 | Figure 4. (a) XGB average results for five CV, (b) ROC, and (c) PRC curve of Morgan-AAC trained XGB. [Click here to view] |
The XGB model produced quite similar results to the RF model, which excels when trained using Morgan-AAC features. The model produced a good score in accuracy (0.82036), recall (0.66946), precision (0.72327), and F-measure (0.69082). While the model produced a higher precision than recall, the same as RF, the XGB model is more balanced, with a recall score higher and a precision score lower than RF. The F-measure of this model also had the highest score compared to the other two models, even though the difference is insignificant. The ROC and PRC curves of the model can be seen in Figure 4b and c.
RF and XGB models produced similar results in terms of accuracy, recall, precision, and F-measure, but XGB scored lower in AUROC (0.85) and AUPR (0.69). These results surpass those produced by SAE-DNN, with the 5th fold being the best among the other folds.
3.1.4 Test data prediction results
The prediction of this test data data aims to predict the potential active compounds contained in C. cinereum as anti-cancer agents. The output of this prediction is the pairs of active compounds from C. cinereum with cancer-associated target proteins. Predictions were made using the best performing model from each machine learning method, namely Morgan-ESM2 SAE-DNN, Morgan-AAC RF, and Morgan-AAC XGB.
Out of 12.978 compound-protein pairings in the test data, 48 pairs emerged from the prediction results of the models, with 35 unique compounds and six unique target proteins. Compound 2-Amino-3-(4-hydroxy-3-methoxyphenyl)propanoic acid exhibited the most interactions among other compounds with five venom proteins, followed by 4-Coumarylalcohol, 4-Hydroxybenzoic acid, and trans-3-Indoleacrylic acid, which interacted with three venom proteins.
Among the prediction results, the pair of 2-Amino-3-(4-hydroxy-3-methoxyphenyl) propanoic acid (CID: 1670) and PLA2 of Crotalus durissus terrificus (South American rattlesnake) had the highest average probability of 0.848. In the top 10 pairs of protein and compounds with the highest average probability, there are three distinct snake species, namely Crotalus durissus terrificus (South American rattlesnake), Daboia russelii (Russell’s viper), and Bothrops asper (Fer-de-lance). The details of these pairs can be seen in Table 2.
Table 2. Top 10 pairs of compound-protein with the highest average probability.
| No | CID | Compound name | Uniprot accession | Protein name | Organism | Probability |
|---|---|---|---|---|---|---|
| 1 | 1670 | 2-Amino-3-(4-hydroxy-3-methoxyphenyl)propanoic acid | P08878 | PLA2 | Crotalus durissus terrificus (South American rattlesnake) | 0.847919 |
| 2 | 5375048 | trans-3-Indoleacrylic acid | A8CG86 | PLA2 | Daboia russelii (Russell’s viper) | 0.843646 |
| 3 | 135602137 | 2-hydroxyphenazine-1-carboxylic acid | Q6JK69 | PLA2 | Bothrops asper (Fer-de-lance) | 0.837415 |
| 4 | 10821 | Phenylacetylene | Q6JK69 | PLA2 | Bothrops asper (Fer-de-lance) | 0.81839 |
| 5 | 5376187 | Sorbicillin | Q6JK69 | PLA2 | Bothrops asper (Fer-de-lance) | 0.802828 |
| 6 | 5280535 | 4-Coumarylalcohol | Q6JK69 | PLA2 | Bothrops asper (Fer-de-lance) | 0.788209 |
| 7 | 5375048 | trans-3-Indoleacrylic acid | Q6JK69 | PLA2 | Bothrops asper (Fer-de-lance) | 0.774928 |
| 8 | 442793 | [6]-Gingerol | Q6JK69 | PLA2 | Bothrops asper (Fer-de-lance) | 0.767288 |
| 9 | 135 | 4-Hydroxybenzoic acid | A8CG86 | PLA2s | Daboia russelii (Russell’s viper) | 0.760714 |
| 10 | 125429 | Hydantocidin | Q6JK69 | PLA2s | Bothrops asper (Fer-de-lance) | 0.759898 |
In addition, there are six distinct proteins from the prediction results of the three models, identified by their Uniprot Accession, namely Q6JK69, A8CG86, P08878, P00605, and P15445. Q6JK69 and Q8AXY1. PLA2s enzyme from Bothrops atrox (Fer-de-lance) had the most interactions of 35 compounds, followed by A8CG86 (PLA2s enzyme from Daboia russelii (Russell’s viper)) and P08878 (PLA2s enzyme of Crotalus durissus terrificus (South American rattlesnake)), which had five interacting compounds (Fig. 5). The prediction results only contain the PLA2 enzyme from various snake species, with each Uniprot Accession corresponding to a protein from a distinct snake species.
 | Figure 5. (a) Number of proteins the compounds interacted with, and (b) proteins with the most interaction with compounds. [Click here to view] |
3.2 Molecular docking
The 3D structures of PLA2 proteins from each organism were successfully retrieved from the PDB database, refined, and compared. The quality assessment and refinement results are summarized in Table 3. Based on the Ramachandran plot, a well-validated model should have over 90% of its residues in the most favored regions. Following refinement, all selected models met this criterion, ensuring their reliability for molecular docking simulation.
Table 3. Structural information of PLA2 proteins used in this study, including their source species, PDB ID, resolution, and percentage of residues in most favored regions.
| Organisms name | Uniprot accession | PDB ID | Resolution | % residues in most favoured regions | |
|---|---|---|---|---|---|
| Before refinement | After refinement | ||||
| Crotalus durissus terrificus | P08878 | 3R0L (Chain D) | 1.35 Å | 88.6% | 94.3% |
| Daboia russelii | A8CG86 | 1CL5 | 2.45 Å | 83.5% | 93.2% |
| Bothrops asper | Q6JK69 | 7M6C | 1.95 Å | 83.5% | 91.5% |
This study evaluated the binding affinities of selected compounds from C. cinereum against PLA2 through molecular docking simulations to assess their potential interactions with the target protein. Known antivenom drugs, varespladib and batimastat, were included as reference controls to provide a comparative framework for binding affinities and interaction effectiveness [72–74]. The docking score, measured in kcal/mol, represents the free energy change (ΔG) during ligand–protein interactions, where a more negative value indicates stronger and more stable binding, suggesting a higher likelihood of ligand efficacy in modulating PLA2 activity [75].
A total of 16 protein–ligand complexes were analyzed through molecular docking to assess their binding affinities and interactions with PLA2 from C. d. terrificus, D. russelii, and B. asper. Detailed interaction analyses and docking figures for the top three selected complexes are presented in the main text, while docking results for the remaining complexes, including binding affinities and key interacting residues, are summarized in Table 4 and provided in the supplementary materials (available from https://doi.org/10.6084/m9.figshare.29609597.v1).
Table 4. Molecular docking results.
| Organism name | Compound name | Binding affinity (kcal/mol) | Amino acid residues |
|---|---|---|---|
| Crotalus durissus terrificus | 2-Amino-3-(4-hydroxy-3-methoxyphenyl) propanoic acid | −6 | ASP48, ASN6, LEU2, HIS47, TRP30, GLY31, GLY29 |
| Varespladib | −6.3 | HIS47, GLY32, GLY31, ASP48, LEU2, TYR51, GLY29, LYS60 | |
| Batimastat | −5.8 | HIS47, GLY29, TRP30, GLY31, ASP48, CYS22, TYR27, GLY32, CYS49, GLY52, TYR51, LYS60, LEU2 | |
| Daboia russelii | trans-3-Indoleacrylic acid | −6.8 | ALA18, TYR22, PHE5, TRP31, GLY30, LYS69, ASP49, HIS48, TYR52, LEU2 |
| 4-Hydroxybenzoic acid | −5.2 | GLY30, TYR22, SER23, ALA18, LEU2, GLY6, PHE5, CYS29 | |
| Varespladib | −6.9 | ILE19, LEU2, TRP31, SER23, CYS29, GLY30, TYR22, PHE5, ALA18, ILE19 | |
| Batimastat | −6.2 | SER23, TRP31, ILE19, LEU2, GLY30, HIS48, TYR22, SER24 | |
| Bothrops asper | 2-hydroxyphenazine-1-carboxylic acid | −6.6 | TYR71, HIS67, ILE29, GLY49, PRO37, VAL38, GLY26, LEU22, CYS48, LYS68 |
| Phenylacetylene | −4.3 | CYS48, ILE29, HIS67, GLY26, VAL38, PRO37, LEU25 | |
| Sorbicillin | −5.9 | TYR41, ILE29, PRO37, GLY26, LEU25, LEU22, LYS80, HIS67, LYS68, TYR71, GLY49 | |
| 4-Coumarylalcohol | −4.8 | VAL38, GLY26, PRO37, ILE29, HIS67, LEU25, LEU22 | |
| trans-3-Indoleacrylic acid | −6 | TRY71, GLY49, LYS68, HIS67, TYR41, GLY26, LEU22, PRO37, VAL38 | |
| [6]-Gingerol | −5.6 | GLY49, CYS48, VAL112, LYS68, HIS67, VAL38, LEU22, TYR41, PRO37, GLY42, LEU25, GLY26, ILE29 | |
| Hydantocidin | −5.9 | LEU22, ILE29, CYS48, GLY42, TYR41, PRO37, HIS67, TYR71, GLY49, GLY26, LEU25 | |
| Varespladib | −6.9 | GLY49, HIS67, TYR71, LYS68, LEU25, VAL38, GLY26, PRO37, LEU22, LEU51, VAL50 | |
| Batimastat | −5.7 | LYS80, HIS67, GLY49, LYS68, TYR71, GLY26, PRO37, LEU25, GLY42, LEU22, VAL38, VAL50 |
Among the tested compounds, trans-3-indoleacrylic acid exhibited the highest binding affinity with PLA2 from D. russelii (−6.8 kcal/mol). The interaction analysis revealed that this compound formed hydrophobic interactions with key residues, including ALA18, TYR22, PHE5, TRP31, GLY30, LYS69, ASP49, HIS48, TYR52, and LEU2 (Fig. 6e). The presence of multiple hydrophobic residues, such as ALA18, TYR22, and PHE5, suggests that trans-3-indoleacrylic acid fits well within the hydrophobic pocket of PLA2, stabilizing the binding through van der Waals forces [76]. Additionally, hydrogen bonds with GLY30 and ASP49 further reinforce the interaction, enhancing binding stability and specificity. The hydrogen bond with ASP49 is particularly significant as this residue plays a crucial role in the enzyme’s catalytic mechanism, potentially contributing to the inhibitory effect of the ligand [77].
 | Figure 6. Molecular docking visualization and protein-ligand interaction of PLA2 (C. d. terrificus), PLA2 (D. russelli), and PLA2 (B. asper) with the top three selected complexes and varespladib. [Click here to view] |
Meanwhile, 2-amino-3-(4-hydroxy-3-methoxyphenyl) propanoic acid exhibited a binding affinity of −6.0 kcal/mol with PLA2 from C. d. terrificus. The docking analysis revealed interactions with key residues, including ASP48, ASN6, LEU2, HIS47, TRP30, GLY31, and GLY29 (Fig. 6a). Notably, hydrogen bonds were formed with ASP48 and ASN6, which likely contribute to the stability and specificity of the binding. Similarly, 2-hydroxyphenazine-1-carboxylic acid displayed a binding affinity of −6.6 kcal/mol with PLA2 from B. asper. The ligand interacted with key residues, including TYR71, HIS67, ILE29, GLY49, PRO37, VAL38, GLY26, LEU22, CYS48, and LYS68 (Fig. 6c). A hydrogen bond with TYR71 was observed, which likely strengthens the stability of the complex. Additionally, interactions with hydrophobic residues such as ILE29, LEU22, and VAL38 suggest van der Waals contributions that further reinforce the ligand’s binding.
Comparative analysis across the three PLA2 isoforms revealed that varespladib, a well-established antivenom drug, exhibited the strongest binding affinity (−6.9 kcal/mol), likely due to its specific role as a PLA2 inhibitor designed to target and neutralize the enzyme’s activity. This role also reflects the overall conservation of the PLA2 binding pockets despite minor structural differences. Additionally, its binding interactions provided valuable insights into the binding pocket, as molecular docking results revealed that compounds from C. cinereum and varespladib interacted with several of the same amino acid residues. This finding is consistent with previous studies, further supporting varespladib’s potential as a reference inhibitor for evaluating new PLA2-targeting compounds [32,78]. However, compounds from C. cinereum displayed variable affinities depending on the isoform, indicating that subtle structural differences among PLA2 proteins [79].
While molecular docking provides an initial prediction of ligand–protein binding affinity, it is limited by its static nature and simplifications such as fixed protonation states and tautomeric forms. In this study, ligands were prepared at physiological pH to assign the dominant protonation state, but multiple tautomeric forms were not explicitly considered. To address this limitation and better capture the dynamic behavior of ligand–protein interactions, MD simulations were performed.
3.3 MD simulation
3.3.1. RMSD
The structural stability of the protein–ligand complexes was assessed using RMSD analysis over a 100 ns MD simulation (Fig. 7). For PLA2 (C. d. terrificus)–2-Amino-3-(4-hydroxy-3-methoxyphenyl) propanoic acid complex, the RMSD trajectory exhibited fluctuations, with a marked increase in instability between 60 and 80 ns, reaching approximately 7 Å. The ligand RMSD remained around 2 Å initially but gradually decreased to ~1 Å by 90 ns, suggesting improved stabilization over time.
PLA2 (D. russelii)–trans-3-Indoleacrylic acid demonstrated less difference between the protein RMSD and the complex RMSD compared to the other systems, indicating minimal structural deviation. The complex maintained the stability below 5 Å until approximately 80 ns, after which a sudden rise in RMSD was observed. The ligand remained stable for the initial phase of the simulation but showed fluctuations beginning around 35 ns and persisting until the end of the trajectory. The cause of this fluctuation remains unclear and may be attributed to ligand conformational rearrangements or transient interactions with the binding site [80].
 | Figure 7. Structural stability and dynamics of three protein–ligand complexes analyzed through (A) RMSD, (B) RMSF, (C) radius of gyration, and (D) hydrogen bond profiles over the simulation time. [Click here to view] |
PLA2 (B. asper) –2-hydroxyphenazine-1-carboxylic acid displayed a larger deviation between the complex and protein RMSD values relative to the other systems. However, the overall RMSD remained below 5 Å throughout the 100 ns simulation, suggesting structural stability. The ligand RMSD remained consistently below 1 Å, indicating a stable binding conformation. Despite PLA2 (D. russelii) –trans-3-Indoleacrylic acid complex exhibiting the lowest binding energy, its MD simulation results suggest reduced stability compared to the other complexes. Additionally, all complexes showed an upward trend in RMSD toward the end of the simulation, implying that a longer simulation may be required to ensure equilibration and assess long-term stability [81,82].
3.3.2 RMSF
The RMSF was calculated to assess the flexibility of individual residues in each protein–ligand complex over the 100 ns MD simulation. Higher RMSF values indicate regions of greater structural flexibility, often corresponding to loop regions, while lower values suggest more stable secondary structures such as α-helices and β-sheets [83].
As shown in Figure 7, PLA2 (C. d. terrificus)–2-Amino-3-(4-hydroxy-3-methoxyphenyl) propanoic acid complex exhibited notable fluctuations around residues 55 and 70, suggesting increased flexibility in these regions. PLA2 (D. russelii) –trans-3-Indoleacrylic acid displayed peaks at residues 30–40, 70, and 80, indicating multiple regions of enhanced mobility. Both complexes showed significantly high fluctuations at the terminal regions, with RMSF values approaching 11 Å, whereas the fluctuations in the remaining residues remained below 3 Å. The fluctuations suggest that the C-terminal regions are highly flexible, as expected for regions that are often solvent-exposed and lack strong structural constraints [84].
PLA2 (B. asper) –2-hydroxyphenazine-1-carboxylic acid followed a similar RMSF pattern to the other two complexes but exhibited a distinct difference in the terminal region, where the fluctuation was limited to approximately 3 Å instead of 11 Å. The fluctuations indicate that this complex may exhibit increased stability at its termini compared to the other two complexes. The differences in RMSF profiles between the complexes may be attributed to slight variations in residue numbering, which could affect how fluctuations are mapped across the sequence.
3.3.3 Rg
The Rg was analyzed to assess the compactness and structural stability of the protein–ligand complexes over the 100 ns MD simulation [85]. Figure 7 presents the Rg trajectories for all three complexes, showing a similar overall pattern between the protein and its respective complex, with the complex consistently exhibiting slightly higher Rg values. However, in PLA2 (B. asper)–2-hydroxyphenazine-1-carboxylic acid, the difference between the complex and protein is more pronounced, indicating a distinct structural behavior.
For PLA2 (C. d. terrificus)–2-Amino-3-(4-hydroxy-3-methoxyphenyl) propanoic acid complex, the Rg initially starts at approximately 15 Å and decreases within the first 25 ns. However, a subsequent rise is observed, with the highest peak reaching around 16.5 Å, suggesting structural fluctuations that may indicate partial unfolding or conformational changes [85]. PLA2 (D. russelii) – trans-3-Indoleacrylic acid complex demonstrates a more stable trajectory compared to PLA2 (C. d. terrificus)–2-Amino-3-(4-hydroxy-3-methoxyphenyl) propanoic acid complex, though a noticeable rise in Rg is observed, peaking around 80 ns. This increase in Rg suggests a potential loosening of the protein–ligand complex, which may correlate with the observed RMSD fluctuations. PLA2 (B. asper)–2-hydroxyphenazine-1-carboxylic acid complex exhibits the most stable Rg profile, maintaining a consistent value around 15 Å throughout the entire simulation, indicating a more rigid and compact structure compared to the other two complexes. The lower fluctuations suggest that the ligand binding in PLA2 (B. asper)–2-hydroxyphenazine-1-carboxylic acid provides a stabilizing effect, contributing to overall structural integrity.
3.3.4 Hydrogen bonds
Hydrogen bonding (H-bond) interactions are essential for the stability of protein–ligand complexes. These interactions contribute to molecular recognition by facilitating the correct positioning of ligands within the binding pocket. The number of H-bonds formed during the 100 ns MD simulation was analyzed for each complex (Fig. 7). PLA2 (C. d. terrificus)–2-Amino-3-(4-hydroxy-3-methoxyphenyl) propanoic acid complex exhibited the highest number of H-bonds, ranging from 1 to 7 throughout the simulation. For most of the trajectory, the complex maintained between 3 to 5 H-bonds, indicating a relatively stable interaction between the ligand and the protein.
In contrast, PLA2 (D. russelii)–trans-3-Indoleacrylic acid complex showed a decline in H-bond formation over time. During the first 30 ns, the complex maintained 2 to 3 H-bonds, but after this period, the number of H-bonds decreased, with only one H-bond persisting for the remainder of the simulation. This pattern indicates a potential weakening of the ligand–protein interaction, possibly due to ligand repositioning, loss of key contacts, or partial dissociation from the active site. Combined with RMSD fluctuations and fewer contacts observed in the final frames, this suggests a dynamic instability in the complex.
PLA2 (B. asper)–2-hydroxyphenazine-1-carboxylic acid complex exhibited the lowest number of H-bonds. Notably, during the initial 10–30 ns, no hydrogen bonds were formed. After this period, the complex showed minimal H-bond formation, with only 1–2 H-bonds observed, reaching a maximum of 3 at isolated time points. This lower frequency of H-bond indicates a weaker interaction between the ligand and the protein compared to the other two complexes [77]. However, this may be compensated by other non-polar interactions, such as π–π stacking or van der Waals forces, which should be explored further in future analyses [86].
3.3.5 Molecular conformational changes before and after MD
The Matchmaker tool in UCSF Chimera was used to compare the conformations of the small molecules before and after the MD simulation (Fig. 8). In the visualization, the orange structure represents the docked conformation, while the blue structure corresponds to the conformation after 100 ns of simulation. The comparison showed a high degree of consistency between the two, indicating that the binding sites remained stable and were not significantly altered during the dynamics.
 | Figure 8. Molecular conformational changes before and after (A) PLA2 (C. d. terrificus)–2-Amino-3-(4-hydroxy-3-methoxyphenyl) propanoic acid; (B) PLA2 (D. russelii) –trans-3-Indoleacrylic acid; (C) PLA2 (B. asper) –2-hydroxyphenazine-1-carboxylic acidE. [Click here to view] |
3.3.6 MM/GBSA calculation
Binding free energy calculations were conducted using the gmx_MMPBSA tool, accounting for both enthalpic and entropic contributions to ligand binding with the proteins. The results of the end-state MMGBSA binding free energy analysis are presented in Table 5.
Table 5. Energy components of the ligand-protein complexes obtained from MMGBSA calculations, including van der Waals energy (ΔVDWAALS), electrostatic energy (ΔEEL), polar solvation energy (ΔEGB), non-polar solvation energy (ΔESURF), gas-phase energy (ΔGGAS), solvation energy (ΔGSOLV), and total binding energy (ΔTOTAL).
| Complex | Energy (Kcal/mol) | |||||||
|---|---|---|---|---|---|---|---|---|
| Protein | Ligand | ΔVDWAALS | ΔEEL | ΔEGB | ΔESURF | ΔGGAS | ΔGSOLV | ΔTOTAL |
| PLA2 (C. d. terrificus) | 2-Amino-3-(4-hydroxy- 3-methoxyphenyl)propanoic acid | −6.26 | −104.84 | 87.54 | −2.81 | −111.09 | 84.73 | −26.37 |
| PLA2 (D. russelii) | trans-3-Indoleacrylic acid | −22.79 | −19.57 | 28.15 | −3.34 | −42.36 | 24.80 | −17.56 |
| PLA2 (B. asper) | 2-hydroxyphenazine-1-carboxylic acid | −28.90 | −5.63 | 17.23 | −3.33 | −34.54 | 13.90 | −20.63 |
The binding free energy analysis using the MM/GBSA method revealed that all three PLA2–ligand complexes exhibited favorable interactions, with negative ΔTOTAL values indicating stable binding. Among them, the PLA2 from Crotalus durissus terrificus complexed with 2-Amino-3-(4-hydroxy-3-methoxyphenyl) propanoic acid showed the most favorable binding energy (ΔG_total = –26.37 kcal/mol), suggesting it forms the most stable complex. This strong affinity is primarily driven by a significant electrostatic component (ΔEEL = –104.84 kcal/mol), although it was primarily offset by high polar solvation energy (ΔEGB = 87.54 kcal/mol). The PLA2 (B. asper) complex with 2-hydroxyphenazine-1-carboxylic acid also showed a strong binding energy (ΔTOTAL = –20.63 kcal/mol), primarily driven by strong van der Waals interactions (ΔVDWAALS = –28.90 kcal/mol). Meanwhile, the D. russelii PLA2 complex with trans-3-indoleacrylic acid also demonstrated stable binding (ΔTOTAL = –17.56 kcal/mol), supported by both van der Waals and electrostatic contributions.
3.4 Molecular insights into PLA2 inhibition by C. cinereum-derived compounds
From a biological standpoint, the exclusive identification of specific PLA2 isoforms as targets among the predicted compound-protein interactions is highly significant. PLA2s are known to be among the most conserved and abundantly expressed components across viperid and elapid venoms, playing central roles in inducing hemotoxic, myotoxic, and inflammatory effects [14,16]. Their catalytic activity on membrane phospholipids leads to rapid disruption of cellular integrity, hemorrhage, and organ dysfunction following envenomation. Therefore, the identification of C. cinereum-derived compounds that selectively bind to PLA2 enzymes aligns with a biologically relevant therapeutic strategy aimed at neutralizing one of the primary drivers of venom toxicity.
Moreover, the specific residues involved in compound binding, particularly His48 and Asp49, are critical to the catalytic mechanism of snake venom PLA2 enzymes. These residues are responsible for coordinating the catalytic site and facilitating the nucleophilic attack on phospholipid substrates. Binding interactions that involve H-bond or hydrophobic stabilization at these key sites may lead to competitive inhibition or allosteric modulation of enzymatic activity, thus mitigating the cytotoxic effects of venom [15,17].
The MD simulation results provide further insights into the dynamic stability of these interactions under physiologically relevant conditions. Stable binding, reflected in lower RMSD and consistent hydrogen bond formation, suggests that compounds like 2-Amino-3-(4-hydroxy-3-methoxyphenyl) propanoic acid could effectively maintain inhibitory interactions over time in vivo. The consistency is crucial, as venom-induced pathologies progress rapidly after envenomation, requiring inhibitors that can act promptly and remain bound under fluctuating biological conditions [87].
Interestingly, the differences observed between the complexes, such as the distinct binding patterns of our compounds across the three PLA2 isoforms and the higher flexibility of PLA2 (D. russelii), may correlate with species-specific differences in venom composition and toxicity profiles. This suggests that the compounds may exhibit isoform-specific inhibitory effects rather than a uniform, broad-spectrum efficacy. It has been reported that minor sequence variations among PLA2 isoforms can significantly affect their enzymatic properties and pathological effects, which may in turn influence the binding dynamics of small molecule inhibitors [79].
Thus, the findings of this study not only validate the computational predictions but also support the biological plausibility of C. cinereum-derived compounds as potential targeted inhibitors of specific venom PLA2 enzymes. Future experimental studies, including enzymatic inhibition assays and in vivo neutralization experiments, are warranted to confirm the therapeutic potential of these candidate compounds.
4. CONCLUSION
This study aimed to identify potential inhibitors of the toxic PLA2 enzyme from the plant C. cinereum using an integrated computational approach. By combining DTI prediction, molecular docking, and MD simulation, several phytochemicals were identified showing strong binding affinity and stable interaction with at least one of the PLA2 isoforms tested. This isoform-specific binding highlights their potential as starting points for the development of targeted anti-snake venom drugs. While the known PLA2 inhibitor varespladib demonstrated a higher binding affinity from molecular docking results, our findings highlight that some C. cinereum phytochemicals exhibited comparable binding scores and favorable MD stability. The focus on these plant-derived compounds is driven by their natural origin, structural diversity, and potential for safer and more sustainable drug discovery.
This computational framework proved effective in screening and characterizing bioactive compounds before experimental testing. Nevertheless, the current findings are limited to in silico analysis; therefore, further validation through in vitro and in vivo assays is essential to confirm the therapeutic efficacy and safety of these compounds against a broader range of PLA2s. Future research should also explore pharmacokinetic and toxicity profiles to advance these phytochemicals towards clinical applications.
5. ACKNOWLEDGMENTS
The authors are grateful to the Advanced Research Laboratory, IPB University, for the HPC (NVIDIA DGX A100) cluster.
6. AUTHOR CONTRIBUTIONS
All authors made substantial contributions to conception and design, acquisition of data, or analysis and interpretation of data; took part in drafting the article or revising it critically for important intellectual content; agreed to submit to the current journal; gave final approval of the version to be published; and agree to be accountable for all aspects of the work. All the authors are eligible to be an author as per the International Committee of Medical Journal Editors (ICMJE) requirements/guidelines.
7. FINANCIAL SUPPORT
This work was supported by the Directorate of Research, Technology and Community Service, Ministry of Education, Culture, Research and Technology, through the Strategic Research Collaboration Program (KATALIS), Grant No. 001/E5/PG.02.00/PL.BATCH.2/2024.
8. CONFLICTS OF INTEREST
The authors report no financial or any other conflicts of interest in this work.
9. ETHICAL APPROVALS
This study does not involve experiments on animals or human subjects.
10. DATA AVAILABILITY
All data supporting the findings of this study, including detailed docking results and supplementary figures, are available on Figshare at https://doi.org/10.6084/m9.figshare.29609597.v1.
11. PUBLISHER’S NOTE
All claims expressed in this article are solely those of the authors and do not necessarily represent those of the publisher, the editors and the reviewers. This journal remains neutral with regard to jurisdictional claims in published institutional affiliation.
12. USE OF ARTIFICIAL INTELLIGENCE (AI)-ASSISTED TECHNOLOGY
The authors declare that they have not used artificial intelligence (AI)-tools for writing and editing of the manuscript, and no images were manipulated using AI.
REFERENCES
1. Afroz A, Siddiquea BN, Chowdhury HA, Jackson TN, Watt AD. Snakebite envenoming: a systematic review and meta-analysis of global morbidity and mortality. PLoS Negl Trop Dis. 2024;18(4):12080. CrossRef
2. Kasturiratne A, Wickremasinghe AR, De Silva N, Gunawardena NK, Pathmeswaran A, Premaratna R, et al. The global burden of snakebite: a literature analysis and modelling based on regional estimates of envenoming and deaths. PLoS Med. 2008;5(11):1591–604. CrossRef
3. Martin G, Yanez-Arenas C, Rangel-Camacho R, Murray KA, Goldstein E, Iwamura T, et al. Implications of global environmental change for the burden of snakebite. Toxicon X. 2021;9:100069. CrossRef
4. Gutiérrez JM, Calvete JJ, Habib AG, Harrison RA, Williams DJ, Warrell DA. Snakebite envenoming. Nat Rev Dis Primers. 2017;3(1):1–21. CrossRef
5. Cruz LS, Vargas R, Lopes AA. Snakebite envenomation and death in the developing world. Ethn Dis. 2009;19(1 SUPPL.1):42–6. Available from: https://www.jstor.org/stable/48667815
6. Sofyantoro F, Yudha DS, Lischer K, Nuringtyas TR, Putri WA, Kusuma WA, et al. Bibliometric analysis of literature in snake venom-related research worldwide (1933–2022). Animals. 2022;12(16):2058. CrossRef
7. Gomez-Betancur I, Gogineni V, Salazar-Ospina A, Leon F. Perspective on the therapeutics of anti-snake venom. Molecules. 2019;24(18):3276. CrossRef
8. Puzari U, Fernandes PA, Mukherjee AK. Advances in the therapeutic application of small-molecule inhibitors and repurposed drugs against snakebite: mini perspective. J Med Chem. 2021;64(19):13938–79. CrossRef
9. Alomran N, Chinnappan R, Alsolaiss J, Casewell NR, Zourob M. Exploring the utility of ssDNA aptamers directed against snake venom toxins as new therapeutics for snakebite envenoming. Toxins (Basel). 2022;14(7):469. CrossRef
10. Chanda A, Salvi NC, Shelke PV, Kalita B, Patra A, Puzari U, et al. Supplementation of polyclonal antibodies, developed against epitope-string toxin-specific peptide immunogens, to commercial polyvalent antivenom, shows improved neutralization of Indian Big Four and Naja kaouthia snake venoms. Toxicon X. 2024;24:100210. CrossRef
11. Silva A, Isbister GK. Current research into snake antivenoms, their mechanisms of action and applications. Biochem Soc Trans. 2020;48(2):537–46. CrossRef
12. Khalek IS, Senji Laxme RR, Nguyen YTK, Khochare S, Patel RN, Woehl J, et al. Synthetic development of a broadly neutralizing antibody against snake venom long-chain Α-neurotoxins. Sci Transl Med. 2024;16(735):1867. CrossRef
13. Aphrodita A, Sentono DN, Yudha DS, Purwestri YA, Nuringtyas TR, Raharjo S, et al. Comparative analysis of hemotoxic, myotoxic, and inflammatory profiles of Calloselasma rhodostoma and Trimeresurus insularis venoms in mice. Narra J. 2025;5(2):e1874. CrossRef
14. Borges RJ, Salvador GHM, Campanelli HB, Pimenta DC, De Oliveira Neto M, Usón I, et al. BthTX-II from Bothrops jararacussu venom has variants with different oligomeric assemblies: an example of snake venom phospholipases A2 versatility. Int J Biol Macromol. 2021;191:255–66. CrossRef
15. De Oliveira LL, Persinoti GF, Giuliatti S, Tinós R. Feature selection via genetic algorithms in the classification of anti-snake venom medicinal plants. J Braz Neural Netw Soc. 2010;8(3):125–34. Available from: https://sbic.org.br/wp-content/uploads/sites/4/2016/07/vol8-no3-art1.pdf
16. Hiu JJ, Yap MKK. Cytotoxicity of snake venom enzymatic toxins: phospholipase A2 and L-amino acid oxidase. Biochem Soc Trans. 2020;48(2):719–31. CrossRef
17. Peggion C, Tonello F. Short linear motifs characterizing snake venom and mammalian phospholipases A2. Toxins (Basel). 2021;13(4):290. CrossRef
18. Tsai IH, Liu HC, Chang T. Toxicity domain in presynaptically toxic phospholipase A2 of snake venom. Protein Struct Mol Enzymol. 1987;916(1):94–9. CrossRef
19. Bala AA, Mohammed M, Umar S, Ungogo MA, Al-kassim Hassan M, Abdussalam US, et al. Pre-clinical efficacy of African medicinal plants used in the treatment of snakebite envenoming: a systematic review. Toxicon. 2023;224:20499361211072644. CrossRef
20. Dharmadasa RM, Akalanka GC, Muthukumarana PRM, Wijesekara RGS. Ethnopharmacological survey on medicinal plants used in snakebite treatments in Western and Sabaragamuwa provinces in Sri Lanka. J Ethnopharmacol. 2016;179:110–27. CrossRef
21. Giovannini P, Howes MJR. Medicinal plants used to treat snakebite in Central America: review and assessment of scientific evidence. J Ethnopharmacol. 2017;199:240–56. CrossRef
22. Adinortey MB. Botanical treatments for snakebite in rural Ghana: a narrative review. J Ethnopharmacol. 2021;280:114432. CrossRef
23. Mogha NG, Kalokora OJ, Amir HM, Kacholi DS. Ethnomedicinal plants used for treatment of snakebites in Tanzania–a systematic review. Pharm Biol. 2022;60(1):1925–34. CrossRef
24. Upasani SV, Beldar VG, Tatiya AU, Upasani MS, Surana SJ, Patil DS. Ethnomedicinal plants used for snakebite in India: a brief overview. Integr Med Res. 2017;6(2):114–30. CrossRef
25. Adrião AAX, Dos Santos AO, De Lima EJSP, Maciel JB, Paz WHP, Da Silva FMA, et al. Plant-derived toxin inhibitors as potential candidates to complement antivenom treatment in snakebite envenomations. Front Immunol. 2022;13:842576. CrossRef
26. Puzari U, Fernandes PA, Mukherjee AK. Pharmacological re- assessment of traditional medicinal plants-derived inhibitors as antidotes against snakebite envenoming: a critical review. J Ethnopharmacol. 2022;292:115208. CrossRef
27. Ravi DA, Hwang DH, Mohan Prakash RL, Kang C, Kim E. Indian medicinal plant-derived phytochemicals as potential antidotes for snakebite: a pharmacoinformatic study of Atrolysin inhibitors. Int J Mol Sci. 2024;25(23):12675. CrossRef
28. Byeon JC, Ahn JB, Jang WS, Lee SE, Choi JS, Park JS. Recent formulation approaches to oral delivery of herbal medicines. J Pharm Investig. 2019;49(1):17–26. CrossRef
29. Alangode A, Rajan K, Nair BG. Snake antivenom: challenges and alternate approaches. Biochem Pharmacol. 2020;181:114135. CrossRef
30. Devy DKV, Ruby S, Kumar M, Kumar PNS, Vijayakumar B. Herbal antivenins as an eye opener towards snake envenomation. Int J Advancement Life Sci Res. 2024;7(4):1–1. CrossRef
31. Suji S, Dinesh M, Keerthi K, Anagha K, Arya J, Anju K. Evaluation of neutralization potential of Naja naja and Daboia russelii snake venom by root extract of Cyanthillium cinereum. Indian J Crit Care Med. 2023;27(11):821–9. CrossRef
32. Yusuf AJ, Bugaje AI, Sadiq M, Salihu M, Adamu HW, Abdulrahman M. Exploring the inhibitory potential of phytochemicals from Vernonia glaberrima leaves against snake venom toxins through computational simulation and experimental validation. Toxicon. 2024;247:107838. CrossRef
33. Kusuma WA, Fadli A, Fatriani R, Sofyantoro F, Yudha DS, Lischer K, et al. Prediction of the interaction between Calloselasma rhodostoma venom-derived peptides and cancer-associated hub proteins: a computational study [Internet]. Heliyon Internet. 2023;9(11):e21149. CrossRef
34. Snake Envenoming Medicines Database. Impact Global Health; [cited 2024 Dec 27]. Available from: https://www.impactglobalhealth.org/data/snakebite-envenoming-medicines-database/dashboard
35. PubChem Substructure Fingerprint V1.3 [Internet]. [cited 2021 Jul 14].
36. Morgan HL. The Generation of a Unique Machine Description for Chemical Structures—A Technique Developed at Chemical Abstracts Service [Internet]. J Chem Doc Internet. 1965;5(2):107–13. CrossRef
37. Lin Z, Akin H, Rao R, Hie B, Zhu Z, Lu W, et al. Evolutionary-scale prediction of atomic-level protein structure with a language model [Internet]. Sci (1979). 2023;379(6637):1123–30. CrossRef
38. Hähnke VD, Kim S, Bolton EE. PubChem chemical structure standardization [Internet]. J Cheminform Internet. 2018;10(1):1–40. CrossRef
39. Cereto-Massagué A, Ojeda MJ, Valls C, Mulero M, Garcia-Vallvé S, Pujadas G. Molecular fingerprint similarity search in virtual screening. Methods. 2015;71(C):58–63. CrossRef
40. Fadli A, Kusuma WA, Annisa, Batubara I, Heryanto R. Screening of potential indonesia herbal compounds based on multi-label classification for 2019 coronavirus disease [Internet]. Big Data Cognit Comput. 2021;5(4):75. CrossRef
41. He S, Zhao D, Ling Y, Cai H, Cai Y, Zhang J, et al. Machine learning enables accurate and rapid prediction of active molecules against breast cancer cells [Internet]. Front Pharmacol Internet. 2021;12:796534. CrossRef
42. Capecchi A, Probst D, Reymond JL. One molecular fingerprint to rule them all: drugs, biomolecules, and the metabolome [Internet]. J Cheminform Internet. 2020;12(1):1–5. CrossRef
43. Chou KC. A novel approach to predicting protein structural classes in a (20–1)-D amino acid composition space [Internet]. Proteins Struct Function Bioinf. 1995;21(4):319–44. CrossRef
44. Bahi M, Batouche M. Drug-target interaction prediction in drug repositioning based on deep semi-supervised learning. In: Amine A, Mouhoub M, Ait Mohamed O, Djebbar B, editors. Computational intelligence and its applications. Cham: Springer; 2018. Vol. 522. pp. 309–19. CrossRef
45. Qi Y. Random forest for bioinformatics. In: Zhang C, Ma Y, editors. Ensemble machine learning. New York, NY: Springer; 2012. pp. 307–23. CrossRef
46. Liu X, Zhao L, Dong Q. Protein remote homology detection based on auto-cross covariance transformation. Comput Biol Med. 2011;41(8):640–7. CrossRef
47. Shahriari B, Swersky K, Wang Z, Adams RP, De Freitas N. Taking the human out of the loop: a review of Bayesian optimization. Proc IEEE. 2016;104(1):148–75. CrossRef
48. LaValle SM, Branicky MS. On the relationship between classical grid search and probabilistic roadmaps. Springer Tracts Adv Rob. 2004;7 STAR(7–8):59–75. CrossRef
49. Blagus R, Lusa L. SMOTE for high-dimensional class-imbalanced data [Internet]. BMC Bioinf. 2013;14(1):1–6. CrossRef
50. Krawczyk B. Learning from imbalanced data: open challenges and future directions [Internet]. Prog Artif Intell. 2016;5(4):221–32. CrossRef
51. Sulistiawan F, Kusuma WA, Ramadhanti NS, Tedjo A. Drug-target interaction prediction in coronavirus disease 2019 case using deep semi-supervised learning model. In: Proceedings of the 2020 International Conference on Advanced Computer Science and Information Systems (ICACSIS); 2020 Oct; Depok, Indonesia. Piscataway, NJ: IEEE; 2020. pp. 83–8. CrossRef
52. Brodersen KH, Ong CS, Stephan KE, Buhmann JM. The binormal assumption on precision-recall curves. In: Proceedings of the 20th International Conference on Pattern Recognition (ICPR); 2010 Aug; Istanbul, Turkey. Piscataway, NJ: IEEE; 2010. pp. 4263–6. CrossRef
53. Jain AN, Nicholls A. Recommendations for evaluation of computational methods [Internet]. J Comput Aided Mol Des. 2008;22(3–4):133–9. CrossRef
54. Tanoli Z, Schulman A, Aittokallio T. Validation guidelines for drug-target prediction methods [Internet]. Expert Opin Drug Discov. 2025;20(1):31–45. CrossRef
55. Berman HM, Westbrook J, Feng Z, Gilliland G, Bhat TN, Weissig H, et al. The protein data bank. Nucleic Acids Res. 2000;28(1):235–42. CrossRef
56. Laskowski RA, Macarthur MW, Moss DS, Thornton JM. PROCHECK: a program to check the stereochemical quality of protein structures. J Appl Crystallogr. 1993;26(2):283–91. CrossRef
57. Lee GR, Won J, Heo L, Seok C. GalaxyRefine2: simultaneous refinement of inaccurate local regions and overall protein structure. Nucleic Acids Res. 2019;47(W1):W451–5. CrossRef
58. Pettersen EF, Goddard TD, Huang CC, Couch GS, Greenblatt DM, Meng EC, et al. UCSF Chimera - A visualization system for exploratory research and analysis. J Comput Chem. 2004;25(13):1605–2. CrossRef
59. Kim S, Chen J, Cheng T, Gindulyte A, He J, He S, et al. PubChem 2023 update. Nucleic Acids Res. 2023;51(1):D1373–1380. CrossRef
60. O’Boyle NM, Banck M, James CA, Morley C, Vandermeersch T, Hutchison GR. Open Babel: an Open chemical toolbox. J Cheminform. 2011;3(10):1–4. CrossRef
61. Eberhardt J, Santos-Martins D, Tillack AF, Forli S. AutoDock Vina 1.2.0: new docking methods, expanded force field, and python bindings. J Chem Inf Model. 2021;61(8):3891–8. CrossRef
62. Trott O, Olson AJ. AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J Comput Chem. 2010;31(2):455–61. CrossRef
63. Wallace AC, Laskowski RA, Thornton JM. Ligplot: a program to generate schematic diagrams of protein-ligand interactions. Protein Eng Design Selection. 1995;8(2):127–34. CrossRef
64. DeLano WL. Pymol: an open-source molecular graphics tool. {CCP4} Newsletter On Protein Crystallogr [Internet]. 2002;40(1):1–8.
65. Abraham M, Alekseenko A, Basov V, Bergh C, Briand E, Brown A, et al. GROMACS 2024.2 Manual [Internet]. Zenodo. 2024 [cited 2024 Dec 27]. Available from: https://doi.org/10.5281/zenodo.11148638
66. Abraham MJ, Murtola T, Schulz R, Páll S, Smith JC, Hess B, et al. Gromacs: high performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX. 2015;1–2(2):19–25. CrossRef
67. Brooks BR, Brooks CL, Mackerell AD, Nilsson L, Petrella RJ, Roux B, et al. CHARMM: the biomolecular simulation program. J Comput Chem. 2009;30(10):1545–6414. CrossRef
68. Huang J, Mackerell AD. CHARMM36 all-atom additive protein force field: validation based on comparison to NMR data [Internet]. J Comput Chem Internet. 2013;34(25):2135–45. CrossRef
69. Zoete V, Cuendet MA, Grosdidier A, Michielin O. SwissParam: a fast force field generation tool for small organic molecules [Internet]. J Comput Chem. 2011;32(11):2359–68. CrossRef
70. Bondi A. Van der Waals Volumes and Radii [Internet]. J Phys Chem. 1964;68(3):441–51. CrossRef
71. Valdés-Tresanco MS, Valdés-Tresanco ME, Valiente PA, Moreno E. Gmx_MMPBSA: a new tool to perform end-state free energy calculations with GROMACS [Internet]. J Chem Theory Comput. 2021;17(10):6281–91. CrossRef
72. Quiroz S, Henao Castañeda IC, Granados J, Patiño AC, Preciado LM, Pereañez JA. Inhibitory effects of varespladib, CP471474, and their potential synergistic activity on Bothrops asper and Crotalus durissus cumanensis Venoms. Molecules. 2022;27(23):8588. CrossRef
73. Lewin M, Samuel S, Merkel J, Bickler P. Varespladib (LY315920) appears to be a potent, broad-spectrum, inhibitor of snake venom phospholipase A2 and a possible pre-referral treatment for envenomation. Toxins (Basel). 2016;8(9):248. CrossRef
74. Escalante T, Franceschi A, Rucavado A, Gutiérrez JM. Effectiveness of batimastat, a synthetic inhibitor of matrix metalloproteinases, in neutralizing local tissue damage induced by BaP1, a hemorrhagic metalloproteinase from the venom of the snake Bothrops asper. Biochem Pharmacol. 2000;60(2):269–74. CrossRef
75. Luque I, Freire E. Structural stability of binding sites: consequences for binding affinity and allosteric effects. Proteins Struct Function Genet. 2000;41(2000):63–71. CrossRef
76. Barratt E, Bingham RJ, Warner DJ, Laughton CA, Phillips SEV, Homans SW. Van der Waals interactions dominate ligand-protein association in a protein binding site occluded from solvent water. J Am Chem Soc. 2005;127(33):11827–234. CrossRef
77. Bulusu G, Desiraju GR. Strong and weak hydrogen bonds in protein–ligand recognition. J Indian Inst Sci. 2020;100(1):31–41. CrossRef
78. Mena-Ulecia K, Tiznado W, Caballero J. Study of the differential activity of thrombin inhibitors using docking, QSAR, molecular dynamics, and MM-GBSA. PLoS One. 2015;10(11):142774. CrossRef
79. Manjunatha Kini R. Excitement ahead: structure, function and mechanism of snake venom phospholipase A2 enzymes. Toxicon. 2003;42(8):827–40. CrossRef
80. Ahmad E, Rabbani G, Zaidi N, Khan MA, Qadeer A, Ishtikhar M, et al. Revisiting ligand-induced conformational changes in proteins: essence, advancements, implications and future challenges. J Biomol Struct Dyn. 2013;31(6):630–48. CrossRef
81. Hollingsworth SA, Dror RO. Molecular dynamics simulation for all. Neuron. 2018;99(6):1129–43. CrossRef
82. Saurabh S, Sivakumar PM, Perumal V, Khosravi A, Sugumaran A, Prabhawathi V. Molecular dynamics simulations in drug discovery and drug delivery. Eng Mater. 2020;9(1):275–301. CrossRef
83. Martínez L. Automatic identification of mobile and rigid substructures in molecular dynamics simulations and fractional structural fluctuation analysis. PLoS One. 2015;10(3):119264. CrossRef
84. Hubbell WL, Altenbach C, Hubbell CM, Khorana HG. Rhodopsin structure, dynamics, and activation: a perspective from crystallography, site-directed spin labeling, sulfhydryl reactivity, and disulfide cross-linking. Adv Protein Chem. 2003;63:243–90. CrossRef
85. Lobanov M, Bogatyreva NS, Galzitskaia OV. Radius of gyration is indicator of compactness of protein structure. Molekuliarnaia Biologiia. 2008;42(4):701–6.
86. Patil R, Das S, Stanley A, Yadav L, Sudhakar A, Varma AK. Optimized hydrophobic interactions and hydrogen bonding at the target-ligand interface leads the pathways of Drug-Designing. PLoS One. 2010;5(8):e12029. CrossRef
87. Lomonte B, Rangel J. Snake venom Lys49 myotoxins: from phospholipases A2 to non-enzymatic membrane disruptors. Toxicon. 2012;60(4):520–30. CrossRef