Exploring the Natural Compounds in Flavonoids for Their Potential Inhibition of Cancer Therapeutic Target MEK1 Using Computational Methods

The Mitogen-Activated Protein Kinase (MAPK) signaling pathway plays an important role in cancer cell proliferation and survival. MAPKs’ protein kinases MEK1/2 serve as important targets in drug designing against cancer. The natural compounds’ flavonoids are known for their anticancer activity. This study aims to explore flavonoids for their inhibition ability, targeting MEK1 using virtual screening, molecular docking, ADMET prediction, and molecular dynamics (MD) simulations. Flavonoids (n = 1289) were virtually screened using molecular docking and have revealed possible inhibitors of MEK1. The top five scoring flavonoids based on binding affinity (highest score for MEK1 is −10.8 kcal/mol) have been selected for further protein–ligand interaction analysis. Lipinski’s rule (drug-likeness) and absorption, distribution, metabolism, excretion, and toxicity predictions were followed to find a good balance of potency. The selected flavonoids of MEK1 have been refined with 30 (ns) molecular dynamics (MD) simulation. The five selected flavonoids are strongly suggested to be promising potent inhibitors for drug development as anticancer therapeutics of the therapeutic target MEK1.


Introduction
Cancer is a leading cause of death in nearly every country in the world, as well as being the most significant barrier to extending life expectancy in the twenty-first century [1]. More than nearly 20 million new cases are predicted to be registered by 2025. Various implications are considered as hallmarks of cancer cells, including proliferative signaling, growth suppressor escape, cell death resistance, immortality, and angiogenesis induced by invasion-metastasis [2]. Since cancer is generally associated with several mutations that affect the main signaling pathways [3], targeted cancer therapy takes advantage of tumor-specific genetic vulnerabilities and mutations and designs therapeutics directly targeting cancer cells. The important therapeutic targets (mostly enzymes) are those whose suppression eventually kills cancer cells and that have minimal effect on normal tissues as they are either mutated or over-expressed in cancer cells [4,5]. One of the pathways that is commonly altered and activated, having a role in oncogenesis, the progression of the tumor, and drug resistance in cancer, is the Mitogen-Activated Protein Kinases (MAPK) signaling pathway [3,6]. The MAPK cascade, which is also known as RAS-regulated RAF-MEK1/2-ERK1/2 or ERK signaling pathway [7], is comprised of a number of kinases that carry out specific cell fate decisions in response to the input signal, explaining why numerous cer types [45][46][47][48][49][50][51]. Myricetin is one of the flavonoids which has been found to inhibit MEK1 in vivo and in vitro, in addition to being a potent ATP-noncompetitive inhibitor of MEK1 [52,53]. Another flavonoid is nobiletin, which has shown an in vitro antitumor effect against MEK in human fibrosarcoma HT−1080 cells [54]. Isorhamnetin is a flavonoid that inhibits MEK1 in the ATP-noncompetitive binding site [55]. A flavonoid that has shown potent and specific ATP-noncompetitive inhibition activity of MEK1 is PD 098059, which is the first MEK inhibitor to be discovered [56] and mostly used in laboratory experiments.
The current study proposed five novel MEK1 inhibitors and anticancer drug candidates by virtual screening of the natural compound flavonoids. A library of all known flavonoids was prepared and the top five screened compounds were proposed as potential MEK1 inhibitors and anticancer agents. The proposed compounds were further checked for key interacting residues, molecular interactions, binding energy, and dissociation constant using various methods. Finally, to show the stability of the protein-ligand, complexes were subjected to MD simulation analysis. This study will provide novel MEK1 inhibitors and anticancer agents with their action mechanism. The proposed flavonoids can further be tested experimentally for their potential use as novel anticancer agents.

