Virtual screening (VS) is a computational method to identify potential bioactive compounds from large chemical libraries using structure-based molecular docking (Shamsi et al., 2024b). In this study, structure-based VS was performed to evaluate the binding affinities of FDA-approved drugs against MEK1. The crystal structure of MEK1 (PDB ID: 7B9L) was retrieved from the Protein Data Bank (PDB) (Berman et al., 2014). Protein preparation included removal of crystallographic water molecules and co-crystallized ligands, protonation of side chains, addition of Kollman charges, and remodeling of missing residues. Missing residues were modeled using PyMod v3 (Janson and Paiardini, 2021) and processed in AutoDock Tools (Huey et al., 2012). A curated library of 3,500 FDA-approved drugs was sourced from the DrugBank database (Knox et al., 2024) in processed format. Ligand flexibility, including rotatable bonds and bond lengths, was fully permitted during docking. Blind docking was conducted using InstaDock (Mohammad et al., 2021) with a grid encompassing the entire protein structure (dimensions: 85 Å × 78 Å × 67 Å; center coordinates: −21.735, −5.746, 20.809 for the X, Y, and Z-axes, respectively). Docking poses were ranked by binding energy, and top-ranked conformations were selected for further analysis. For each compound, the top-ranked pose with the most favorable binding energy was extracted and analyzed further. Ligands showing the lowest (most negative) docking scores were prioritized for downstream evaluation.

The biological activities of the selected compounds were predicted using the PASS (Prediction of Activity Spectra for Substances) web server (Filimonov et al., 2014). This tool employs structure-activity relationship (SAR) models trained on diverse biological activity datasets to calculate the probability of a molecule exhibiting specific pharmacological effects, expressed as “probability to be active” (Pa) and “probability to be inactive” (Pi). Compounds with Pa > Pi were prioritized as high-confidence candidates for further investigation. However, PASS predictions do not confirm molecular target specificity and should be interpreted as preliminary, hypothesis-generating tools. The PASS results require subsequent validation through experimental assays or mechanistic modeling. To investigate binding site interactions, molecular docking poses of MEK1–ligand complexes were analyzed using visualization tools. Polar contacts, hydrogen bonds, and hydrophobic interactions were visualized using PyMOL (version 2.5.4) (DeLano, 2002). Discovery Studio Visualizer (v2023) (Visualizer, 2005) was employed to map interactions with residues critical to MEK1’s ATP-binding and catalytic sites for in-depth analysis of binding modes. Key binding residues responsible for ligand stabilization were identified based on distance thresholds and chemical compatibility (e.g., ≤3.5 Å for hydrogen bonds), facilitating insight into potential MEK1 inhibition mechanisms.

Molecular dynamics (MD) simulations were performed to study the dynamic behavior, conformational stability, and flexibility of MEK1 alone and in complex with Alectinib, Radotinib, and Selumetinib. All simulations were carried out with GROMACS 2022.4 (Van Der Spoel et al., 2005) with a CHARMM36m force field (Huang and MacKerell, 2013) using the TIP3P (Mark and Nilsson, 2001) water model. Ligand topology files were created with the CGenFF web server (Zhu, 2019), submitted for geometry optimization, and the charges checked with Avogadro 1.2.0 and our own Python scripts. Four different systems were set up: apo-MEK1, MEK1-Alectinib, MEK1-Radotinib, and MEK1-Selumetinib. All systems were minimized in energy via the steepest descent algorithm (5,000 steps) to remove steric clashes. During this phase, equilibration was conducted in two stages: (1) NVT ensemble (constant number of particles, volume, and temperature) ensemble for 100 ps at 300 K (Berendsen thermostat), and then (2) NPT ensemble (constant number of particles, pressure, and temperature) ensemble for 100 ps at 1 bar (Parrinello-Rahman barostat). Using periodic boundary conditions and a 2-fs time step, production runs were performed for 500 ns. Trajectory analysis performed using GROMACS utilities: the root mean square deviation (RMSD) and the root mean square fluctuation (RMSF) were calculated with gmx rms and gmx rmsf, respectively, to evaluate the structural stability and residue flexibility. The radius of gyration (Rg; gmx gyrate) and solvent-accessible surface area (SASA; gmx sasa) were used to calculate the compactness and solvent accessibility, respectively. We examined the dynamics of the hydrogen bond using gmx hbond. All plots were generated using the Grace (XMGRACE) plotting software. These analyses together provided a comprehensive picture of protein-ligand complex behavior, conformational adaptability, and structural stability throughout the simulation timeframe.

