This article is mentioned in:

Introduction

Spinal cord injury (SCI), an event associated with permanent neurologic deficit (1,2), is reported to affect 10–40/million individuals annually in developed countries (3). In the United States, there are ~253,000 cases of SCI, with ~11,000 new cases diagnosed annually (3). The economic cost of this disease in the United States alone is estimated to be billions of dollars (4). Violence and accidents are the two causes of SCI, with 26% of cases associated with violence and the remaining associated with accidents (5).

The first mention of SCI is recorded to be in 1700 BC in the Edwin Smith Papyrus, an Egyptian medical text, describing this disease as an ‘ailment not to be treated’, which demonstrates the severity of SCI as it was considered an untreatable ailment (6). SCI is a severe condition, with mortality occurring in the majority of patients prior to any primary hospital care, while patients treated at the hospital are prone to morbidity and mortality (5). Neurological deficit and disabilities resulting from SCI not only affect the sensory and motor capabilities, but also has a strong impact on the physiological condition of the patient (7). SCI is categorized into primary SCI, defined as the damage inflicted at the time of trauma, and secondary SCI, which is defined as the body's response to the initial trauma (8).

Due to the significant decrease in the quality of life of patients following SCI, investigation of the pathophysiology and treatment of this disease is urgent. To date, studies have focused on preventing secondary SCI, promoting regeneration and replacing the destroyed spinal cord tissue (9,10). Methylprednisolone is the only current agent approved for the treatment of primary SCI that is widely used, although it presents selective/limited efficacy and significant side effects (11–15). In order to overcome the side effects of this glucocorticoid, the use of poly (lactic-co-glycolic acid) nanoparticle formulations for the encapsulation of methylprednisolone has been used (16). Novel therapeutic strategies are being developed based on the basic approach involving the targeting of the cascading mechanism that leads to secondary SCI (10). In this regard, altering the neuro-inflammation (17), reducing free-radical damage (18), reducing excitotoxic damage to neurons (19), improving the blood flow to the primary injury site (20) and countering the effect of ionic exchange at the primary injury site (21) are the most commonly studied therapeutic strategies.

Estrogen and its analogs have been widely investigated as secondary SCI inhibitors, with previous studies suggesting their anti-apoptotic effects, as well as their role in decreasing the activation of cysteine protease, in the attenuation of vascular endothelial growth factor and in aquaporin upregulation (21–24). Estrogen has been reported to aid in the improvement of SCI and promotion of the repair of the damaged spinal cord by modulating the immune response (25). Thus, studies have been conducted on boosting and/or modulating the immune response, and one such approach involved the inhibition of the Ras/phosphatidylinositol 3-kinase (PI3K)/RAC α serine/threonine-protein kinase (Akt)/serine/threonine-protein kinase mTOR (mTOR) signaling pathway (22,26,27).

The involvement of Ras in the sympathetic neuronal survival via the PI3K and dual specificity mitogen-activated protein kinase kinase mek-1 signaling pathway is well established, while its role in suppressing apoptosis by regulating the cellular tumor antigen p53 pathway is known to benefit SCI treatment (28–30). The majority of studies to date have examined the Ras/PI3K/Akt/mTOR pathway as the neuroprotective and neuroregenerative functions in SCI are regulated by this pathway, and rapamycin inhibition of mTOR has been reported to exhibit a positive effect on mice with SCI (31–33). Although the role of the Ras/Raf/extracellular signal-regulated kinase 1/2 (ERK1/2) signaling pathway in cell proliferation, migration, differentiation and death is well documented (34–36), its involvement in the treatment of SCI has not been widely investigated. Studies have revealed that the U0126 inhibitor functions as a potential drug for restoring SCI by affecting spinal cord neurons (SCNs) via altering the Raf/ERK signaling pathway (26,37–39).

The present study attempted to identify a novel natural compound (Lead) involved in the alteration of the Raf/ERK pathway for the regeneration of SCNs following SCI. The study used the Interbioscreen (IBS) natural compound database (www.ibscreen.com) consisting of 48,531 natural compounds, known to be largest collection of natural compounds, their derivatives and mimetics worldwide. A multistep structure-based virtual screening technique, followed by docking and simulation, was employed for identifying a promising lead compound, which is a widely used method for identifying such compounds (40–42).

