Pharmacoscreening, molecular dynamics, and quantum mechanics of inermin from Panax ginseng: a crucial molecule inhibiting exosomal protein target associated with coronary artery disease progression

Collection and optimization of natural compounds

The bioactive compounds of Panax ginseng were identified and its structures were obtained from the Traditional Chinese Medicine Systems Pharmacology (TCMSP) database (https://tcmsp-e.com/tcmsp.php). Then, the retrieved compound structures were optimized with a OPLS3force-field using the LigPrep tool in the Schrödinger suite (Maestro 11.2). Likewise, the known CAD drugs (amlodipine, furosemide, digoxin, clopidogrel, and atorvastatin) were downloaded as structure data files (SDF) format from PubChem (https://pubchem.ncbi.nlm.nih.gov/) for comparative analysis. Then the optimized natural compounds were screened with QikProp (Maestro, Schrödinger suite) (Ntie-Kang et al., 2013) to predict their absorption, distribution, metabolism and excretion (ADME) property based on molecular descriptors such as molecular weight, donor hydrogen bond, acceptor hydrogen bond, water partition coefficient, Caco-2 cell permeability, central nervous system (CNS) activity, blood/brain partition coefficient, cell permeability, oral absorption, serum albumin binding, Lipinski rule of five, and Jorgensen’s rule (Ntie-Kang et al., 2013; Rajagopal et al., 2020).

CAD genes and exosome

A gene expression array dataset from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) was selected with the following criteria that includes, (1) dataset containing expression profile of more than three samples in a group, (2) dataset containing CAD participants with the appropriate control executed in the CAD sample from the cardiac associated tissue, (3) the dataset providing human gene expression profile that uses microarray or RNAseq technique was considered, and (4) experiment carried out in animal or tissue culture were eliminated/excluded. To the selected dataset, limma R-program was implemented to identify the differentially expressed genes (DEGs) in CAD with log2fold-change (FC) and p-value < 0.05. Among DEGs, only the over expressed genes were selected for further assessment.

Intersection of DEGs with exosomal proteins cargo

