In Silico Screening of Potential Inhibitors against dPLA2 from Named Chinese Herbs for Identification of Compounds with Antivenom Effects Due to Deinagkistrodon acutus Snake Bites

Xingyang Xiao¹, Mengyi Lai¹, Zechang Rao¹, Jianzhong Huang¹, Yiwei Xie², Hongbin Zhang^1*

¹Fuzhou Medical College of Nanchang University, Fuzhou, China.
²Medical College, Jinggangshan University, Ji’an, China.

DOI: 10.4236/ajmb.2024.143009                

Abstract

Phospholipase A2 (PLA2) is the key enzyme to the venom from Deinagkistrodon acutus which is one of the highly venomous snakes in China. In addition to being a catalyst for the hydrolysis of phospholipases A2 from snake venom, it’s well known that it possesses a broad spectrum of pharmacological activities, such as myotoxicity, neurotoxicity, cardiotoxicity, and hemolytic, anticoagulant and antiplatelet activities. However, snakebites are not efficiently treated by conventional serum therapy. Acute wounds can still cause poisoning and death. In order to find effective inhibitors of Deinagkistrodon venom acid phospholipase A2 (dPLA2), we obtained 385 compounds in 9 Chinese herbs from the TCMSP. These compounds were further performed to virtual screen using in silico tools like ADMET analysis, molecular docking and molecular dynamics (MD) simulation. After Pharmacokinetics analysis, we found 7 candidate compounds. Besides, analysis of small molecule interactions with dPLA2 confirmed that the amino acid residues HIS47 and GLY29 are key targets. Because they bind not only to the natural substrate phosphatidylcholine and compounds known for having inhibitory functions, but also for combining with potential antidote molecules in Chinese herbal medicine. This study is the first to report experience with virtual screening for possible inhibitor of dPLA2, such as the interaction spatial structure, binding energy and binding interaction analysis, these experiences not only provide reference for further experimental research, but also have a guideline for the study of drug molecular mechanism of action.

Keywords

Chinese Herbal Medicine, Phospholipase A2 Inhibitor, Molecular Docking, Molecular Mechanism

1. Introduction

Venomous snakebite is a worldwide neglected public health emergency affecting predominantly poor residents in rural areas, especially in many tropical and subtropical countries [1] [2] [3] . According to estimation, between 1.2 million and 5.5 million snakebites happen annually all over the world. Moreover, approximately 400,000 victims have been left with permanent sequela, and 20,000 deaths occur each year [4] [5] .

Deinagkistrodon acutus, also named as Agkistrodon acutus, is one of the hypertoxic snakes in China [6] . The venom of Deinagkistrodon acutus are heterogeneous mixtures of proteins or peptides that act on prey synergistically [7] . Jin H and coworkers found that the most abundant constituent of the Deinagkistrodon acutus venom was phospholipases A2 (30.0%), snaclec (21.0%), antithrombin (17.8%), thrombin (8.1%) and metalloproteinases (4.2%) [8] . Deinagkistrodon acutus phospholipases A2 (dPLA2) that calcium ion dependent plays an essential role in the pathological process on the patient [7] [9] . Since this protease catalyzes the hydrolysis of the fatty acid ester at the sn-2 position of phospholipids releasing lysophospholipids and fatty acids [10] [11] , and exerts a wide variety of relevant toxic actions such as myotoxic, neurotoxic, cardiotoxic, anticoagulant, cytotoxic, edematogenic, and proinflammatory effects [12] [13] . Therefore, compounds capable of neutralizing dPLA2 may be effective in suppressing snake venom attacks.

Currently, the main treatment of snake envenomation is serum antivenom immunotherapy [14] . However, this approach experiences challenges with respect to poor antibody specificity, limited para-specificity, high incidence of adverse reactions, poor affordability and low availability to those who need them. Furthermore, it is difficult to neutralize the local tissue injury, since the rapid diffusion and reaction of the toxins present in Deinagkistrodon acutus venom reduce the antiserum efficiency [15] [16] [17] [18] . Undoubtedly, research and development for simple and efficient antivenom to complement this therapy is a promising study field, and a combination of functional and structural assays will be used to test candidate ligands against dPLA2.

Compared with classical antivenom serum, medicinal plants used in folk medicine, are a potential source of antivenom natural bioactive inhibitory compounds [19] [20] [21] [22] . Many studies indicate that some single-ingredient herbal medicines have potential therapeutic effects and value in the management of snakebites. Hedyotis diffusa Willd, Lobelia chinensis Lour, Houttuynia cordata Thunb, Sedum sarmentosum Bunge, Andrographis paniculata Nees, Pinellia ternata have demonstrated potential in the development of antivenom complementary therapies [23] [24] [25] . Pithayanukul et al. reported the ability of butanolic extracts from Eclipta prostata to neutralize the lethal effect, hemorrhagic effect and PLA2 activity of Calloselasma rhodostoma venom [26] . Phenolic compounds, isolated from the hexane extract of Curcuma longa roots, was shown to be a potent neutralizer of both the hemorrhagic effect induced by Bothrops jararaca and the lethal effect of Crotalus d. terrificus venoms [27] . Recently, the investigations of curative monomers in Chinese herbal medicines also made great progress [28] [29] [30] . Linsheng Zeng et al. screened 7 potential active ingredients from Cynanchum paniculatum against the venom of Bungarus multicinctus by Network Pharmacological method [31] . However, there are few reports about the treatment of Deinagkistrodon acutus venom by Chinese herbal monomer, dPLA2 have not been explored as a molecular target so far [32] .