Materials and methods

Protein and compounds of natural origin dataset preparation

The atomic coordinates for the Ras [Protein Data Bank (PBD) ID, 4LPK] structure required for the development of a lead compound were obtained from the PDB (rcsb.org/pdb/home/home.do). This crystal structure of GTPase KRas is GDP-bound, and the coordinates of this Ras protein were used. The energy-minimization of the structure was performed by using the Swiss-PdbViewer software (version 4.1; Swiss Institute of Bioinformatics, Lausanne, Switzerland). The root mean square deviation (RMSD) was monitored using the GROMOS96 43b1 force field (43). Approximately 48,531 natural compounds from the IBS database were used for targeting the GTP binding site of Ras.

Protein network reconstruction

The Ras/PI3K/Akt/mTOR protein-protein interactions networks were constructed using Cytoscape (version 3.5.1; National Resource of Network Biology, San Francisco, CA, USA). The tool was used for compiling protein functional classification, and physical and genetic interaction networks.

Virtual screening and drug-likeliness prediction

The ArgusLab suit was used for virtual screening of compounds of natural origin (44). In total, >1% (414) compounds were shortlisted according to the binding energy (∆G) calculations. A ∆G of −7 kcal/mol was set as the cut off to obtain the initial subset of compounds. This cut off value was a limiting bias set in order to limit the number of compounds, since a ∆G of −8 kcal/mol provided very few compounds for further investigation and a value of −6 kcal/mol provided a very large number of compounds. Thus, ∆G was set at −7 kcal/mol to provide an acceptable number of compounds for further validation (40–42). The number of selected compounds was then further reduced by subjecting to the rules set by Lipinski (45), which assesed the molecular weight (MW), hydrogen bond acceptor (HBA), hydrogen bond donor (HBD) and octanol/water partition coefficient (cLogP) of each compound. Lipinski's rule of five (RO5) parameters provided five compounds that were subjected to further analysis. Absorption, distribution, metabolism and excretion (ADME) analysis was performed using the PreADME/T online server (preadmet.bmdrc.kr/adme/).

Molecular docking analysis

A structure-based drug designing method was used by employing the AutoDock 4.2 tool for molecular docking analysis (46). This tool calculates energy values by classification into internal and torsional free energy. The internal energy is the sum of the desolvation, hydrogen bonding, van der Waals and electrostatic energies. The Lamarckian genetic algorithm default parameters (46) were used for calculating the ∆G of each shortlisted compound. A grid box (40×40×40 Å) was built around the GTP binding site. Energy values were generated and the binding mode with the GTP binding site was used to limit the investigated compounds to a single molecule.

Molecular dynamics simulation (MDS)

The single compound shortlisted was simulated for 100 nsec using MDS, which was conducted using the GROMACS suit (47). From the compounds initially identified, molecular docking identified the optimal compound as lead3, which was then studied using MDS. To mimic in vitro conditions, the Ras-Lead3 complex was kept in an in silico neutral (no charge) water bath, where a water molecule was represented by a simple point charge (SPC216). GROMOS 43a1 force field for Ras was used for simulation, and the force field for the compound was calculated using the PRODRG server (48). Energy minimization procedures for 1 nsec were performed using canonical [number of particles, volume and temperature (NVT) are held constant] and isothermal-isobaric [N, pressure and T (NPT) are held constant] ensembles. In NVT and NPT ensembles the coupling scheme of Berendsen, SHAKE algorithm and particle mesh Ewald method were used (49–51). Next, the energy-minimized Ras-Lead3 complex was simulated for 100 nsec, and the trajectories generated were subjected to Molecular Mechanics Poisson-Bolzmann Surface Area (MM-PBSA) calculations. The g_mmpbsa tool developed for GROMACS, which was used for principal component analysis (PCA) (52). Snapshots of the coordinates and the total energies were obtained after 500 psec, while 501 snapshots of the RAS-Lead3 complex were also subjected to the MM-PBSA calculation. The binding free energy (ΔGbind) was composed of the following species:

In the aforementioned equations, Gcomplex, Gprotein and Gligand represent the free energy of the respective species. ΔEMM refers to the gas phase energy, ΔGsol is the solvation energy, TΔS represents an entropy term, ΔEval is the sum of the internal energy of bonds, angle and torsion, ΔEele is the electrostatic interaction, ΔEvdw is the van der Waals interaction energy, ΔGp is the polar salvation free energy, and finally ΔGnp refers to the nonpolar salvation free energy.

Molecular visualization analysis

The RAS-Lead3 complex was investigated using the visualization tools Pymol (version 1.8.6.2) (53) and Discovery Studio (version 16.1.015350) (54). For graph construction, the Grace Program (version 5.1.22) (55) and Gnuplot (version 5.0.6) (56) were used.

Results and Discussion

A database comprised of natural compounds, their derivatives and mimetics were considered in order to identify compounds that would inhibit the Ras/Raf/ERK1/2 pathway by inhibiting the GTP binding site of Ras (Fig. 1). The role of the Ras/Raf/ERK signaling pathway in promoting apoptosis in neural cells in SCI is well established (27). Ras inhibitors modulate the intracellular signaling of Ras to exhibit a pro-apoptotic or neuroprotective effect (31). The U0126 inhibitor demonstrates a neuroprotective activity by suppressing ERK1/2 activation (57). The consistent evidence on the role of ERK1/2 in neuronal apoptotic cell death provided an opportunity to identify novel small molecule inhibitors using in silico approaches (37). This signaling cascade is triggered by the membrane receptor, which allows Ras to swap GDP for a GTP and to become active. This activated Ras then activates the kinase activity of Raf, which in turn phosphorylates and activates ERKs (ERK1 and ERK2). In the present study, in order to identify a novel Ras inhibitor, virtual screening, drug-likeliness, molecular docking analysis and MDS methods were used, and this methodology is demonstrated in Fig. 2. Virtual screening helped to limit the number of compounds from 48,531 natural products to 414 compounds using a limiting bias of ∆G of −7 kcal/mol.

Figure 1. Protein subnetworks of Ras (the query protein) with Raf and ERK1/2. The subnetworks were generated by Cytoscape software. ERK1/2, extracellular signal-regulated kinase 1/2. KRAS, GTPase KRas; BRAF, serine/threonine-protein kinase B-raf; MAP2K, dual specificity mitogen-activated protein kinase kinase; SHOC2, leucine-rich repeat protein SHOC-2; RAF1, RAF proto-oncogene serine/threonine-protein kinase; ARAF, serine/threonine-protein kinase A-Raf; SOS1, son of sevenless homolog 1; RGL2, DELLA protein RGL2.

Figure 2. Methodology of the present study, attempting to identify the inhibitor of Ras protein for the downregulation of Ras/Raf/extracellular signal-regulated kinase 1/2 signaling pathway. IBS, Interbioscreen; ∆G, binding energy.

In order to focus on compounds that may be promising for further development, each of the top identified compounds was examined for drug-likeliness. The drug-likeliness of the shortlisted compounds was defined according to the mutagenic and carcinogenic properties, Lipinski's RO5 and total polar surface area (TPSA). The RO5 properties included the number of HBDs and HBAs, the MW and cLogP, with the permissible range being HBD≤5, HBA≤10, MW≤500 Da and cLogP≤5. Fig. 3 represents the three-dimensional point plot of the HBA, HBD and cLogP values. The coloring in the figure is according to the TPSA.

Figure 3. Three-dimensional point plot of the IBS natural product database. The point plot represents the natural products at the virtual screening, drug-likeliness, molecular docking and molecular dynamics simulation stages of the study, demonstrating the final selection of five lead compounds according to Lipinski's rule of five parameters. IBS, Interbioscreen; tpsa, total polar surface area; cLogP, octanol/water partition coefficient.