Principal component analysis (PCA) is a powerful dimensionality reduction and multivariate analysis technique widely used to extract meaningful large-scale motions from MD simulations of biomolecular systems. In the context of this study, PCA was employed to identify and quantify the dominant collective motions and structural fluctuations of MEK1 in its apo and ligand-bound states, thereby elucidating ligand-induced stabilization or flexibility. By filtering out high-frequency random noise, PCA allows the capture of functionally relevant, low-frequency conformational transitions that can be crucial for understanding binding mechanisms and structural stability. The primary reason for employing PCA in this analysis was to explore whether Radotinib and Alectinib modulate MEK1’s conformational dynamics more effectively than Selumetinib. Specifically, PCA helped characterize large-scale domain motions and determine whether these repurposed drugs stabilize MEK1 in a restricted conformational subspace, which may correlate with improved binding affinity and reduced entropic penalties (Papaleo et al., 2009). For PCA calculations, the covariance matrix of atomic positional fluctuations was generated using the Cα atoms from 500-ns MD trajectories, employing the gmx covar module in GROMACS. The covariance matrix captures correlated displacements of atom pairs and is defined as: C ij = <   x i   ‐   <   x i >   x j   ‐   <   x j >   > where x_i/x_j represents the coordinate of the i^th/j^th atom, and < - > is the ensemble average over the simulation time. Diagonalization of this matrix yields eigenvectors representing directions of motion and corresponding eigenvalues that quantify the variance along those directions. Singular value decomposition (SVD) was applied for matrix diagonalization. This approach enables the identification of ligand-induced effects on protein flexibility by comparing the conformational sampling across systems.

Gibbs free energy landscapes (FELs) were constructed to quantify the thermodynamic stability of MEK1 conformations in apo and ligand-bound states. Projections of the MD trajectories onto the first two principal components (PC1 and PC2), obtained via PCA, were used as reaction coordinates for FEL construction (Papaleo et al., 2009). The conformational probability distribution P was computed using the gmx sham module in GROMACS, and the corresponding free energy (ΔG) was estimated using the Boltzmann relation: Δ G = ‐ k B T ⁡ ln P where P is the probability density, k _B is the Boltzmann constant, and T is the simulation temperature. Energy basins (low ΔG) represent thermodynamically stable states, while peaks correspond to high-energy transition states. FELs were visualized as contour plots, with colors scaled from blue (low energy) to red (high energy). This approach elucidated ligand-specific stabilization patterns and identified dominant conformational clusters, revealing how Alectinib and Radotinib modulate MEK1’s energy landscape compared to Selumetinib.

Molecular docking serves as a predictive tool to model the interaction between a ligand and its target protein, allowing for identifying molecules with high binding affinity and appropriate conformational orientation within the binding pocket (Muhammed and Aki-Yalcin, 2024). In this study, a virtual screening of 3,500 FDA-approved drugs, curated from the DrugBank repository, was performed against the crystal structure of MEK1 (PDB ID: 7B9L) using the InstaDock platform (Mohammad et al., 2021). The objective was to identify potential MEK1 inhibitors among existing therapeutics that can be repurposed for cancer treatment. Following the docking process, compounds were ranked based on their binding energies, and the top 10 candidates were selected for further analysis (Table 1). Docking scores ranged from −10.1 to −10.8 kcal/mol, indicative of strong binding potential. Notably, all candidates outperformed Selumetinib (−7.2 kcal/mol), underscoring their superior affinity. Importantly, all selected candidates demonstrated significantly better docking scores than the reference MEK1 inhibitor Selumetinib, which exhibited a binding energy of −7.2 kcal/mol. This difference in binding energy suggests that the shortlisted compounds may engage MEK1 more effectively than Selumetinib, potentially translating to improved inhibitory activity. The findings from this docking-based screening highlight the potential of drug repurposing strategies to uncover alternative MEK1 inhibitors that may overcome limitations associated with current therapeutics, such as resistance or toxicity. However, molecular docking represents only the initial filtering step. Further dynamic and energetic assessments, such as MD simulations and binding free energy calculations, are essential to validate these interactions under more physiologically relevant conditions.