Thus, in the present work, virtual screening, molecular docking, and molecular dynamics simulation were performed to test the potential inhibitors against dPLA2. Furthermore, we supplemented the molecular mechanism of interaction mode of complex, structural stability and flexibility, binding free energy et al.

2. Materials and Methods

2.1. The Acquisition of Software and the Use of Databases

The software used in this study includes AutodockTools-4.2.6, Autodock vina-1.1.2, Pymol-2.4.0 and GROMACS-2020.6, OpenBabel-3.1.1, Discovery studios Visualization-19.1.0 (DS), APBS-3.0.0 (Adaptive poisson-Boltzmann Solver) and gmx_mmpbsa.bsh; The database used in this study includes TCMSP-Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (tcmsp-e.com), Protein Data Bank(RCSB PDB: Homepage).

Among them, Autodock vina, used for molecular docking, is downloaded from AutoDock (scripps.edu); Pymol, which is installed using the whole file provided by UC Irvine Laboratory, was employed to process molecular docking image; Discovery Studio was applied for protein and small molecule preparation, virtual screening (ADMET), batch molecular docking and protein ligand complex interaction display and 2D mapping; TCMSP database was used to obtain small molecules of target Chinese herbal medicine and preliminary screening; OpenBabel is able to convert files, it is downloaded from https://openbabel.org/; GROMACS is available in Molecular Dynamics module (https://www.gromacs.org/); APBS and gmx_mmpbsa.bsh are used for calculating the binding free energy, they were downloaded from https://github.com/Electrostatics/apbs/releases and https://github.com/Jerkwin/gmxtool respectively.

2.2. Ligand and Receptor Preparation

2.2.1. Ligand Preparation

Based on TCMSP, 9 Chinese herbal medicines (Eclipta prostrata L, Hedyotis diffusa Willd, Lobelia chinensis Lour, Houttuynia cordata Thunb, Sedum sarmentosum Bunge, Andrographis paniculata Nees, Cynanchum paniculatum, Curcuma longa L, Pinellia ternata), which have been claimed to have a therapeutic effect on snakebite, were selected for this study. They were retrieved in TCMSP and then screened according to the criterion of drug-like properties (DL) ≥ 0.18 and oral bioavailability (OB) ≥ 30% [33] . Then those compounds which satisfied the standard, were downloaded as ligands and imported into OpenBabel for building files as “sd” or “sdf” format. Subsequently, these files were inputted into DS Minimize Ligands module for hydrogenation, adding charge and structural optimization under CHARMm force field.

2.2.2. Receptor Preparation

The important clinical manifestations of patients who were bitten by Deinagkistrodon acutus are tissue edema and bleeding, which are closely associated with acidic PLA2 in snake venom. So the inhibition of acidic PLA2, has important significance for relieving tissue edema and bleeding. Based on it search for the acidic PLA2 in the PDB database as a receptor for molecular docking, the results showed that there is only one acidic PLA2 of known structure (PDB ID: 1IJL).The crystal structure of dPLA2 (PDB ID: 1ijl) was accessed from Protein Data Bank (PDB) of RCSB (http://www.rcsb.org/) as receptor. Subsequently, the protein was modified using DS, including removal of water, adding hydrogen, and replenishing of lost residues. By means of CHARMm force field (OPLS3e), the protein structure was optimized to assure the structural accuracy of final receptor [10] .

2.3. Pharmacokinetics Analysis

ADME screening and toxicological prediction were used for virtual screening of the identified small molecules. The absorption, distribution, metabolism, excretion and toxicity (ADMET) properties were analyzed and assessed using DS Small Molecules-ADMET online server [34] . The evaluation criterion of ADMET was that Absorption_level ≤ 1, 2 ≤ Solubility_level ≤ 4, BBB_level ≤ 3. The prediction model we use was Developmental Toxicity Potential DTP, Rat Oral LD50 and Ames Mutagenicity. Finally the toxicity prediction results were derived.

2.4. Molecular Docking

2.4.1. Batch Molecular Docking

The compounds through pharmacokinetics screened was performed to batch molecular docking by the Dock Ligand modules. The binding sites were Y27, G29, G31, H47, D48, Y51, and D98, so the active center was X = 6.826, Y = 12.735, Z = 50.709, the grid was set as 15 × 15 × 15. Docking type was set as Libdock. Moreover, the docking score prediction of the molecules was performed, each molecule retained the highest score and arranged these scores in descending order.

2.4.2. Precision Molecular Docking

The molecules with the top score conformations were further redocked using Autodock vina with the same parameters as batch molecular docking to generate binding energy, hydrogen bonding and hydrophobic interaction results. The docking results were saved in pdb format respectively. Then we imported the results into PYMOL in pdbqt format for visualization processing and into DS in pdb format to obtain 2D structure diagram.

2.5. Molecular Dynamics (MD) Simulations

Molecular dynamics simulation (at 100 ns) was employed to explore the stability of ligand-protein combinations and analyze the molecular mechanisms of protein-target interactions. The top docking score of compound conformations from the molecular docking process were subjected to the GROMACS 2020 (Groningen Machine for Chemical Simulations Calculations) [35] and Amber99sb-ildn force field for Molecular dynamics simulation.

2.5.1. Preparation Topology

Initial configurations of ligands were produced with ACPYPE (http://bio2byte.be/acpype/) Server to prepare the topology structures, simultaneously, the Pdb2gmx was applied to generate the protein topology and structure files.

2.5.2. Preparation Reaction System

Firstly, the protein was centered on the cube box. Then the simple point charge (SPC) water model was solved. For the purpose of energy minimization of the system, we eliminated collisions at atomic positions through implement the fastest descent method. Due to the exposed residual charge of the protein given the system -12 net charge, so 12 Na⁺ counter ions were added to conform electroneutrality ahead of the actual simulation. Next, the band constraint of LINCS algorithm was applied to the system, the long range electrostatic interaction was calculated by particle grid Ewald (PME) algorithm. Afterwards, the NVT and NPT ensembles (constant number of atoms, volume, pressure, and temperature) were performed to the system for the MD simulation studies. In addition, at 1.4 nm, the cutoff values of Van der Waals interaction and short range Coulomb were set respectively in the system, periodic boundary conditions were applied to the whole three-dimensional spaces. Finally, 100 ns intervals of MD simulations were carried out for all the complexes with the time step 2fs, as well as the following data were investigated, including RMSD, RMSF and potential energy of protein-ligand complexes, radial distribution of ligand, complexes free energy profile and hydrogen bonds.

2.6. MD Data Processing

It was discovered that the protein or the ligand-protein Root Mean Square Deviation (RMSD) and Root Mean Square Fluctuation (RMSF) plots, as well as patterns of protein-ligand interaction, are markers of ligand stability.

2.6.1. RMSD (Root Mean Square Deviation)

The root mean square deviation reflects the degree of data deviation from the mean value. In the dynamic simulation of protein-ligand complexes, it represents the degree of fit or deviation between the conformation after simulated “nt” time and the initial conformation. The greater the RMSD value, the greater the difference between the conformation after “nt” time simulation and the original conformation, that is, the worse the conformation stability. On the contrary, the more stable it is. The calculation formula is below:

$\text{RMSD}=\sqrt{\frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\left({\delta }_i-{\delta }_0\right)}^2}}$