Table I demonstrates the drug-likeliness properties of the top five identified compounds. The drug-likeliness values of these compound conformed to the values expected from typical drugs. These five compounds where subjected to ADME analysis. In silico methods were used to calculate the Caco-2 cell permeability (Caco-2 p), MDCK cell permeability, human intestinal absorption, plasma protein binding (PPB) and blood-brain barrier (BBB), and the results generated are shown in Table II. The ADME properties of the shortlisted compounds were within the acceptable range for known drugs available in the market.

Table I. Drug-likeliness of top five lead compounds.

Table I.

Drug-likeliness of top five lead compounds.

IBS no.	MW	cLogP	HBD, n	HBA, n	TPSA	Mutagenicity	Carcinogenicity
05678	314.419	4.41	1	3	38.69	No	No
46780	405.401	3.35	2	7	109.57	No	No
49515	413.183	3.65	2	7	109.75	No	No
49817	405.151	3.94	3	8	129.28	No	No
64118	405.121	3.12	2	6	96.28	No	No

[i] cLogP, octanol/water partition coefficient; TPSA, total polar surface area; IBS, Interbioscreen; HBA, hydrogen bond acceptor; HBD, hydrogen bond donor.

Table II. Absorption, distribution, metabolism and excretion properties calculated by in silico approach.

Table II.

Absorption, distribution, metabolism and excretion properties calculated by in silico approach.

Name	IBS no.	Definition	Caco-2 p (nm/sec)	MDCK p (nm/sec)	HIA (%)	BBB (c.b/c.bl)	PPB PPB (%)
Lead1	64118	N-(3-hydroxypropyl)-3-(5-methyl-7-oxo-3-phenyl-7H-furo[3,2-g]chromen-6-yl)propanamide	20.79	0.14	94.97	0.052	91.18
Lead2	46780	2-(3-(5-methyl-7-oxo-3-phenyl-7H-furo[3,2-g] chromen-6-yl)propanamido)acetic acid	19.14	0.09	97.82	0.012	93.02
Lead3	05678	(13S)-13-methyl-6,7,8,9,11,12,13,14,15,16-decahydrospiro[cyclopenta[a]phenanthrene-17,2′-[1,3]dioxolan]-3-ol	43.61	56.29	95.62	3.620	100.00
Lead4	49515	4-(2-(3-(tert-butyl)-5,9-dimethyl-7-oxo-7H-furo [3,2-g]chromen-6-yl)acetamido)butanoic acid	20.02	0.12	96.89	0.012	89.00
Lead5	49817	(S)-2-(2-(3-(tert-butyl)-5-methyl-7-oxo-7H-furo [3,2-g]chromen-6-yl)acetamido)-3-(4-hydroxyphenyl)propanoic acid	21.23	0.05	95.97	0.019	89.61

[i] Absorption, distribution, metabolism and excretion properties of top five identified compounds were calculated by the PreADME/T online server. Caco-2 p, Caco-2 cell permeability; MDCK p, MDCK cell permeability; HIA, human intestinal absorption; PPB, plasma protein binding; BBB, blood-brain barrier; IBS, Interbioscreen; c.b., concentration in brain; c.bl., concentration in blood.

In the current study, the AutoDock tool was used for molecular docking simulations, selecting the top binding pose based on the ∆G for further analysis (Fig. 4). Each binding pose was analyzed using the Discovery Studio software, and the default parameters were used to calculate all the possible interactions. The results generated are listed in Table III, where the ∆G, binding pocket and number of hydrogen bonds formed are demonstrated. The interactions investigated included the van der Waals, conventional hydrogen bond, carbon hydrogen bond, π-cation, π-donor hydrogen bond, alkyl and π-alkyl interactions. The Lead1-Ras complex presented a ∆G of −6.12 kcal/mol and the interaction is shown in Fig. 5, where Lead1 is forming two conventional hydrogen bonds with ALA146 and ASN116 of the Ras GTP binding site. The binding pocket of Lead1 comprises following amino acids ALA59, GLY12, PRO34, GLY13, VAL29, PHE28, LEU120, LYS117, ASN116, LYS147, GLY15, ALA18, SER145, ALA146, SER17, ASP57, LYS16 and THR58. The O23 atomic site of Lead1 shows hydrogen bond interaction with ALA146 and ASN 116, with a distance between the Lead and Ras of 2.14 and 1.97 Å, respectively.

