A subset of 137,972 natural products obtained from the ZINC12 database (Sterling and Irwin, 2012) underwent preparation procedures at a pH of 7.4, mirroring the physiological pH found in the human body. This procedure included the generation of every possible tautomer and ionized state (Release Schrödinger Suite, 2012). Following this, a meticulous screening protocol was implemented to assess each compound based on predefined criteria associated with drug-likeness properties, such as molecular weight (MW) not surpassing 500 Da, a LogP value of 5 or lower, a maximum of 10 hydrogen bond acceptors (HBA), no more than 5 hydrogen bond donors (HBD), a limitation of 10 rotatable bonds (RB), and a total polar surface area (TPSA) constrained to 140 Å² or less. Following this screening, 15,000 compounds that met these criteria were screened.

The ligands used as input for the docking investigation were obtained from the virtual screening hits of the natural databases. Subsequently, ligands were integrated into the workstation, and energy was minimized by utilizing the OPLS3e force field in the LigPrep module of the program. This minimization aids in assigning bond order, adding hydrogens to ligands, and converting 2D structures to 3D structures for performing docking studies. The generated output file, including the best conformations of the ligands, was then used for docking studies (Schrödinger Release, 2017).

The MAPK^ERK protein’s crystal structure (PDB ID: 4qte) underwent refinement via Schrödinger’s protein preparation wizard (Protein Preparation Wizard, 2017). This meticulous procedure ensured an elevated level of confidence in the accuracy of the protein structure. By converting it from its raw state to a refined form, the protein was ensured for subsequent molecular docking and dynamic studies, enhancing its suitability and reliability for comprehensive studies. The protein preparation process consisted of correcting bond orders, removing water molecules and other non-specific components from the crystal structure, and adding hydrogen atoms to the protein structure to modify the tautomeric and ionization states of amino acid residues. After adding the missing hydrogen, the system underwent restricted energy reduction using the OPLS 2005 force field to ensure high accuracy.

The MAPK^ERK protein and its prepared ligand underwent docking using Schrödinger’s GLIDE module (Glide, 2017). This process aimed to improve binding affinity by aiding in the identification of specific structural motifs. The docking results were analyzed using XP descriptors, which provided useful details on a variety of intermolecular interactions. The Glide XP descriptor data facilitated the determination of energy associated with each pose, aiding in the optimization of the ligand conformation within the ligand-receptor complex. Before initiating the docking protocol, Schrödinger’s Glide Grid generation was employed to create a specific Glide Grid surrounding the co-crystallized ligand. Subsequently, this grid served as the framework for the docking procedure, facilitating the precise positioning and interaction analysis between the prepared ligand and the protein structures within the designated spatial coordinates.

The physico-chemical parameters were utilizing the Swiss ADME service (Daina et al., 2017) and ADMET lab 2.0 (ADMETlab 2.0, 2021). The anticipated parameters encompass several key molecular attributes: MW, HBD, HBA, the octanol/water partition coefficient (log P), aqueous solubility (Log S), predicted apparent Caco-2 cell permeability, human intestinal absorption propensity (HIA), plasma protein binding (PPB), and the tally of RB. These parameters collectively serve as crucial indicators for assessing the physicochemical properties, drug-likeness, and potential bioavailability of the compounds under investigation.

Validation of the docking study is necessary to ensure precise molecule docking. The precision of the molecular docking process is assessed using the Root Mean Square Deviation (RMSD). In order to calculate the RMSD, the co-crystallized ligand of 4qte was removed from it and then subjected to XP glide docking using the matching receptor grid. The obtained RMSD value, which was less than 2.0 Å, indicates the reliability of the docking study (Gheidari et al., 2023).