To gather the information on exosomal cargo proteins, the EXOCARTA database was used (http://www.exocarta.org/) which contains the exosomal cargo of DNA, mRNA, miRNA, lipids, and proteins from multiple species. All human exosomal cargo proteins were downloaded and duplicates were removed to obtain an unique protein list. These exosomal unique lists were mapped to the over expressed DEGs of CAD to discover a common protein list (SET-A). The proteins in SET-A were demonstrated to have qualities (1) over expressed in CAD and (2) carried via exosome. Then the SET-A was used to build an inter-protein network with the STRING database (https://string-db.org/). Thereby, the constructed network contains both SET-A and its interacting proteins, resembling that SET-A on exosomal delivery into a target cell can readily interact and activate the cellular proteins to perform functions.

Enrichment analysis

The function of these network proteins was assessed by gene ontology (GO) and KEGG pathway using ShinyG0 0.77 (http://bioinformatics.sdstate.edu/go/). The GO analysis was carried out to identify the biological process, molecular function, and cellular component of the network proteins with a false discovery rate (FDR) less than 0.05. Similarly, the molecular pathways analysis were explored using KEGG pathway module (FDR < 0.05).

Molecular docking

From the protein interaction network, the SET-A proteins with a high degree of interaction were selected as target and their 3D protein structure was retrieved from the Protein Data Bank (PDB) database. Using a protein preparation module (Schrödinger suite), the target protein was optimized, which removes water molecules, heteroatoms, ions, and unwanted protein chains. Similarly, the LigPrep module (Schrödinger suite) was used for compound preparation such as adding hydrogen bonds. The binding pockets of protein were identified by P2Rank tool (https://prankweb.cz/) for grid generation around the pockets (Umashankar et al., 2021). Simultaneously, the OPLS3 force-field was applied to protein and small molecules (natural compounds + CAD drugs). A Glide suite (Schrödinger suite) was used to perform docking, to find the molecular interaction between the protein with natural compounds and also with the known CAD drugs. Finally, the small molecules with least binding energy was considered for molecular dynamic simulation.

Molecular dynamics simulation

After gaining the understanding of the ligand–protein complex’s through docking, the molecular dynamics (MD) simulation was performed for the target protein with the compound from Panax ginseng and CAD drug having highest binding affinity. For MD simulations, the Desmond tool included in the Schrödinger Drug Design Suite was utilized. To investigate the stability of the protein-ligands complex was carried out for 100 nanosecond (ns). Each simulation was carried out in three phases involving system construction, minimization, and the simulation run. Following the determination of the docked ligand–protein complex, the system was modeled using the SPC solvent system within the orthorhombic boundary and ions neutralized then subjected to energy minimization until it reached a gradient threshold of 25 kcal/mol at a pressure of 1 bar and a temperature of 300 K using the NPT ensemble class. The OPLS2005 force field optimization was used to minimize the system by carrying 2,000 iterations (Convergence scale of 1 kcal/mol/Å). A Nose–Hoover Chain thermostat was used to maintain the temperature, while a Martyna-Tobias-Klein barostat was utilized to sustain pressure with the equilibrium was set 200 ps. Finally, root-mean-square-deviation (RMSD), RMSF root-mean-square-fluctuation (RMSF) trajectories of the protein and ligand, protein–ligand interactions, and contacts with a variety of amino acids were used to determine the stability of the complex during the MD simulation.

Quantum computation calculation

The DFT calculation was preceded by using the Jaguar suite in Schrödinger (Yele et al., 2021). The predictions of chemical compound properties were calculated with the help of the integrated function of the B3LYP technique along with 6-311G** ⁺⁺ basis set (Yele et al., 2021). The geometrical parameters have been calculated theoretically including bond angle, bond length, and dihedral angles (Deghady et al., 2021). Similarly, highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) values were utilized to calculate other parameters (such as chemical potential, optical softness, chemical hardness, net electrophilicity, electrophilicity index, nucleophilicity index, and other quantum chemical parameters) by using Formulas 1–8. Likewise, the B3LYP/6-311++G (d,p) method in DFT was employed to identify the molecular electrostatic potential (ESP) and Mulliken charge analysis of the lead compound (Deghady et al., 2021).

(1)${\displaystyle \Delta \mathrm{E}={\mathrm{E}}_{\text{LUMO}}-{\mathrm{E}}_{\text{HOMO}}}$ (2)${\displaystyle \mathrm{A}=-{\mathrm{E}}_{\text{LUMO}}}$ (3)${\displaystyle \mathrm{I}=-{\mathrm{E}}_{\text{HOMO}}}$ (4)${\displaystyle \chi =\frac{\mathrm{I}+\mathrm{A}}{2}}$ (5)${\displaystyle \eta =\mathrm{I}-\mathrm{A}}$ (6)${\displaystyle \omega =\frac{{\mu }^2}{2\eta }}$ (7)${\displaystyle \mu =-\frac{\left(\mathrm{I}+\mathrm{A}\right)}{2}}$ (8)${\displaystyle \sigma =\frac{1}{\eta }}$

Compound collection and ADMET prediction

A total of 190 compounds of Panax ginseng reported in the TCMSP database (Accessed on: 14 April 2023) were collected and optimized with LigPrep module. Then QikProp was used to screen 190 compounds that predict the ADME properties. On ADME screening, 27 compounds were noticed to fulfill most drug-likeness properties (File S1) and selected for molecular docking.

Differential expression analysis

Meanwhile, our search in NCBI, GEO datasets (Accessed on: 30 April 2023) allowed us to identify the suitable dataset GSE120774 based on the selection criteria. The GSE120774 dataset consists of 19 human epicardial adipose tissue (EAT) from ten controls and nine CAD that experimented in microarray platform GPL6244 ((HuGene-1_0-st) Affymetrix Human Gene 1.0 ST Array (transcript (gene) version)). By employing limma programming, a total of 3,499 DEGs (1,408 up and 2,091 down-regulated) with log2FC and p < 0.05 were identified. However, only the over-expressed 1,408 genes in CAD were selected for the subsequent assessment.

Protein–protein interaction and enrichment analysis

Using the data retrieved from EXOCARTA, 6,394 unique exosomal proteins were identified and mapped with the 1,408 over-expressed CAD genes. A total of 489 proteins identified to be common were termed as SET-A, as mentioned in the methodology section. Then, the protein network that was constructed with the 489 (SET-A) proteins was extended to have 1,487 (489+998) proteins to form a complex network (Fig. 1), suggesting that these 489 proteins on exosomal delivery to a target cell can readily interact with the 998 cellular proteins that might influences molecular function of the target cell. Thus, the sub-network was constructed to capture the proteins from SET-A that have a high degree of interaction with the proteins (998) of the extended network. Notably, Mothers against decapentaplegic homolog 2 (SMAD2) from the SET-A showed high degree of interaction in the protein network (Fig. 2). Consequently, the gene enrichment analysis was carried out with 1,487 proteins. The top 20 enriched GO terms were shown in Fig. S1. Likewise, the top 20 enriched molecular pathways were illustrated in Fig. S2. On the assessment of sub-network (Fig. 2), SMAD2 protein was noticed with the highest interacting proteins (n = 174) that represented as putative target for molecular docking.

Molecular docking with SMAD2

The protein structure of SMAD2 was retrieved from PDB (ID: 5XOD) and optimized to perform a docking with 27 natural compounds along with five references CAD drugs (amlodipine, furosemide, digoxin, clopidogrel, and atorvastatin). The active site for the SMAD2 was determined as GLU270, PRO271, PHE273, TRP274, TYR340, GLY342, GLY343, VAL345, PHE385, ASN387, LEU442, GLY444, and PRO445 using P2Rank tool. The compound inermin (MOL003648) has the highest binding score (−5.024 kcal/mol). Likewise, furosemide showed highest affinity of −4.472 kcal/mol with SMAD2 compared other analyzed CAD drugs. These results indicate that the compound from Panax ginsenghas a higher affinity towards SMAD2 than the CAD drugs. Top 10 molecular docking results were represented in Table 1 and the interaction of the inermin and furosemide with SMAD2 were displayed as Figs. 3A & 3B. Overall docking results of 27 natural compounds and five references CAD drugs were given in the File S2.

Molecular dynamic simulation of SMAD2—inermin complex

The docked Panax ginseng compound inermin and the known CAD drug furosemide with SMAD2 protein were used for MD simulation. The RMSD plot shows the inermin-SMAD2 and furosemide-SMAD2 complexes were tend to deviate initially and stabilized at the end of the simulation (Figs. 4A & 4B). The average RMSD value of the inermin-SMAD2 complex was 2.67 Åand the average of furosemide-SMAD2 complex was 3.69 Å. From the RMSD trajectory, the inermin has value less than 3 Å, which indicates inermin-SMAD2 formed a stable complex. Based on the RMSF plots, the inermin-SMAD2 complex (Fig. 5A) shows lower fluctuation than the furosemide -SMAD2 complex (Fig. 5B) during the 100 ns MDS run. Inermin (Figs. 6A & 7A) and furosemide (Figs. 6B & 7B) had 35 and 22 ligand contacts with the SMAD2 target. Further, the inermin average values of Rg, SASA and molecular surface area were 3.44 Å, 258.7 Å2, and 247.75 Å2, respectively (File S3) (Fig. 8A). Likewise, the average Rg, SASA, and molecular surface area trajectories for furosemide properties were 3.82 Å, 211.91 Å2, and 269.55 Å2 (Fig. 8B) (File S3).

DFT quantum analysis—molecular geometry parameters

The optimized structure of inermin was shown in Fig. 9A. The bond length shows the distance between the nuclei of one atom and the other atoms. The angle between the two atomic bonds was determined by the bond angle analysis. The highest distance between the nuclei was noted in C10-C11 (1.56 Å), while O21-H33 (0.95 Å) has the shortest distance. The carbon-hydrogen bond has the lowest bond length compared to other atoms of inermin. In contrast, the carbon–carbon bond has the highest bond length, which shows that carbon interaction withholds high binding capacity. The lowest and the highest bond angle of inermin were observed between O21-H33 (0.95°) and C10 -C11 (1.56°). The lower torsion angle noted in C1-C2-C9-H25 (−179.91°), while greatest torsion angle observed in C15-C17-C18-O21 (179.94°). DFT analysis output data was provided in the File S4.

DFT quantum analysis—HOMO and LUMO

The highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) are calculated with DFT-B3LYP6-311++ (dp). The frontier molecular orbital of inermin was illustrated in Fig. 9B. The HOMO orbital of inermin was partially covered without hydrogen atom, while the LUMO orbital fully dispersed except for the hydrogen atoms. The chemical activity of the molecule was calculated by correlating bond energies and their associated quantum parameters. The molecular behavior of the chemical compound was determined by the quantum parameters (ionization potential, electrophilicity (ω), electronegativity (χ), softness (S), hardness (η), and chemical potential (μ)). The electrophilicity scale shows that the atom has the ability to take up electrons. Likewise, the electronegativity scale describes the ability of the atom to attract electrons toward it. The hardness and softness values represent the electron charge transfer resistance and acceptance of an atom. The HOMO value was found to be −5.40 eV, while the LUMO was −0.45 eV, and has an energy gap (ΔE) of 4.95 eV (Fig. 9B). Furthermore, the calculated quantum parameters showed ionization potential (I) was 5.40 eV, electron affinity (A) was 0.45 eV, electronegativity (χ) was 2.93 eV, chemical hardness (η) was 4.95 eV, chemical potential (μ) was −2.93 eV, electrophilicity index (ω) was 1.46 eV, and chemical softness (σ) was 0.20 eV. Overall, the results suggest that inermin has high reactivity by quantum chemical analysis.

DFT quantum analysis—electrostatic potentials

The B3LYP/6-31G**++ in the Jaguar suite was utilized for ESP analysis. Inermin was analyzed to identify the reactivity, hydrogen bonding, and structural activity. The color variation denotes the electrostatic potential (Fig. 9C). Electrorophilic and nucleophilic attack regions were identified by the electrostatic potential. Particularly, blue, red, and white color indicates the positive, negative, and neutral electrostatic region. Interestingly (Fig. 9C), the major hydrogen atoms with high positive charge are prone to nucleophilic attack, which might be favorable for the protein-ligand interaction. The oxygen atoms present in inermin are negatively charged, which are expected to have nucleophilic attack. The higher electron density (deep red) was noted in the oxygen atom O21, which is likely to form a hydrogen bond with the protein (Fig. 9C).

DFT quantum analysis—Mulliken charges

The molecules capacity of the binding was determined by the electronic charges. The Mulliken charges of the representative atoms of inermin were shown in Fig. 9D. Surprisingly, all the hydrogen atoms are positively charged and the ones facing adjacent side to oxygen are highly positive compared to other hydrogen atoms. The electronegativity of the oxygen atom contribute strong interaction with the hydrogen atoms that ultimately supports the lead molecule to bind to protein efficiently (Fig. 9D). Likewise, the carbon atom in phenol ring interacted with oxygen atom that has a positive charge, while the others were negatively charged. The top three hydrogen atoms O21, O12, and O14 each had charge of −0.55, −0.52, and −0.51, which shows a negative value that might attract other atoms. In contrast, the C18 has the greater positive charge of 0.34 (Fig. 9D).