Drug profiling and predicting other biological activities are key elements of any computational drug repurposing workflow. After the initial molecular docking, the top 10 candidate compounds were subjected to biological activity prediction using the PASS server. The PASS algorithm uses structure‐activity relationship (SAR) models to predict pharmacological effects (Filimonov et al., 2014). The drug profiling for the 10 docked molecules was carried out to evaluate their pharmacological potency in MEK1 inhibition and relevant anticancer properties (Supplementary Table S1). Specifically, Radotinib and Alectinib were identified to have appropriate drug profiles and higher potential for anticancer-related indications. Though other compounds exhibited strong docking potential, their lack of predicted anticancer activity excluded them from further study. The PASS analysis for the selected molecules was carried out to evaluate their pharmacological potency other than MEK1 inhibition (Table 2). Radotinib and Alectinib exhibited greater anticancer potential, with high Pa values across diverse anticancer categories. Importantly, Pa values for both compounds were substantially higher than their respective Pi values, indicative of biological activity. For instance, Radotinib had Pa values particularly high as a growth factor agonist (Pa = 0.797), protein kinase inhibitor (Pa = 0.790), and Bcr-Abl kinase inhibitor (Pa = 0.748), etc. (Zabriskie et al., 2015). Such profiles are consistent with its known mechanism as a tyrosine kinase inhibitor and compatibility with MEK1 targeting.

Both Radotinib and Alectinib had high Pa values (≫ Pi) for multiple anticancer activities (Table 2). For instance, Radotinib’s top predicted activities include growth factor agonism (Pa≈0.80) and kinase inhibition (Pa≈0.79), consistent with its tyrosine kinase profile. Alectinib likewise showed high Pa for antineoplastic activities. These profiles suggest promising anticancer potential. We further screened Alectinib, an FDA-approved drug for non-small cell lung cancer. We provided favorable probabilities for antineoplastic and antimetastatic activities, demonstrating the program’s utility in onco-logic settings. These results imply that Radotinib and Alectinib show a high binding affinity toward MEK1, coupled with their complementary biological activity profiles consistent with MEK1 inhibition and the treatment of their indicated cancers (Herden and Waller, 2018). Hence, integrating docking data with the biological activity prediction improves the confidence of these hits as potential drug candidates for repurposing. However, PASS predictions are only indicative and based on in silico SAR models. Although PASS has a high reported cross-validation accuracy (∼95%) high Pa values do not guarantee an actual biological effect. Thus, over-reliance on PASS is cautioned–experimental assays are needed to confirm the predicted activities.

Radotinib and Alectinib were prioritized for further analysis due to their strong docking scores and favorable biological activity profiles. A detailed interaction analysis was conducted for Radotinib, Alectinib, and the reference MEK1 inhibitor Selumetinib using PyMOL and Discovery Studio Visualizer to gain mechanistic insight into their binding modes. Protein–ligand interactions are fundamental to understanding ligand efficacy, as they influence enzymatic activity, signal transduction, protein stability, and drug specificity. Visual representation of the docked complexes revealed that all three compounds, Radotinib (cyan), Alectinib (magenta), and Selumetinib (orange), bound within the MEK1 allosteric pocket located adjacent to the ATP-binding site (Figure 1). Figure 1A shows a cartoon representation of MEK1 highlighting the binding locations of all three ligands. Zoomed-in views in Figures 1B–D illustrate the binding orientations of Alectinib, Radotinib, and Selumetinib, respectively. The ligands predominantly interacted with key regions involved in ATP and inhibitor binding, including the Leu74–Val82 cleft, Lys97 (ATP/inhibitor-binding site), Asp190 (proton acceptor active site), and the Asp208–Val211 segment, indicating their potential to modulate MEK1 activity through direct interaction with catalytically relevant residues.