${\delta }_i$ represents the position of the conformation of the I-th frame, and ${\delta }_0$ represents the position of the initial conformation.

2.6.2. RMSF (Root Mean Square Fluctuation)

The root-mean-square fluctuation indicates the position change of an atom over a period of time. It also reflects the degree of flexibility of the atom. The calculation was derived as follows:

$\text{RMSF}=\sqrt{\frac{1}{T}{\displaystyle \underset{t_j=1}{\overset{T}{\sum }}{\left(x_i\left(t_j\right)-x_0\right)}^2}}$

$x_i\left(t_j\right)-x_0$ represents the position of an atom at “t” time subtraction the position at the initial time.

2.7. MM-PBSA Calculation (Molecular Mechanics-Poisson-Boltzmann Superficial Area)

The expected binding energy and the standard deviation were calculated for the complete 100 ns time interval using gmx_mmpbsa.bsh. The working principle is to count the difference between the binding free energy of two solvated molecules in bound and unbound states, or to compare the free energy of the same molecule in different solvated conformations. That is, decomposing the binding free energy into molecular mechanical terms and solvation energy, they were computed separately. The binding free energy was calculated using the following formula:

$\Delta G_{\text{bind}，\text{solv}}^0=\Delta G_{\text{bind},\text{vacuum}}^0+\Delta G_{\text{solv}，\text{complex}}^0-\left(\Delta G_{\text{solv}，\text{ligand}}^0+\Delta G_{\text{solv}，\text{receptor}}^0\right)$

$\Delta G_{\text{bind}，\text{solv}}^0$ denotes the binding free energy under solvent; $\Delta G_{\text{bind},\text{vacuum}}^0$ expresses the binding energy of the interaction between molecules and proteins with the addition of solvent; $\Delta G_{\text{solv}，\text{complex}}^0$ represents the energy of the complexes in the solvent; $\Delta G_{\text{solv}，\text{complex}}^0$ represents the energy of the receptor in the solvent.

3. Results and Discussion

3.1. Screening Preliminary

3.1.1. TCMSP Analysis

Approximately, 385 chemical components from 9 kinds of herbs (Eclipta prostrata L, Hedyotis diffusa Willd, Lobelia chinensis Lour, Houttuynia cordata Thunb, Sedum sarmentosum Bunge, Andrographis paniculata Nees, Cynanchum paniculatum, Curcuma longa L, Pinellia ternata), with potential antivenomous activity,were carried out to virtual screening by TCMSP. Under the above conditions, 89 compounds were selected for further assessment.

3.1.2. Pharmacokinetics Analysis

To exclude these molecules that may be harmful to human body, it is critical to screen the candidates to fit for drug-like properties. After ADME analysis,the distribution of 89 molecules are as shown in Figure 1(a). After ADME screening, 7 molecules were identified with a greater potential as the inhibitors against dPLA2 as shown in Figure 1(b). They are presented in Table 1 as, 3-[(2S,3R)-2-(4-hydroxy-3-methoxy-phenyl)-3-(hydroxymethyl)-7-methoxy-2, 3-dihydrobenzofuran-5-yl] propyl acetate; Ruvoside-qt; 14-deoxy-12-methoxyandrographolide; 14-deoxy-11-oxo-andrographolide; (3S,6S)-3-(benzyl)-6-(4-hydroxybenzyl) piperazine-2,5-quinone (MOL006957); tomentogenin (MOL005622) and Isoramanone. Their chemical structures are shown in Figure 1(c).