MD simulation was conducted utilizing Desmond through the Schrödinger-Maestro interface (Desmond, 2017). This enables the precise identification and prediction of ligand-receptor interactions and the validation of molecular docking results. MD simulations were performed by solving a substance in an explicit orthorhombic water box using the SPC water model. Sufficient Cl ions were added to the system in order to counterbalance the overall charge of the complex. The duration of the simulation was 100 ns. The NPT ensemble was utilized, maintaining a constant number of atoms, a pressure of 1.01325 bar, and a temperature of 300 K. The default thermostat was the 1.0‐ps interval Nose-Hoover chain method, and the default barostat was the 2.0‐ps interval Martyna-Tobias-Klein. The maestro simulation interaction diagram was used to assess the MD simulation.

The DFT is used to ascertain the electron’s density and energy properties. The Gaussian 09W program performs calculations that illustrate the estimated structure of atoms, molecules, crystals, and surfaces, as well as their interactions (Frisch et al., 2009). The wavenumbers of the vibrations were determined by the use of the B3LYP method and a 6-31++G (d,p) basis set. The B3LYP functional is a valuable approach for accurately characterizing harmonic vibrational frequencies in molecules of small to medium size. The output verification files were analyzed using GuassView 6.0. Molecular orbital (MO) analysis is essential in quantum chemistry and has been used to thoroughly characterize chemical behavior. The highest occupied molecular orbital (HOMO) and lowest unoccupied molecular (LUMO) of a molecule are employed to characterize chemical properties, encompassing reactivity, stability, kinetics, hardness, softness, and electronegativity.

The optimization of ligands (L ₁ -L ₁₈ ) was initially conducted using the B3LYP/6-31++G (d,p) basis set in the gas phase, and the results obtained from this optimization are presented in Table 1.

The geometries of the selected ligands underwent optimization until reaching the lowest energy gradient, confirming their status as true local minima, as no imaginary frequencies were detected. Figure 2 illustrates the optimized structures of the selected ligands.

The study of MO is very important in quantum chemistry, significantly enhancing our understanding and knowledge of chemical behavior. The HOMO and LUMO that are the chief molecular orbitals in a ligand are shown in Figure 3. The red and green color distributions represent the positive and negative phases, respectively, in the MO wave function.

The parameter hardness (h) quantifies the degree of hardness or softness exhibited by a molecule. Greater molecular softness correlates with enhanced reactivity. The atom’s electronegativity (X) indicates its ability to attract electron pairs towards its nucleus. The ligand L ₄ , exhibiting the smallest HOMO–LUMO energy gap value of 2.4884 eV, suggests a heightened potential for chemical reactivity. Also, L ₄ demonstrates a low hardness value of 1.2442, positioning it as the softest ligand among all ligands. The greater electronegativity value of ligand L ₄ indicates that L ₄ has a strong ability to attract electrons and operates as a superior electrophile compared to other ligands. Ligands L ₃ and L ₁₀ exhibited favorable reactivity after L ₄ , with energy gap values of 3.7621 eV and 4.7434 eV, respectively. Similarly, after L ₃ , the ligands L ₁₁ , L ₄ , and L ₅ exhibit significant polarizability, with corresponding values of 236.316, 239.453, and 240.025. The energetic parameters of ligands (L ₁ -L ₁₈ ) are shown in Table 2.

Schrödinger’s Ligand docking module utilizes grid-based ligand docking with energetics (GLIDE) (Glide, 2017) in three specific docking modes: high-throughput virtual screening (HTVS), standard precision (SP), and extra precision (XP) for docking and scoring. 15,000 compounds underwent screening by HTVS docking and scoring. A total of 50 ligands that had the highest G-Score with HTVS were further evaluated using SP docking and scoring. This was done to ensure that the docking postures were reliable and proper. To decrease false-positive findings, 18 ligands that were completely thriving were screened using XP docking and scoring. An analysis was conducted on the interactions of these protein-ligand complexes, and the outcomes are shown in Table 3. MAPK^ERK’s active site’s amino acid residues, both bonding and non-bonding interactions with ligand L ₁ –L ₁₈ , included Asn154, Ser153, Asp111, Tyr113, Gly32, Glu33, Gly34, Asp167, Leu156, Ile31, Cys166, Val39, Lys54, Tyr30, Lys114, Gln105, Ile84, Lys151, Asp111, Asp149, Ala35, Tyr36, Glu71, Ala52, and Lys117.