The two-dimensional interaction maps provided a detailed account of the molecular interactions (Figure 2). Alectinib displayed a stable and meaningful interaction pattern. It formed a hydrogen bond with Arg189 and engaged in pi-pi stacking with His239. A pi–sulfur interaction was observed with Arg234, while hydrophobic alkyl and pi–alkyl interactions involved Arg189, Ile216, and Met219. Van der Waals interactions further anchored the ligand via contacts with residues such as Gly79, Lys97, Ile99, His188, Asp190, Asp208, Phe209, Gly210, Leu215, Asp217, Ala220, Asn221, Gly237, Thr238, and Tyr240 (Figure 2A). Notably, Alectinib interacted with Gly79, Lys97, and Asp190, which are part of the catalytic core, further suggesting its therapeutic relevance in MEK1 inhibition (Di Fruscia et al., 2021). At the same time, Radotinib also demonstrated a rich interaction profile, forming hydrogen bonds with residues Gly79, Lys97, Ser150, Ser194, Asn195, and Asp208. Additionally, it exhibited halogen interactions (e.g., fluorine) with Asp152, Gln153, and Ser194, along with pi–sulfur interactions involving Lys97, Asp190, Asp208, and Met143. Radotinib’s halogen bond with Asp152 may enhance binding specificity by mimicking ATP’s phosphate interactions. Hydrophobic contacts such as alkyl and pi–alkyl interactions were observed with Leu74, Val82, Ala95, Ile141, Leu197, and Cys207, while van der Waals forces stabilized the complex through interactions with Gly75, Ala76, Gly77, Asn78, Gly80, Ile99, Leu115, Met146, Lys192, Phe209, Leu215, Met219, and Asn221 (Figure 2B).

These extensive interactions, especially with residues in the ATP-binding cleft and active site, underscore Radotinib’s potential as a robust MEK1 inhibitor. In contrast, Selumetinib, the reference compound, is a non-ATP-competitive MEK inhibitor, showed comparatively limited interactions. It formed hydrogen bonds with Lys97 and Ser212 and engaged in pi–sigma and pi–sulfur interactions with Leu215 and Met143, respectively. Hydrophobic contacts included interactions with Leu118, Val127, Ile141, Ile216, and Leu215, while van der Waals interactions occurred with Leu115, Cys207, Asp208, Gly210, Val211, and Met219 (Figure 2C). Among the residues in the ATP-binding region, Selumetinib directly interacted with Lys97, indicating a more constrained binding profile compared to Radotinib and Alectinib. Both Radotinib and Alectinib demonstrated broader and more diverse interaction networks within the MEK1 binding pocket than the reference inhibitor. Their ability to engage multiple key residues, particularly at the ATP-binding and allosteric sites, reinforces their potential as strong MEK1 inhibitors. These results provided a structural rationale for their selection and prompted further assessment using MD simulations to evaluate binding stability and dynamic behavior under physiological conditions.

MD simulations have become an indispensable tool in molecular biology and drug discovery, enabling high-resolution insights into the atomic-level behavior of proteins and their interactions with ligands (Singh and Singh, 2020). In this study, MD simulations were employed to investigate the dynamic behavior of MEK1 in its apo form and in complex with three compounds: Radotinib, Alectinib, and Selumetinib. The simulations were conducted using GROMACS 2022.4 on a Linux platform, with a total simulation time of 500 nanoseconds for each system. Trajectory analysis of these simulations provided multiple parameters to assess the dynamic behavior and stability of MEK1 and its complexes. Three energy components were initially evaluated: potential, kinetic, and total energy. The potential energies (in kJ/mol) for the MEK1, MEK1-Alectinib, MEK1-Radotinib, and MEK1-Selumetinib systems were found to be −1,058,300, −742,093, −743,103, and −742,440, respectively. The corresponding kinetic energy values were 202,973, 145,759, 145,871, and 145,765 kJ/mol, while the total energy values were −855,325, −596,335, −597,231, and −596,675 kJ/mol, respectively. These energy profiles suggest that all three complexes exhibit stable energy states, with MEK1 retaining its structural integrity even upon drug binding. Notably, the drug-bound systems showed consistently lower total and potential energy values than the apo form, indicating enhanced thermodynamic stability upon ligand binding.