The toxicity of seven molecules was predicted using three models of DS: Mutagenicity, Developmental_Toxicity_Potential and Rat_Oral_LD50. A prediction score less than 0 indicates relative nontoxicity, a score greater than 0 indicates toxicity, and a larger value indicates greater toxicity. Toxicity prediction revealed that of the seven compounds, only 14-deoxy-11-oxo-andrographolide had no potential developmental toxicity. From the perspective of oral half-lethal dose in rats, the highest value was found on tomentogenin (8.12 g/kg), with the potential developmental toxicity of 0.47699, indicating that it was less toxic. From the perspective of mutagenicity, none of the 7 molecules were teratogenic, and the lowest teratogenic was Ruvoside-qt (−37.5246). Moreover the oral lethal dose of Ruvoside-qt was 6.38, which was relatively large, suggesting that it was less toxic.

3.2. Molecular Docking

We learned from the PDB database that the amino acid residues (Y27, G29, G31 and D48) of the calcium ion binding ring are the active sites of dPLA2. It indicated that dPLA2 is a calcium-dependent protease. However, the enzyme activity of dPLA2 is determined by three main factors, including integrity of the active site (His48, Asp49, Tyr52, Asp99), coordination of Ca²⁺ cofactors (Tyr28, Gly30, Gly32 and Asp49), and Enzyme adsorption on outer surface of membrane lipid [36] . As a result, the regions composed of Y27, G29, G31, H47, D48, and D89 were set as binding sites. The center was X = 6.826, Y = 12.735, Z = 50.709, and the grid was 15 × 15 × 15. The above parameters were used for the following docking experiments.

3.2.1. Batch Molecular Docking Results

Libdock batch molecular docking was performed to predict the binding possibility of the seven ligands to the active site of dPLA2. As shown in Table 1, the top docking score is Ruvoside-qt with 140.465, which indicated that Ruvoside-qt was the most likely to combine and function with dPLA2. The other compounds were sorted by score for 3-[(2S,3R)-2-(4-hydroxy-3-methoxy-phenyl)-3-(hydroxylmethyl)-7-methoxy-2,3-dihydrobenzofuran-5-yl] propyl acetate; tomentogenin; 14-deoxy-11-oxo-andrographolide; 14-deoxy-12-methoxyandrographolide et al.

Figure 1. Virtual screening results in (a) and (b) show the details of ADMET screening, (a) shows the distribution of all the small molecules after ADMET analysis. (b) shows the distribution of the small molecules with Absorption_level ≤ 1, 2 ≤ Solubility_level ≤ 4, BBB_level ≤ 3, meanwhile EXT_CYP2D6_Prediction (cytochrome P450 2D6), EXT_PPB_Prediction and EXT_Hepatotoxic_Prediction were set as false. (c) shows the chemical structural formula of the seven small molecules screened. The top three are MOL012225, MOL004350, MOL008213 from left to order, and the bottom four are MOL008203, MOL005622, MOL006957, MOL003851 from left to right.

Table 1. The virtual screening results.

Note: DTP is the potential for developing toxicity, the smaller the value indicates the smaller the potential for developing toxicity, positive value has the potential for developing toxicity (Toxic), negative value indicates no potential for developing toxicity (non-Toxic); Ames score is mutagenicity, the smaller the value indicates the lower mutagenicity, positive value indicates mutation (Mutagenicity), and negative value indicates no mutation (non-Mutagenicity).

3.2.2. Positive Control-Phosphatidylcholine

Phosphatidylcholine (1,2-Diacyl-SN-glycero-3-phosphocholine), the natural substrate of dPLA2, is the most abundant phosphatidylcholine component in cell membranes. The hydrolase dPLA2 breakdowns membrane phospholipids, which generates inflammatory mediators, such as arachidonic acid, prostaglandin and platelet activating factor, causing inflammation in humans. Our use of phosphatidylcholine-dPLA2 complex as a positive control to screen for competitive inhibitors is highly convincing. It can also verify the accuracy of the active site.

The experiments in vitro incubation found that Gallic acid and pyrogallol shown obvious inhibitory effects on dPLA2 in a concentration-dependent manner [37] [38] . For this reason, we employed Gallic acid and pyrogallol to interact with dPLA2 respectively, to verify the feasibility of the active site, laying a foundation for the accuracy of molecular docking and the feasibility of screening.

3.2.3. Precision Molecular Docking Results

Autodock vina was considered as an important computational tool in search of the interaction mode of the complex. It can more precisely predict the binding interactions between the compounds and the active site of the target protein. In this way, AutoDock Vina generated a docked complex for each ligand with various conformation and affinity scores. The best protein-ligand complex docked pose were analyzed and graphically visualized by DS.

As shown in Table 2, molecular docking revealed that phosphatidylcholine

Table 2. Comparison of the molecular docking interactions.

Note: Table 2 shows the docking results of the seven small molecules with dPLA 2 molecules (top 7), and the docking results of dPLA 2 with phosphatidylcholine, gallic acid and catechol (last 3). In Table, MOL006957 to MOL008213 are the potential active components obtained in virtual screening according to the size of docking binding energy; MOL008290 to MOL000106 are the small molecules that inhibit PLA2 and used as a positive control. The abbreviated VDW applied in the table represents the van der Waals.