Figure 4. Solid ribbon structure of the GTP binding pocket of Ras with the top five lead compounds identified.

Figure 5. Two-dimensional representation of the binding pocket of the top five identified compounds with Ras.

Table III. Binding pose analysis of top five identified compounds, including calculation of the ligand-binding pocket and the hydrogen bond formation using Discovery Studio.

Table III.

Binding pose analysis of top five identified compounds, including calculation of the ligand-binding pocket and the hydrogen bond formation using Discovery Studio.

Name	IBS no.	ΔG (kcal/mol)	Ligand binding pocket	H-bonds
Lead1	64118	−6.27	ALA59, GLY12, PRO34, GLY13, VAL29, PHE28, LEU120, LYS117, ASN116, LYS147, GLY15, ALA18, SER145, ALA146, SER17, ASP57, LYS16 and THR58	ALA146:HN-Lead1:O23 (2.14 Å) ASN116:HD21-Lead1:O23 (1.97 Å)
Lead2	46780	−7.72	PRO34, TYR32, GLU31, ASP30, PHE28, ALA18, ASN116, LEU120, LYS117, GLY15, ASN85, VAL14, SER17, LYS16 and GLY13	TYR32:HN-Lead2:O12 (2.02 Å) ASN85:HD21-Lead2:O25 (1.88 Å)
Lead3	05678	−8.58	GLY60, THR58, SER17, GLY15, TYR32, PHE28, VAL29, ASP30, GLU31, ALA18, LYS117, GLY13, LYS16, PRO34, GLY12 and ALA59	GLY60:HN-Lead3:O21 (2.11 Å)
Lead4	49515	−5.47	TYR32, GLU31, PRO34, GLY12, ASN86, ASN85, GLY13, LYS117, LYS116, SER17 and ASP33	ASN85:HD21-Lead4:O20 (2.07 Å)
Lead5	49817	−6.08	ASP119, LEU120, LYS117, ASN116, GLY15, ALA18, TYR32, PRO34, SER17, GLY13, VAL29, PHE28, GLU31, LYS147, PHE28 and ASP30	Lead5:O32-SER17:OG (2.58 Å) SER17:HN-Lead5:O32 (1.68 Å)

[i] The atoms involved in hydrogen bond formation are mentioned with their position in subscript. Amino acids shown in bold represent those involved in the formation of hydrogen bond with the compounds. IBS, Interbioscreen; OG, Oxygen group; HD, Hydrogen donor.

Lead2 exhibited two conventional hydrogen bond interactions with the GTP binding pocket of Ras. Two atoms of Lead2, O12 andO25, formed a bond with TYR32 and ASN85, with a bond length of 2.02 Å and 1.88 Å, respectively. In addition, the binding pocket of Lead2 was comprised of 15 amino acids, as follows: PRO34, TYR32, GLU31, ASP30, PHE28, ALA18, ASN116, LEU120, LYS117, GLY15, ASN85, VAL14, SER17, LYS16 and GLY13. Furthermore, Lead3 exhibited a ΔG of −8.58 kcal/mol, which was the best reported value among the top five shortlisted compounds. The binding pocket of the Lead3 molecule with a GTP binding site of Ras comprises of the following amino acids: GLY60, THR58, SER17, GLY15, TYR32, PHE28, VAL29, ASP30, GLU31, ALA18, LYS117, GLY13, LYS16, PRO34, GLY12 and ALA59. Out of these, Lead3 formed hydrogen bond with GLY60, while the interaction of interest in the present study was formed by the Lead3 O21 position.