To illustrate the most dominant motion of the MEK1 protein after binding of Alectinib, Radotinib, and Selumetinib drugs during MD simulation, PCA was utilized on the coordinate covariance matrix extracted from trajectories of a 500 ns MD simulation. Most of the structural motions are captured through the first two principal components; thus, we performed PCA analysis of MEK1, MEK1-Alectinib, MEK1-Radotinib, and MEK1-Selumetinib complexes. Figure 8 displays the superimposed PCA plot for MEK1, MEK1-Alectinib, MEK1-Radotinib, and MEK1-Selumetinib complexes of the first two principal components. As the PCA plot shows, each complex occupied a different range of vibrational space, which indicates different motion patterns. The MEK1-Alectinib and MEK1-Radotinib complexes covered less vibrational space, while the reference MEK1-Selumetinib complex was dispersed and occupied a larger area. The area of motion covered by MEK1 at PC1 −4 nm to 12.9 nm at PC2 −5.2 nm to 7.0 nm, MEK1-Alectinib complex at PC1 −9.1 nm to 10.3 nm at PC2 −8.5 nm to 7.4 nm, MEK1-Radotinib complex at PC1 −5.3 to 7.2 at PC2 5.5 nm–8.0 nm and MEK1-Selumetinib complex at PC1 −8.7 nm to 14.7 nm at PC2 15.5 nm–10.9 nm (Figure 8A).

The results showed that the first PC1 carried most of the motion. For example, PC1 had an eigenvalue of ∼15.2 (accounting for roughly 71% of the total variance), and the second component (PC2) had an eigenvalue of ∼3.4 (∼16% of variance). Together, PC1 and PC2 explained 87% of the total conformational fluctuation (dominant motions). These values indicate that nearly all collective motion is captured by the first two modes, which is consistent with typical MD PCA results where PC1 dominates. The calculated values and superimposed PCA plot show that MEK1-Alectinib and MEK1-Radotinib complex display a cluster and compact type of motion compared to MEK1-Selumetinib complex. The time‐dependent eigenvector traces similarly show that projections onto PC1 and PC2 remain relatively stable for the Alectinib and Radotinib cases, whereas the Selumetinib complex exhibits larger fluctuations (Figure 8B). Together, the high percentage of variance in PC1+PC2 and the compact clusters for Alectinib/Radotinib indicate that most conformational variability is confined to a few dominant motions.

Further to describe the structural dynamics of MEK1, MEK1-Alectinib, MEK1-Radotinib, and MEK1-Selumetinib complexes, the metastable conformations state, which were dominant during the simulation, were fetched from the FEL. Different energy states are indicated by different colors, from highest to lowest, denoted by red to blue. The FEL of MEK1, MEK1-Alectinib, MEK1-Radotinib, and MEK1-Selumetinib complexes are given in Figure 9. The MEK1 FEL map had one long blue energy minimum and multiple energy funnels linked (Figure 9A). The MEK1-Alectinib complex map shows a broader area of blue basin with combined energy funnels (Figure 9B). In contrast, the MEK1-Radotinib complex had two blue energy minima and two separate energy funnels, which indicate two different folding states (Figure 9C). The reference MEK1-Selumetinib complex map shows multiple energy funnels, which are separated from each other, demonstrating different folding states (Figure 9D). The resulting PCA analysis and FEL calculation suggested that MEK1-Alectinib and MEK1-Radotinib complexes were in stable form during simulations. These observations imply that Radotinib and Alectinib modulate the conformational flexibility of MEK1 more effectively than Selumetinib. Their distinct energy basins and stable low-energy states suggest a reduced likelihood of unfavorable conformational transitions. This thermodynamic stability further supports their potential as robust MEK1 inhibitors. Collectively, the FEL and PCA analyses reinforce the structural reliability and inhibitory promise of these repurposed compounds.

Importantly, both Radotinib and Alectinib are approved kinase inhibitors with known off-target profiles. Radotinib is a BCR-ABL1 tyrosine kinase inhibitor that also inhibits DDR, EPHB, LYN, and PDGFR kinases at low-nanomolar levels (Zabriskie et al., 2015; Liu et al., 2022). Alectinib is an ALK inhibitor with potent activity against RET kinase (IC50 ≈ 1.9 nM for ALK, 4.8 nM for RET (Kodama et al., 2014). These multi-kinase activities mean that repurposing them for MEK1 could bring unintended effects via their original targets. Therefore, future work should include broad kinase profiling to confirm MEK1 selectivity and assess off-target risks. A key limitation of this study is its purely in silico nature. All findings must be validated experimentally; in future work, we plan in vitro kinase assays, cell-based MEK1 activity tests, and in vivo studies to confirm efficacy. Additionally, biochemical assays will be used to determine the selectivity profiles of Radotinib and Alectinib against a panel of kinases.