L ₁₅ ranked second out of the three top-scoring ligands due to its strong binding interactions with the active site of the MAPK^ERK protein. Ligand L ₁₅ has six carbon-hydrogen bonds, corresponding to the amino acid residues Ser153 and Glu33, with bond lengths of 2.42 and 2.74, respectively. Additionally, it interacts with Asp167 at a bond length of 3.02 and with Asn154 at bond lengths of 2.18 and 2.23. Furthermore, five hydrophobic interactions have been observed between the ligand L15 and the MAPK^ERK protein. These include bond lengths of 4.88 Å for Lys54, 3.82 Å for Cys166, 4.56 Å for Leu156, and two band lengths of 4.60 and 4.85 Å for Val39. Also, ligand L ₁₅ interacts with Gly134, Lys151, Gln105, Ile84, Ala52, Ile31, Gly32, and Asp111 in five van der Waals interactions. Figure 5 displays the 2D and 3D binding interactions of L ₁₅ within MAPK^ERK’s active pocket.

Ligand L ₁₇ has five carbon-hydrogen bonds, two of which are associated with amino acid Ser153 and have bond lengths of 2.21 and 2.79, respectively; the other two are related to amino acid Asn154 and have bond lengths of 2.31 and 2.72, respectively; and another is related to Glu33 and has a bond length of 2.48. In addition, it has four alkyl interactions with the residues Leu156, Cys166, and Val39, and nine van der Waals interactions with Gly34, Asp167, Gln105, Asp111, Lys54, Gly32, Lys151, and Ile31. Figure 6 displays the 2D and 3D binding interactions of L ₁₇ within MAPK^ERK’s active pocket. In summary, Table 3 illustrates that the two main interactions in all complexes are hydrophobic intermolecular interactions and carbon-hydrogen bonds.

In order to enhance the reliability of the obtained results and confirm the stability of the system, we combined the docking methods with MD simulation. This allowed us to investigate the changes in the ligand-receptor complex’s conformation during the simulation time (Caporuscio et al., 2011). Simulations offer a detailed investigation of the exact movement of each atom over time, allowing us to investigate changes and fluctuations in protein patterns. MD simulation was conducted on the best-scoring L ₁₂ complex for 100 ns. The dynamic stability of the complexes was assessed by measuring the RMSD of the complex backbone atoms during the whole trajectory to confirm their stability (Figure 7). The RMSD curves of the ligand (Lig fit Prot) and the protein backbone (Cα) showed that our simulation had converged and equilibrated after an initial period of disturbance. Throughout the 100 ns simulation, variations in the RMSD of L ₁₂ in the binding pocket ranged from 1.70 to 4.20 Å. The MAPK^ERK altered its maximum RMSD from 0.85 Å in the starting frame to 3.20 Å. Higher fluctuations were observed between 22 and 30 ns, after which the RMSD remained constant until the last 100 ns. The convergence of RMSD values indicated that L ₁₅ and MAPK^ERK maintained their interaction through the simulation.

Root Mean Square Fluctuation (RMSF) serves as a metric quantifying the average deviation of each atom’s position from its mean position within a specified simulation or ensemble of structures. The RMSF values for the catalytic region of the MAPK^ERK of all the complexes exhibited stability. There were no considerable fluctuations observed where the ligand binds to the protein.The comprehensive analysis of the interaction between L ₁₂ and the binding site residues of MAPK^ERK is illustrated in Figure 8. The residues that interact with L ₁₂ are as follows: Ile31,Glu33,Ala35,Val39,Ala52,Lys54,Ile84,Gln105,Leu107,Leu112,Lys114,Lys151,Ser153,Asn154,Leu156,Cys166,Asp167. The residue interacting with L ₁₅ is shown in green, while the protein’s secondary structures, helices and β-strands, are represented by orange and blue bands, respectively. The RMSF values for the residues in the binding site were found to be less than 2 Å.