Virtual Screening of Natural Compound Class Flavonoids for MEK1 Inhibition
The capacity of molecular docking studies to anticipate the right binding conformations of small molecules as ligands to the appropriate target binding site is an important feature of in silico drug discovery. In this view, the goal of our study was to do virtual screening of flavonoids for possible MEK1 inhibition, which is an essential protein target in the MAPK pathway in cancer cells. (Figure 1) shows the histogram of dock scores for the screened compounds (Supplementary File S1) where the selected top five ranked compounds are highlighted in the left. Most of the compounds are providing docking affinity (>−7.0 kcal/mol), but the highest scoring and best fit compounds were selected and the redundant compounds were omitted to give the best and most diverse structure compounds. The top five selected flavonoids ( Figure 2) against MEK1 in (Table 1) show the highest docking score of (−10.8 kcal/mol) and the fifth rank with a score of (−10.4 kcal/mol) compared to the docking score of native inhibitor (−9.0 kcal/mol). The native inhibitor of the 3D structure is PD318008 (5-bromo-N-(2,3-dihydroxypropoxy)-3,4dofluoro−2-[(2-fluoro-4-iodophenyl)amino]benzamide), an analog of PD184352 which is a highly selective, ATP-noncompetitive, and potent MEK1 and MEK2 inhibitor [16]. These results suggest that the selected flavonoids are probably better inhibitors of MEK1. Further, the molecular docking of trametinib, a specific and highly selective inhibitor of MEK1/2, was performed and gave the score of (−9.7 kcal/mol), which is also lower than those of the selected flavonoids, providing credence to the selected compounds [69].    better the binding.    as required for adequate inhibition ( Table 1). The molecular interaction (Table 2) showed 2 hydrogen bonds and 62 (16) non-bonded contacts (hydrophobic interactions). The key interacting residues were Leu-118, Ile-141, Asp-208, Phe-209, and Ile-216 with 5, 6, 5, 9, and 9 non-bond contacts, respectively. The Met-219 turned out to have the maximum ∆ASA (loss in Accessible Surface Area) (42.37 Å2) followed by Asp-208 (39.18 Å2). Hydrogen bonds were formed by Val-127 and Gly-128 measure 2.42 Å and 3.21 Å, respectively. Observing amino acid residues from the flavonoid binding site within the target protein seeks to predict the interactions that occur and that are thought to contribute to the flavonoid compound's pharmacological activities, such as MEK kinase inhibition. Inhibition of enzymatic activity of a protein by a compound is largely due to non-covalent bonds, including nonbonded contacts and hydrogen bonds [70]. The binding of native inhibitor with interacting residues and their interactions is also provided for comparison ( Figure 4A). Comparing to the binding of the native inhibitor, there are six residues, Leu-118, Ile-141, Asp-208, Phe-209, Leu-215, and Met-219, common to the lists of interacting residues of this compound and the native inhibitor, and they include the key residues Asp-208 and Phe-209. The fact that the proposed compound binds to the similar group of residues in the catalytic site validates our prediction that it is comparable to the native inhibitor.  Figure 4A), including the key residues Lys-97, Asp-208, and Phe-209. The proposed compound bound to the same set of residues in the catalytic site, making it promising as the native inhibitor.    Figure 4C). The strength of binding scores' docking affinity of (−10.6 Kcal/mol), binding energy (−10.96 Kcal/mol) and dissociation constant (pK d , 8.03) were reasonably high, as required for good inhibition ( Table 1). The protein-ligand complex was stabilized by non-bonding interactions (Table 3) through 78 (20) non-bonded contacts and 2 hydrogen bonds. The key interacting residues were Ile-99, Asp-208, Lys-97, and Phe-209 with 15,9,8, and 6 nonbonding interactions of each, respectively. Met-219 turned out to have the maximum ∆ASA with (43.97 Å 2 ), followed by Asp-208 with (39.24 Å 2 ) and Lys-97 with (31.09 Å 2 ). Hydrogen bonds in Gly-80 and Lys-97 measured 2.80 Å and 2.99 Å, respectively. Of the 19 interacting residues, 7 residues were common with those of the binding of the native inhibitor: Lys-97, Leu-118, Ile-141, Asp-208, Phe-209, Leu-215, and Met-219 ( Figure 4A), including the key residues Lys-97, Asp-208, and Phe-209. The proposed compound bound to the same set of residues in the catalytic site, making it promising as the native inhibitor.   Figure 4D). The predicted binding strength scores for the protein-ligand complex were docking affinity of (-10.5 Kcal/mol), binding energy (−9.49 Kcal/mol), and dissociation constant ( , 6.95) and showed quality binding as required for good inhibition ( Table 1). The molecular interaction showed one hydrogen bond through Ser-212 measures 3.17 Å and 54 (12) non-bonded contacts (hydrophobic interactions). The key interacting residues (Table 4) are Asp-208 The residues forming non-bonding interactions are shown as red bristles, while residues forming hydrogen bond and the bound ligand are shown as ball-and-stick representations. The carbon atoms are shown as black balls, nitrogen atoms as blue balls, oxygen atoms as red balls, fluorine atoms as green balls, bromine atom as pink balls, and iodine atom as a yellow ball. The interacting residues common with those of the native inhibitor are shown in circles. The hydrogen bonds are shown as green dashed lines labeled with bond length (in Å). The third rank flavonoid (CID: 10991656) bound to the MEK1 binding pocket ( Figure 3) and showed interactions with 12 residues: Leu-115, Leu-118, Ile-141, Met-143, Asp-190, Asp-208, Phe-209, Val-211, Ser-212, Leu-215 Ile-216, and Met-219 ( Figure 4D). The predicted binding strength scores for the protein-ligand complex were docking affinity of (−10.5 Kcal/mol), binding energy (−9.49 Kcal/mol), and dissociation constant (pK d , 6.95) and showed quality binding as required for good inhibition ( Table 1). The molecular interaction showed one hydrogen bond through Ser-212 measures 3.17 Å and 54 (12) nonbonded contacts (hydrophobic interactions). The key interacting residues (Table 4) are Asp-208 and Phe-209, with 13 and 10 non-bonding interactions of each, respectively. The maximum ∆ASA is with Met-219 residue (43.97 Å 2 ) followed by Asp-208 with (40.1 Å 2 ). The interacting residues and their interactions with the native inhibitor are also provided for comparison ( Figure 4A). Out of the 12 interacting residues, 7 residues were common with those of the native inhibitor: Leu-118, Ile-141, Met-143, Asp-208, Phe-209, Leu-215, and Met-219, which included the key residues Asp-208 and Phe-209. This also suggested that the proposed compound was blocking the same set of residues as the native inhibitor, thereby inhibiting the protein's action.   Table 1). The key interacting residues are Asp-208 and Phe-209 with 11 and 14 non-bonding interactions of each, respectively. Met-219 turned out to have the maximum ∆ASA with (52.28 Å 2 ), then Asp-208 with (38.96 Å 2 ) and next Asp-190 with (30.23 Å 2 ) ( Table 5). Compared to the native inhibitor ( Figure 4A), there are seven common amino acids: Leu-118, Ile-141, Met-143, Asp-208, Phe-209, Leu-215, and Met-219, and they share the same key residue of Asp-208 and Phe-209 and thus, they may inhibit MEK1 kinase activity in the same way as the native inhibitor does. This compound is a chiral compound and possesses R and S stereoisomeric configurations, as shown in Supplementary Figure Table 1). The molecular interaction shows 3 hydrogen bonds and 44 (11) non-bonded contacts (hydrophobic interactions). The three hydrogen bonds with Asp-190, Arg-234, and ATP measure 2.73 Å, 3.04 Å, and 3.10 Å, respectively. The key interacting residues are Asp-208 and Phe-209, with 9 and 11 non-bonding interactions, respectively. Met-219 turned out to have the maximum ∆ASA with (50.93 Å2), and then Asp-208 with (40.01 Å2) ( Table 6). Of the 11 interacting residues, 6 residues were common among the interacting residues of the native inhibitor: Leu-118, Ile-141, Met-143, Asp-208, Phe-209, and Met-219 ( Figure 4A), including the same key residues Asp-208 and Phe-209; thus, they might inhibit the MEK1 kinase activity similar to the native inhibitor.