combine well with the protein in the position of TRP30, GLY31, HIS47 and LYS60 through hydrogen bonds as shown in Figure 2(b), hydrophobic interactions with CYS44, HIS47 and PRO120, anion interactions with TRP30, charge attraction interactions with LYS60, and C-H bonds with GLY32. These binding

Figure 2. The molecular docking results in Fig. On the left is the molecular binding site and hydrogen bond interaction, where green indicates the ligand small molecule, blue indicates the receptor amino acid, and gray indicates the receptor surface structure diagram, showing the situation around the ligand binding site. The 2D structure diagram of the molecules interacting with the surrounding amino acids is shown on the right. (a) shows the interaction between Gallic acid and dPLA2, (b) shows the interaction between 1,2-diacyl-sn-glycero-3-phosphocholine and dPLA2, (c) shows the interaction between tomentogenin and dPLA2, and (d) represents the interaction between (3S, 6S)-3-(benzyl)-6-(4-hydroxybenzyl) piperazine-2,5-quinone and dPLA2.

sites are similar to the dPLA2 active site. Meanwhile comparing with phosphatidylcholine, gallic acid and pyrogallol have similar dPLA2 binding sites, both gallic acid and pyrogallol form hydrogen bonds with GLY29 and HIS47 as shown in Figure 2(a). Besides, all three form hydrogen bonds and hydrophobic interactions with HIS47. It then suggested that HIS47 may perform an important part in the catalytic activity of dPLA2.

In terms of binding energy, pyrogallol and phosphatidylcholine have the same binding energy of −5.7 kJ/mol, gallic acid −6.3 kJ/mol and combine with Ca²⁺ by metallic bond. Thus, according to the interaction of three controls with dPLA2, the catalytic active center of the receptor is a calcium-ion binding ring composed of TYR27, GLY29, GLY31, HIS47 and ASP48. Moreover HIS47 and GLY29 may be key targets in the catalytic activity.

After docking, 7 molecules are able to combine with the target protein of dPLA2. The binding energies are all less than −5.7 kJ/mol. The lowest binding energy is MOL006957 (−9.0 kJ/mol). The second least binding energy for Isoramanone being −8.1 kJ/mol. As it is known that the lower the energy, the stronger the bond, they were the ligands with the most binding affinity to the receptor. In terms of interaction, MOL006957 forms hydrogen bonds with GLY29 and GLY31 as shown in Figure 2(d), hydrophobic bonds with TRP30, metallic bonds with Ca²⁺, cations interaction with LYS60, and CYS 44 forms the Pi-sulfur bond and the C-H bond with HIS47, which is most similar to phosphatidylcholine in general. Isoramanone forms hydrogen bond interaction with HIS47 and ASP48 as shown in Figure 2(c), hydrophobic interaction with LEU2, TRP30 and TYR51, and metallic bond with Ca²⁺. Consequently, MOL006957 is most likely to be the ligand of dPLA2, which was further verified by molecular dynamics simulations.

3.3. Molecular Dynamics (MD) Simulations

3.3.1. RMSD and RMSF Analysis

Figure 3(a) shows RMSD of the backbone of the seven protein-ligand complexes which changed with the time over 100 ns. A conclusion can be drawn from the picture that after simulation, the conformations of these seven complexes tended to converge. Comparing with the initial conformation, the deviation of all the complexes were less than 0.4 nm. In particular, the RMSD of MOL005622-dPLA2 and MOL006957-dPLA2 complexes were lower than other compounds. The result indicated that the two compounds were relatively stable over 100 ns. Then, in comparison to dPLA2, RMSD of MOL005622-dPLA2 and MOL006957-dPLA2 complexes was lower overall and it suggested that the structural flexibility of dPLA2 were reduced after binding ligands, as shown in Figure 3(b).

RMSF reflected the vibration of each complex compared with the initial conformation, as shown in Figure 3(c). In comparison to other compounds, there was the smallest molecular vibration in MOL006957-dPLA2 complex, which reflects that it is more stable than other compounds at the molecular level. Simultaneously comparing with dPLA2, MOL006957-dPLA2 molecular vibration was

Figure 3. Represents the results of RMSD analysis and RMSF analysis after molecular dynamics simulation. (a) indicates the seven complexes curve of RMSD value as a function of simulation time. (b) represents the comparative plots of the RMSD fluctuations between MOL006957-dPLA 2 and dPLA2. (c) represents the RMSF values of each atom of the seven complexes. (d) represents the comparative plots of the RMSF values between MOL006957-dPLA 2 and dPLA2.

significantly weakened, as shown in Figure 3(d).

3.3.2. Free Energy Topography (FEL) Analysis and Hydrogen Bond Analysis

FEL is able to evaluate comprehensively the free energy of different conformations in simulation process through the following two parameters: root mean square error (RMSD) and gyrate. As shown in Figure 4(a), the lower the energy (the darker the blue), the more stable the conformation. Through the parameter files, it was found that the lowest energy conformation was frame 3670 and frame 3680, and compared the lowest energy conformation with the initial conformation, as shown in the Figure 4(b). Furthermore, hydrogen bonds play an important role in ligand-receptor interactions, and the number of hydrogen bonds is also an important index to judge the binding stability of protein-ligand complexes. As shown in the Figure 4(c), the overall number of hydrogen bonds increased with the extension of simulation time, and finally remained at 3 to 4.