Based on the ΔG and the number of interactions, Lead4 was ranked last among the top five shortlisted compounds. Its binding pocket was comprised of 11 amino acids, namely: TYR32, GLU31, PRO34, GLY12, ASN86, ASN85, GLY13, LYS117, LYS116, SER17 and ASP33. Lead4 presented a ΔG of −5.47 kcal/mol, and had a single hydrogen bond interaction between the Lead4 O20 position and Ras ASN85 (2.07 Å). The last Lead compound identified in the current study was termed Lead5, which demonstrated a ΔG of −6.08 kcal/mol and formed a maximum number of interactions in the GTP binding pocket of Ras. The binding pocket of Lead5 was comprised of ASP119, LEU120, LYS117, ASN116, GLY15, ALA18, TYR32, PRO34, SER17, GLY13, VAL29, PHE28, GLU31, LYS147, PHE28 and ASP30.

The analysis in the present study observed that Lead3 was compound with the best BBB and PPB, as well as the lowest ΔG value, and was thus selected for further analysis. In order to explore the stability and binding mode of Lead3 with the GTP binding site of Ras, MDS was performed for a 100 nsec run under GROMOS 43a1 force field on the energy-minimized Lead3-Ras complex. The comparative RMSD plot was analyzed to examine the stability of this complex over time, and this plot revealed a substantial decrease in the RMSD of the complex throughout the simulation (Fig. 6). Furthermore, the dominant motions of the Lead3-Ras complex were also investigated using the PCA tool. The free energy landscape was plotted in order to retrieve the lowest energy conformer (t=67,324 psec). Subsequently, the retrieved structure of the complex was used investigate the interactions at that time point, as shown in Fig. 7. The MM-PBSA calculations revealed that, in the GTP binding site, GLY60 was the best binding site with Lead3 over time, thus revealing the inhibitory property of Lead3. Table IV demonstrates the ∆G of the simulated complex.

Figure 6. RMSD plot of Ras, the Ras-Lead3 complex and Lead3. Ras (black), Ras-Lead3 complex (red) and Lead3 compound (green). RMSD, root-mean-square deviation.

Figure 7. Free energy landscape (inverted) projection of the principal component analysis of the first (PC1) and second (PC2) principal eigenvectors for the backbone of Ras-Lead3 complex. The lowest energy conformer is presented, and the structure retrieved demonstrates a respective hydrogen bond pattern with atomic details. PC, principal component; OXT, terminal oxygen atom; CA, carbon atom.

Table IV. Molecular mechanics poisson-bolzmann surface area calculations.

Table IV.

Molecular mechanics poisson-bolzmann surface area calculations.

Binding energy	Values (kJ/mol)
Van der Waal	−142.32±16.21
Electrostatic	−6.51±1.31
Polar solvation	44.87±7.89
Solvent accessible volume	−102.32±18.61
Binding	−189.33±28.22

In patients with SCI, the Ras/Raf/ERK1/2 signaling pathway has been associated with the promotion of apoptosis in neural cells (39); small molecule inhibitors and other methods of modulating the Ras/Raf/ERK1/2 signaling pathway have become a promising target for the development of novel treatments for traumatic injuries to the nervous system. To the best of our knowledge, the present study is first to use structure-based drug design to identify a novel lead drug of natural origin as a small molecular inhibitor for modulating Ras/Raf/ERK signaling pathway. Inhibitors, including PD98059 and U0126, have been used in in vivo studies targeting the Ras/Raf/ERK1/2 signaling pathway (58).

In conclusion, the in silico approach used in the present study investigated the IBS database of ~50,000 compounds and identified a single compound that bound to the GTP binding site of Ras, thus downregulating of the Ras/Raf/ERK1/2 signaling pathway. The inhibition of this pathway is able to restore the neuronal migration, adhesion and dendritic spine development. A chemical-computational approach was adopted in the current study, focusing on the development of Ras inhibitors of a natural origin that may prevent neuronal apoptosis in SCI, using used computer-aided drug designing approaches. The Ras GTP inhibitor proposed in the present study can be developed into a drug in further studies; thus, these findings provide a starting point for the development of a drug that may have a potential effect in the treatment of SCI. However, further in vitro and in vivo studies are required to translate the proposed lead compound into a potential drug.

References