Ligand-protein interactions can be monitored throughout the simulation. These interactions can be categorized as hydrogen bonds, hydrophobic, ionic, and water bridges, and summarized. As depicted in Figure 9, Ala52, Ile84, and Leu156 engaged in hydrophobic interactions with the ligand for approximately 24%, 4%, and 7% of the simulation duration. Moreover, over at least 3% of the simulation period, Lys54, Gln105, and Asp111, Lys151, Asn154, and Asp167 formed water bridge interactions with the ligand. Significantly, Ser153 demonstrated a variety of interactions, encompassing water-bridged and hydrogen-bond interactions with the ligand. Consequently, this residue experienced numerous interactions throughout the simulation.

Throughout the 100-ns MD simulation, ligand properties were examined. The RMSD of a ligand with respect to its reference conformation (time t = 0) was calculated to be between 0.56 and 0.85 Å. The radius of gyration (Rg) provided a means to evaluate the variation in compactness of the ligand-protein complex within the range of 2.97–3.14 Å. The molecular surface area (MolSA) was computed using the van der Waals surface area between 250 and 260 Å. No intramolecular hydrogen bond (intra-HB) was detected. The solvent accessible surface area (SASA) analysis examines the interaction between a ligand and solvents during a 100 ns MD simulation, yielding values ranging from 40 to 158 Å. The analysis also included the examination of the polar surface area (PSA), which provides information about the solvent’s ability to access the surface area given by oxygen and nitrogen atoms. The PSA value ranged from 29 to 40 Å (Figure 10).

The physicochemical, pharmacokinetic, and medicinal chemistry features mentioned for our top 18 candidates are shown in Table 4. The above 18 lead ligands were found to have desirable drug-like properties based on Lipinski’s rule of five. The pharmacokinetic parameters revealed that all compounds are highly absorbed after oral administration through the gastrointestinal tract (GI). Structural alarms and pan-assay interference (PAINS) have been utilized in medicinal chemistry to forecast the presence of unstable, reactive, and toxic fragments in a compound’s structure (Baell and Holloway, 2010). All ligands in PAINS descriptors have zero alarms. The synthetic accessibility score (SA score) is a metric used to evaluate the ease of synthesizing drug-like molecules. It was observed that all the compounds possess a favorable SA score, indicating that they can be readily synthesized. The use of CaCo-2 cells, which are generated from human colon epithelial cells, is a widely accepted approach for studying the intestinal absorption of medicines in people. The CaCo-2 cell permeability findings showed that all ligands except L ₂ , L ₉ , L ₁₀ , L ₁₁ , L ₁₂ , L ₁₃ , L ₁₄ , L ₁₅ , L ₁₇ , and L ₁₈ were in the satisfactory range, indicating that these ligands have favorable membrane permeability. In terms of plasma glycoprotein (PGP) substrate, it was observed that, except for ligands L ₃ , L ₄ , L ₅ , and L ₁₁ , other ligands have inhibitory effects on PGP. As well, all ligands exhibit PGP substrate activity. The computed values for HIA indicate that all ligands except L ₁₂ , L ₁₄ , and L ₁₈ possess a high likelihood of being effectively absorbed through the intestinal membrane. The assessment of PPB has significance in the evaluation of drug safety. Drugs exhibiting a high PPB value (>90%) are associated with a narrow therapeutic index, whereas those with a low PPB value are considered to be comparatively safer. In the current investigation, it was shown that all ligands exhibited low PPB values. This finding suggests that these particular compounds have a wide therapeutic index, indicating a favorable safety profile. A substance with a positive blood-brain barrier value has a higher lipophilicity profile and may be easily absorbed from plasma membranes. It was noticed that ligand L ₁ had a higher lipophilicity profile, with values as high as BBB+++. Ligands L ₁₁ , L ₁₂ , L ₁₃ , L ₁₄ , and L ₁₆ are shown to be carcinogenic, while the other ligands are non-carcinogenic. The AMES toxicity profile showed that all ligands have the potential to be toxic. Overall, all ligands demonstrated better ADMET profiles; all values are shown in Table 5.