Figure 4. Free energy topography analysis and hydrogen bond analysis results of MOL006957-dPLA 2. (a) shows the free energy morphology diagram, the other left side is the 2D diagram, the right side is the 3D structure diagram, the deeper the blue indicates that the smaller the free energy of the structure conformation and the more stable the conformation is. (b) shows the comparison of the initial conformation of the complex and frame 424, with the initial conformation in green and the other 424 frames, in which frame 424 is the lowest conformation obtained by the analysis of the free energy topography. (c) shows the results of the hydrogen bond analysis showing the number of hydrogen bonds during simulation with the simulation time.

3.4. Binding Free Energy Calculation

The Poisson-Boltzmann superficial area approach from molecular mechanics (MM/PBSA) was used to determine the binding free energy of protein-ligand complexes. In actual fact, the amount of spuriously positive outcomes in molecular docking study is decreased owing to this accurate quantum mechanical computation.

The binding free energy reflects the stability of the ligand-protein complexes in the solvent, and the influence of the solvent on the process of ligand combining with protein. The trajectory of two complexes (MOL006957-dPLA2 and MOL005622-dPLA2) with the best molecular dynamics simulation results, were performed to compute the binding free energy. Then the stability of the complexes in the solvent was checked. The results are shown in Table 3. The average binding energy of MOL006957-dPLA2 and MOL005622-dPLA2 in solvent respectively are −99.334 kJ/mol and −73.242 kJ/mol, which expresses that the solvent has little effect on the stability of the complex, and the ligand binding is quite stable. However, the average binding energy of others are relatively low, revealed that the solvent has a great influence on the stability of the complex, and the ligand binding is unstable. Thus, we provide stronger evidence for the potential inhibitory effect of MOL006957 on the dPLA2 venom of Agkistrodon acutus.

4. Conclusions

Pinellia ternata is one of the most commonly used traditional Chinese medicinal plants. It has been demonstrated to be highly effective for treating cough, vomiting, infection, and inflammatory diseases. Modern pharmacological investigations have demonstrated its multiple activities, such as antitussive, expectorant, antiemetic, antitumor, antibacterial, and sedative-hypnotic activities [39] [40] .

The compound MOL006957 screened in this study is one of the active potential drug like molecule of Pinellia [41] . It has been reported that the active

Table 3. Comparison of the binding free-energy results of compounds.

Note: dH indicates the enthalpy difference, to present the enthalpy difference in the presence of the protein-ligand complex and the protein and ligand alone; MM: van der Waals and electrostatic energy (COU + VDW); PB: polar solvation; SA: nonpolar solvation; COU: electrostatic energy; VDW: van der Waals; the dG indicates the free energy difference, to assess the binding forces between the ligand receptors, the more negative the dG, Me means stronger binding.

ingredient (MOL006957) of Pinellia in Chinese prescription [42] has potential medicinal value. In the course of this experiment, it was found that MOL006957 showed outstanding advantages compared with other compounds in the screening. The pharmacokinetic analysis, molecular docking or molecular dynamics simulation, indicated that the ligand binding is quite stable, and can play an important role in the treatment of snake venom and neutralization of dPLA2. However, whether this molecule can be used to inhibit snake venom inside the human body remains unclear, further in vivo experiments are needed to validate the potential clinical utility of the molecules identified.

NOTES

*Corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	World Health Organization (2007) Rabies and Envenomings. A Neglected Public Health Issue. Geneva. https://www.who.int/publications/i/item/9789241563482
[2]	Gutiérrez, M.J., Warrell, A.D. and Williams, J.D. (2017) The Need for Full Integration of Snakebite Envenoming within a Global Strategy to Combat the Neglected Tropical Diseases: The Way Forward. PLOS Neglected Tropical Diseases, 7, e2162.
[3]	Gutiérrez, M.J., León, G. and Burnouf, T. (2011) Antivenoms for the Treatment of Snakebite Envenomings: The Road Ahead. Biologicals, 3, 129-142. https://doi.org/10.1016/j.biologicals.2011.02.005
[4]	Kasturiratne, A., Wickremasinghe, R.A., Silva, D.N., et al. (2008) The Global Burden of Snakebite: A Literature Analysis and Modelling Based on Regional Estimates of Envenoming and Deaths. PLOS Medicine, 11, e218. https://doi.org/10.1371/journal.pmed.0050218
[5]	María, J.G., Thierry, B., Harrison, R.A., et al. (2015) A Call for Incorporating Social Research in the Global Struggle against Snakebite. PLOS Neglected Tropical Diseases, 9, e0003960. https://doi.org/10.1371/journal.pntd.0003960
[6]	Nie, X.K., He, Q.Y. and Zhou, B. (2021) Exploring the Five-Paced Viper (Deinagkistrodon acutus) Venom Proteome by Integrating a Combinatorial Peptide Ligand Library Approach with Shotgun LC-MS/MS. The Journal of Venomous Animals and Toxins Including Tropical Diseases, 27, e20200196. https://doi.org/10.1590/1678-9199-jvatitd-2020-0196
[7]	Tasoulis, T. and Isbister, G.K. (2023) A Current Perspective on Snake Venom Composition and Constituent Protein Families. Archives of Toxicology, 97, 133-153. https://doi.org/10.1007/s00204-022-03420-0
[8]	Jin, H., Minrui, Z., Chu, X., Liang, J.Q. and Huang, F. (2022) Analysis of the Composition of Deinagkistrodon acutus Snake Venom Based on Proteomics, and Its Antithrombotic Activity and Toxicity Studies. Molecules, 27, Article 2229.
[9]	Kong, Y., Sun, Q., Zhao, Q. and Zhang, Y. (2018) Purification and Characterization of a Novel Antiplatelet Peptide from Deinagkistrodon acutus Venom. Toxins, 10, Article 332. https://doi.org/10.3390/toxins10080332
[10]	Gu, L.C., Zhang, H.L., Song, S.Y., Zhou, Y. and Lin, Z. (2002) Structure of an Acidic Phospholipase A2 from the Venom of Deinagkistrodon acutus in a New Crystal Form. Acta Biochimica et Biophysicas Inica, 34, 266-272. https://doi.org/10.2210/pdb1ijl/pdb
[11]	Zou, Z., Zeng, F., Zhang, L., Niu, L., Teng, M. and Li, X. (2012) Purification, Crystallization and Preliminary X-Ray Diffraction Analysis of an Acidic Phospholipase A2 with Vasoconstrictor Activity from Agkistrodon halys pallas Venom. Acta Crystallographica, 68, 1329-1332. https://doi.org/10.1107/S1744309112038523
[12]	Chinnasamy, S., Selvaraj, G., Selvaraj, C., et al. (2020) Combining in Silico and in Vitro Approaches to Identification of Potent Inhibitor against Phospholipase A2 (PLA2). International Journal of Biological Macromolecules, 144, 53-66. https://doi.org/10.1016/j.ijbiomac.2019.12.091
[13]	Costa, S.K.P., Camargo, E.A. and Antunes, E. (2017) Inflammatory Action of Secretory Phospholipases A2 from Snake Venoms. In: Cruz, L. and Luo, S., Eds., Toxins and Drug Discovery, Springer, Dordrecht, 35-52. https://doi.org/10.1007/978-94-007-6452-1_10
[14]	Timothy, O., Sarah, K., Abraham, O., et al. (2020) Antivenin Plants Used for Treatment of Snakebites in Uganda: Ethnobotanical Reports and Pharmacological Evidences. Tropical Medicine and Health, 48, Article No. 6.
[15]	Williams, J.D., Gutiérrez, J., Calvete, J.J., et al. (2011) Ending the Drought: New Strategies for Improving the Flow of Affordable, Effective Antivenoms in Asia and Africa. Journal of Proteomics, 74, 1735-1767. https://doi.org/10.1016/j.jprot.2011.05.027
[16]	Maria, J.G., et al. (2017) Snakebite Envenoming. Nature Reviews Disease Primers, 3, Article No. 17063.
[17]	Stuart, A., Julien, S., Nessrin, A., et al. (2018) The Paraspecific Neutralisation of Snake Venom Induced Coagulopathy by Antivenoms. Communications Biology, 1, Article No. 34. https://doi.org/10.1038/s42003-018-0039-1
[18]	Harrison, R.A., Casewell, N.R., Ainsworth, S.A. and Lalloo, D.G. (2019) The Time Is Now: A Call for Action to Translate Recent Momentum on Tackling Tropical Snakebite into Sustained Benefit for Victims. Transactions of the Royal Society of Tropical Medicine and Hygiene, 113, 835-838. https://doi.org/10.1093/trstmh/try134
[19]	Suveena, S., Saraswathy, V., et al. (2022) In Silico Screening of the Phytochemicals Present in Clitoria ternatea L. as the Inhibitors of Snake Venom Phospholipase A2 (PLA2). Journal of Biomolecular Structure Dynamics, 41, 7874-7883.
[20]	Upasana, P., Alexandrino, P.F. and Mukherjee, A.K. (2022) Pharmacological Re-Assessment of Traditional Medicinal Plants-Derived Inhibitors as Antidotes against Snakebite Envenoming: A Critical Review. Journal of Ethnopharmacology, 292, Article ID: 115208. https://doi.org/10.1016/j.jep.2022.115208
[21]	Silva, D.P.D., Ferreira, S.D.S., Torres-Rêgo, M., et al. (2022) Antiophidic Potential of Chlorogenic Acid and Rosmarinic Acid against Bothrops leucurus Snake Venom. Biomedicine & Pharmacotherapy, 148, Article ID: 112766. https://doi.org/10.1016/j.biopha.2022.112766
[22]	Carvalho, B.M.A., Santos, J.D.L., et al. (2013) Snake Venom PLA2s Inhibitors Isolated from Brazilian Plants: Synthetic and Natural Molecules. BioMed Research International, 2013, Article ID: 153045. https://doi.org/10.1155/2013/153045
[23]	Chen, J.J., Shi, D.J. and Li, K.H. (2005) Effect of Mohagan on the Inflammation and Bleeding Induced by Venom of Four Species of Pit Viper. Ophistoma, 2, 65-68.
[24]	Liu, Z.K., Li, Q.P. and Zhou, W.Z. (2013) Effect of Clearing Heat, Cooling Blood and Detoxification on Coagulation Function of Patients Bitten by Five-Step Snake. Journal of Nanjing University of Chinese Medicine, No. 1, 12-15.
[25]	Dey, A. and De, J.N. (2012) Phytopharmacology of Antiophidian Botanicals: A Review. International Journal of Pharmacology, 8, 62-79. https://doi.org/10.3923/ijp.2012.62.79
[26]	Pithayanukul, P., Laovachirasuwan, S., Bavovada, R., Pakmanee, N. and Suttisri, R. (2004) Antivenom Potential of Butanolic Extract of Eclipta prostrata against Malayan Pit Viper Venom. Journal of Ethnopharmacology, 90, 347-352. https://doi.org/10.1016/j.jep.2003.10.014
[27]	Soares, A.M., Ticli, F.K., et al. (2005) Medicinal Plants with Inhibitory Properties against Snake Venoms. Current Medicinal Chemistry, 12, 2625-2641. https://doi.org/10.2174/092986705774370655
[28]	Ma, Y.J. (2021) Research Progress on the Biology and Anti-Tumor of Artemisinins. Modern Distance Education of Chinese Medicine, 24, 192-196.
[29]	Ju, T.T., Liu, J.F. and Li, X.D. (2021) Research Progress on Extraction Method and Application of Salicylic Acid. Tianjin Chemical Industry, No. 4, 7-9.
[30]	Fan, H.H., Wang, L.Q., Liu, W.L., et al. (2020) Repurposing of Clinically Approved Drugs for Treatment of Coronavirus Disease 2019 in a 2019-Novel Coronavirus-Related Coronavirus Model. Chinese Medical Journal, 133, 1051-1056. https://doi.org/10.1097/CM9.0000000000000797
[31]	Zeng, L., Hou, J.J., Ge, C.H., et al. (2022) Network Pharmacological Study on the Mechanism of Cynanchum Paniculatum (Xuchangqing) in the Treatment of Bungarus Multicinctus Bites. BioMed Research International, 2022, Article ID: 3887072. https://doi.org/10.1155/2022/3887072
[32]	Quach, N.D., Arnold, R.D. and Cummings, B.S. (2014) Secretory Phospholipase A2 Enzymes as Pharmacological Targets for Treatment of Disease. Biochemical Pharmacology, 90, 338-348. https://doi.org/10.1016/j.bcp.2014.05.022
[33]	Ru, J., Li, P., Wang, J., et al. (2014) TCMSP: A Database of Systems Pharmacology for Drug Discovery from Herbal Medicines. Journal of Cheminformatics, 6, Article No. 13. https://doi.org/10.1186/1758-2946-6-13
[34]	Daina, A., Michielin, O. and Zoete, V. (2017) SwissADME: A Free Web Tool to Evaluate Pharmacokinetics, Drug-Likeness and Medicinal Chemistry Friendliness of Small Molecules. Scientific Reports, 7, Article No. 42717. https://doi.org/10.1038/srep42717
[35]	Abraham, M.J., Murtola, T., Schulz, R., et al. (2015) GROMACS: High Performance Molecular Simulations through Multi-Level Parallelism from Laptops Tosupercomputers. SoftwareX, 1, 19-25. https://doi.org/10.1016/j.softx.2015.06.001
[36]	Cedro, A.C.R., Menaldo, L.D., Costa, R.T., et al. (2018) Cytotoxic and Inflammatory Potential of a Phospholipase A2 from Bothrops jararaca Snake Venom. Journal of Venomous Animals and Toxins Including Tropical Diseases, 24, Article No. 33. https://doi.org/10.1186/s40409-018-0170-y
[37]	Hikari, M.T., Airam, R., de Ferreira, L.L.M., et al. (2022) Gallic Acid as a Non-Selective Inhibitor of α/β-Hydrolase Fold Enzymes Involved in the Inflammatory Process: The Two Sides of the Same Coin. Pharmaceutics, 14, 368.
[38]	Costa, T.R., Francisco, A.F., Cardoso, F.F., et al. (2021) Gallic Acid Anti-Myotoxic Activity and Mechanism of Action, a Snake Venom Phospholipase A2 Toxin Inhibitor, Isolated from the Medicinal Plant Anacardium humile. International Journal of Biological Macromolecules, 185, 494-512. https://doi.org/10.1016/j.ijbiomac.2021.06.163
[39]	Bai, J., Qi, J., Yang, L., et al. (2022) A Comprehensive Review on Ethnopharmacological, Phytochemical, Pharmacological and Toxicological Evaluation, and Quality Control of Pinellia ternata (Thunb.) Breit. Journal of Ethnopharmacology, 298, Article ID: 115650. https://doi.org/10.1016/j.jep.2022.115650
[40]	Mao, R. and He, Z. (2020) Pinellia ternata (Thunb.) Breit: A Review of Its Germplasm Resources, Genetic Diversity and Active Components. Journal of Ethnopharmacology, 263, Article ID: 113252. https://doi.org/10.1016/j.jep.2020.113252
[41]	Wang, L.C., Wang, D.G., Lin, T.D., et al. (2022) Medication Rule and Potential Mechanism of the Core Drug Pairs in Treating Reflux Esophagitis by LIN Tiandong, the National Master of TCM. Journal of Hunan University of Chinese Medicine, 42, 1493-1501.
[42]	Luo, W., Ding, R., Guo, X., et al. (2022) Clinical Data Mining Reveals Gancao-Banxia as a Potential Herbal Pair against Moderate COVID-19 by Dual Binding to IL-6/STAT3. Computers in Biology and Medicine, 145, Article ID: 105457. https://doi.org/10.1016/j.compbiomed.2022.105457