Drug-Likeness and Pharmacokinetics Prediction
The strength of the ligand binding on the target protein is not the only factor in the discovery of novel medications. To assess the degree of effectiveness and therapeutic efficacy, it is also studied in terms of drug-likeness, pharmacokinetics, and toxicity. Lipinski's rule of five predicts the drug-likeness of the selected compounds. Moreover, pharmacokinetics including Absorption, Distribution, Metabolism, Excretion, and Toxicity (ADMET) play an important role in medicinal chemistry, which describes how drugs move through the body. The prediction of drug-likeness for the selected inhibitor for MEK1 is presented in (Table 7). Considering the desired values: molecular weight < 500, H-bond donors < 5, H-bond acceptors < 10, Rotatable bonds < 10, and lipophilicity (logP) < 5, all the selected compounds follow the desired values, except for (CID: 10813589) and (CID: 10524567), where the lipophilicity is slightly higher than 5. All the ADMET predicted properties of the top selected inhibitors for MEK1 are presented in (Table 8). The efficacy of the selected compounds as oral medicine was determined using two models for measuring absorption properties, including CaCO2 permeability and intestinal absorption. Where the desired CaCO2 permeability is >0.90 and intestinal absorption >30%, the prediction of the screened compounds shows the CaCO2 permeability values all in positive integers, with exception of (CID: 10575055). While intestinal absorption shows a high percentage, all are higher than 70%, which is considered as good absorption. Skin permeability is the next important factor in absorption. The ideal skin permeability is >−2.5 log Kp and all of the compounds under study have permeability values of less than −2.5 log Kp, indicating poor skin permeability. The ATP-binding cassette (ABC) transporter, which is necessary for efficient molecular transport across cell membranes, contains P-glycoprotein. P-glycoprotein substrates, and inhibitors of P-glycoprotein I and II which were examined in all of the substances that were screened. Except for (CID: 10813589) and (CID: 10524567), all of the compounds were found to be substrates, indicating that they can be transported through the cell membrane through the ABC transporter. The (CID: 10575055) was found to be ineffective as an inhibitor of P-glycoprotein I transporter, suggesting that they could be incapable of inhibiting these drug efflux pumps. The distribution of the substances in the body was determined using four distinct assays: volume of distribution (VDss), fraction unbound, BBB permeability, and central nervous system (CNS) permeability. To begin with, in the VDss assay, which is used to evaluate the total amount of drugs needed for uniform drug distribution in the bloodstream, readings less than −0.15 log are considered negative, while values greater than 0.45 log are considered good diffusion. Thus, (CID: 10991656) and (CID: 10524567) show average VDss values, while other compounds have low distribution volume. The potential of a drug to reach the brain is determined by the permeability of the blood-brain barrier (BBB). If the logBB values are more than 0.3, they will cross BBB. The logBB value of the screened compounds are less than 0.3, meaning that none of them will be able to cross BBB except for (CID: 10991656). The desired value of the CNS permeability is >−2 and the screen shows good results except for (CID: 10575055), where it was less than the desired value. Seven distinct cytochrome models were used to examine the test drug's metabolism in the body. All of the compounds were tested for their capacity to inhibit CYP1A2, CYP2C19, CYP2D6, CYP2C9, and CYP3A4 as well as their ability to function as a substrate for CYP2D6 and CYP3A4. The total clearance rates of all of the examined compounds were varied, and none of them appeared to be a substrate for organic cation transporter 2 (OCT2). They also failed to anticipate AMES toxicity, showing that these chemicals are neither carcinogenic nor mutagenic, except for: (CID: 129696793), and (CID: 10813589). Three of the selected flavonoids predicted to be negative for hepatotoxicity were (CID: 10813589), (CID: 10991656), and (CID: 10524567), whereas none of the substances tested positive for skin sensitization. Overall, the selected compounds proposed to be safe drug-candidates for human cancer therapy. Table 7. Structures and chemical properties of the selected flavonoids against MEK1 to predict the drug-likeness: molecular weight, lipophilicity (LogP), number of (#) rotatable bonds, hydrogen acceptors and hydrogen donors, and polar surface area. positive for skin sensitization. Overall, the selected compounds proposed to be safe drugcandidates for human cancer therapy. positive for skin sensitization. Overall, the selected compounds proposed to be safe drugcandidates for human cancer therapy. positive for skin sensitization. Overall, the selected compounds proposed to be safe drugcandidates for human cancer therapy. positive for skin sensitization. Overall, the selected compounds proposed to be safe drugcandidates for human cancer therapy. positive for skin sensitization. Overall, the selected compounds proposed to be safe drugcandidates for human cancer therapy. Table 7. Structures and chemical properties of the selected flavonoids against MEK1 to predict the drug-likeness: molecular weight, lipophilicity (LogP), number of (#) rotatable bonds, hydrogen acceptors and hydrogen donors, and polar surface area.

MD Simulation
Molecular Dynamic (MD) simulation was performed to refine and assess the stability of protein after adding the missing loop (Supplementary Figure S3) and the binding stability of the protein-ligand complex system. The MD simulation evaluates and delineates the dynamic behavior of the ligand and the binding site residues. The best conformation flavonoids obtained from virtual screening for MEK1 inhibition advanced to MD simulation. The MD results were examined for RMSD of backbone, RMSD of heavy ligand atoms, RMSF values, hydrogen bonds number, and the radius of gyration to assess the stability of the protein-ligand complex. The RMSD value is a measure of how far a protein molecule deviates from its initial conformation over the course of the simulation. The RMSD values for the first rank were found to be within 0.2 nm for the backbone and 0.05-0.1 nm through the simulation, and the protein-ligand complex was considered to be stable ( Figure 5A). The second rank RMSD value of backbone began at 0.2 nm then rose to 0.3 nm in 20 ns, then went back to stabilize at 0.2 nm. For the heavy atoms of ligand RMSD ( Figure 5B computed to define the flexible areas of the protein during the course of the simulation. The RMSF is a metric that measures how much ligand binding affects the flexibility of MEK1 residues. The RMSF values ( Figure 5C) show that the maximum fluctuation was in the amino acid residue region 276-320, which is a proline rich loop and highly flexible region, and the fluctuation was within the range of 1.25 nm. From this perspective also, owing to low fluctuation in RMSF, the protein-ligand complex seems stable ( Figure 5C). The radius of gyration is a measure of protein compactness, and the low fluctuations in its value indicate stability of the protein backbone. The radius of gyration fluctuates in a range of around 2.0 nm for all the four simulations and during the entire simulation, which points towards stability of the protein ( Figure 5D). The presence of hydrogen bonds is critical for a protein complex's stability. Analysis of the number of hydrogen bonds in (Figure 6) shows them appearing and disappearing during the course of simulation. Looking more closely at (Figure 6), the maximum number of hydrogen bonds and pairs reaches eight for first rank compound, five for the second, two for the third and fourth ranks, and eight for the fifth rank compound. Due to the isolated hydrogen bonds and low average hydrogen-bond number per time frame, the hydrogen-bond network in the complexes appeared to be weak. Other interactions were thought to hold hydrogen bonds in places where they had disappeared. As a result, no notable conformational changes in the complexes were observed across the simulated time period. These findings suggested that the MD simulation trajectory for the complex after equilibrium was reliable enough for future investigation.

Data Retrieval and Preparation
The three-dimensional structure coordinate of MEK1 complexed with the native inhibitor is retrieved from the protein data bank (PDB, https://www.rcsb.org/, accessed on 19 January 2022). The structure of MEK1 (PDB code: 1S9J) with 2.40 Å resolution [16] was selected and used for the study. MODLLER was used to model the missing prolin-rich loop in the crystallized structure (residues from 276 to 305), which is a high flexible region located after the catalytic site [71]. The unique allosteric site for MEK1 was verified [72], and the coordinates for the grid box covering the catalytic site were prepared using Au-toDockTools−1.5.6 [73]. Other preparations included: deleting water, checking for missing atoms, removing heteroatoms, adding polar hydrogens, computing and adding charges, and finally converting the protein into a (pdbqt) file for the molecular docking, also performed using AutoDockTools−1.5.6 software package. A drug-like library prepared from PubChem (www.ncbi.nlm.nih.gov, accessed on 4 April 2021) and 2630 flavonoids were filtered to 1289 by 3D structure availability and the Lipinski rule of five [74]. Ligands were prepared for virtual screening using Open Babel command line [75] and converted from (sdf) file to (pdbqt) after adding charges and hydrogens. The 2-D illustrations for the chemical compounds were prepared using MarvinSketch v18.4, ChemAxon (http://www.chemaxon.com/products/marvin, accessed on 6 June 2021).

Molecular Docking
Docking was carried out using AutoDock Vina [76] after preparing the configuration file with the details of the grid box coordinates, with energy range of 4 and maximum

Data Retrieval and Preparation
The three-dimensional structure coordinate of MEK1 complexed with the native inhibitor is retrieved from the protein data bank (PDB, https://www.rcsb.org/, accessed on 19 January 2022). The structure of MEK1 (PDB code: 1S9J) with 2.40 Å resolution [16] was selected and used for the study. MODLLER was used to model the missing prolinrich loop in the crystallized structure (residues from 276 to 305), which is a high flexible region located after the catalytic site [71]. The unique allosteric site for MEK1 was verified [72], and the coordinates for the grid box covering the catalytic site were prepared using AutoDockTools-1.5.6 [73]. Other preparations included: deleting water, checking for missing atoms, removing heteroatoms, adding polar hydrogens, computing and adding charges, and finally converting the protein into a (pdbqt) file for the molecular docking, also performed using AutoDockTools-1.5.6 software package. A drug-like library prepared from PubChem (www.ncbi.nlm.nih.gov, accessed on 4 April 2021) and 2630 flavonoids were filtered to 1289 by 3D structure availability and the Lipinski rule of five [74]. Ligands were prepared for virtual screening using Open Babel command line [75] and converted from (sdf) file to (pdbqt) after adding charges and hydrogens. The 2-D illustrations for the chemical compounds were prepared using MarvinSketch v18.4, ChemAxon (http://www.chemaxon.com/products/marvin, accessed on 6 June 2021).

Molecular Docking
Docking was carried out using AutoDock Vina [76] after preparing the configuration file with the details of the grid box coordinates, with energy range of 4 and maximum exhaustiveness of 24. Best mode, least RMSD, and highest docking affinity results were taken and ranked for MEK1. For further analysis of the docking, protein-ligand interaction plots of selected flavonoids with MEK1 was performed using Ligplot+ v1.4.5 [77] and illustrations of the docking were prepared using PyMOL v2.4.0 [78]. Further, calculations of binding energy and dissociation constant were performed by XScore v1.2.11 [79]. The degree of ligand filling the binding site was evaluated by loss in accessible surface area (ASA).
∆ASA i = ASA protien i − ASA protien−ligand i A residue is said to be taking part in filling the binding site if it loses more than 10 A 2 ASA due to binding [80]. All the ASA calculations of the protein-ligand complexes and the unbound proteins were performed by Naccess v2.1.1 [81].

Molecular Dynamic Simulation
The docked protein-ligand complexes were subjected to energy minimization using Gromacs v2020.5 [83] with the CHARMM36 all atom force field. Ligand and protein were separated to add ligand hydrogen atoms using Avogadro v2018 (https://avogadro.cc/, accessed on 6 January 2022) and then converted for topology using CHARMM force field (https://cgenff.umaryland.edu/, accessed on 6 January 2022), then wrote back with the complex topology file. The models were solvated with a water model in a cubic periodic box with 1 nm distance from the edge of the complex atoms. The solvated system was neutralized by five sodium ions. Energy minimization was carried out through 50,000 steps. An equilibration was conducted by number of particles, volume, and temperature (NVT), and number of particles, pressure, and temperature (NPT) temperature was coupled for ligand, protein, solvent, and ions, separately. Then the system proceeded to the actual MD simulation. The final models obtained at the end of MD were validated and illustrated by VMD (https://www.ks.uiuc.edu/Research/vmd/, accessed on 5 September 2021). For the analysis, Gromacs and Xmgrace (https://plasma-gate.weizmann.ac.il/Grace/, accessed on 7 September 2021) were used.

Conclusions
In conclusion, the molecular docking results suggest that some flavonoids could be better inhibitors of MEK1 compared to the native inhibitor based on the binding affinity and ligand interactions. The selected flavonoids could be potential drug candidates after re-engineering to improve the pharmacokinetic properties. Further, MD simulation studies with 100 ns time scale confirm the stability of the first rank flavonoid and MEK1 complex by root-mean-square deviation, root-mean-square fluctuation, and the radius of gyration. Our findings suggest that natural flavonoids are a promising and readily available source of anticancer targeted therapy in the future. However, these interpretations need further confirmatory analysis and validations for the screened molecules to ascertain their efficacy in the illness treatment.