In silico drug screen reveals potential competitive MTHFR inhibitors for clinical repurposing

Abstract MTHFR (Methylenetetrahydrofolate reductase) is a pivotal enzyme involved in one-carbon metabolism, which is critical for the proliferation of cancer cells. In line with this, published literature showed that MTHFR knockdown caused impaired growth of multiple types of cancer cells. Moreover, higher MTHFR expression levels were linked to shorter overall survival in hepatocellular carcinoma, adrenocortical carcinoma, and low-grade glioma, bringing the need to design MTHFR inhibitors as a possible treatment option. No competitive inhibitors of MTHFR have been reported as of today. This study aimed to identify potential competitive MTHFR inhibitor candidates using an in silico drug screen. A total of 30470 molecules containing biogenic compounds, FDA-approved drugs, and those in clinical trials were screened against the catalytic pocket of MTHFR in the presence and absence of cofactors. Binding energy and ADMET analysis revealed that Vilanterol (β2-adrenergic agonist), Selexipag (prostacyclin receptor agonist), and Ramipril Diketopiperazine (ACE inhibitor) are potential competitive inhibitors of MTHFR. Molecular dynamics analyses and MM-PBSA calculations with these compounds particularly revealed the amino acids between 285-290 for ligand binding and highlighted Vilanterol as the strongest candidate for MTHFR inhibition. Our results could guide the development of novel MTHFR inhibitor compounds, which could be inspired by the drugs brought into the spotlight here. More importantly, these potential candidates could be quhickly tested as a repurposing strategy in pre-clinical and clinical studies of the cancers mentioned above. Communicated by Ramaswamy H. Sarma


Introduction
One-carbon metabolism is an omnipresent process, which plays a pivotal role in many diseases including cancer and neurodegenerative as well as cardiovascular diseases.During this cellular process, single-carbon methyl substituents are transported to facilitate the synthesis of many key metabolites (Newman & Maddocks, 2017).DNA synthesis (thymidine and purines), amino acid homeostasis (serine, glycine, and methionine), redox balance, and epigenetic maintenance are all example outcomes of this critical metabolic process (Mattaini et al., 2016;Newman & Maddocks, 2017).Even though the relevance of one-carbon metabolism has long been acknowledged, new findings have further emphasized the prominent role of this pathway's homeostasis in cancer (Jain et al., 2012;Mattaini et al., 2016).Especially, the folate and methionine cycles, which are subsets of one-carbon metabolism, are known to adapt to cancer in a way that allows rapid synthesis of methyl groups that are consequently used in biosynthesis to maintain highly proliferative and persistent characteristics (Mazat, 2021;Tibbetts & Appling, 2010).
One of the major regulator proteins of this metabolism is methylenetetrahydrofolate reductase (MTHFR) (Chen et al., 2005).MTHFR is an enzyme that connects the folate and methionine cycles in one-carbon metabolism by catalyzing the conversion of 5,10-methylenetetrahydrofolate (CH 2 -THF) to 5methyltetrahydrofolate (CH 3 -THF) (Zheng et al., 2019).The methyl tetrahydrofolate is produced and then used to convert homocysteine to methionine (Froese et al., 2018;Figure 1a).As a result, MTHFR gene is indeed a critical hub gene in one-carbon metabolism (Zheng et al., 2019).MTHFR protein (a.k.a.5,10-methylenetetrahydrofolate) consists of 656 amino acids, having a molecular weight of 74597 Da (Goyette et al., 1998).It has a catalytic and regulatory domain and is predominantly found in a homodimer form (Bezerra et al., 2021).Some of its known ligands are flavin adenine dinucleotide (FAD) as a cofactor and nicotinamide adenine dinucleotide phosphate (NADP) as reducing agent, both of which bind to the catalytic domain (Burda et al., 2015).
Due to its emerging role in the metabolism of cell proliferation, folate and methionine metabolisms, MTHFR has been studied within the concept of potential cancer treatment strategies, with a focus on inhibiting DNA synthesis during uncontrolled cell proliferation (Stankova et al., 2008).Intriguing research showed that MTHFR knockdown caused remarkable inhibitory effects on the growth of human colon, prostate, lung, neuroblastoma, and breast tumor cells (Stankova et al., 2005) as well as gastric cancer cells (Sun et al., 2008).MHFTR was also shown to be the Achilles' heel of methionine-dependent cancer cells and it modulates the sensitivity of MYC-targeting therapies (Su et al., 2020).Various polymorphisms on the MTHFR gene have also been found to affect drug sensitivity and resistance (Kim, 2009;Yang et al., 2016).Moreover, in vivo studies on lung carcinoma demonstrated that reducing MTHFR expression shrank the tumor size nearly by half, where the apoptotic effects of cisplatin were enhanced (Stankova et al., 2005).In line with these experimental observations, MTHFR gene expression was shown to be associated with overall survival in liver hepatocellular carcinoma (LIHC) patients (Liu et al., 2019).When we performed further analysis on GEPIA web server (Tang et al., 2017), which utilizes The Cancer Genome Atlas (TCGA) datasets, we found significant links between overall survival and MTHFR gene expression not only in LIHC but also in adrenocortical carcinoma (ACC) and low-grade glioma (LGG) (Figure 1a).Additionally, there was a similar yet insignificant trend for acute myeloid leukemia (LAML).
Although the efficacy and therapeutic potential of MTHFR in cancer has been well demonstrated, MTHFR inhibiting compounds are yet to be identified.S-Adenosyl Methionine (SAM) was identified as the natural allosteric MTHFR inhibitor along with three sinefungin analogues (Bezerra et al., 2021;Bhatia et al., 2020).None of these reached clinics yet and no competitive inhibitors of MTHFR have yet been identified.The aim of this study is to find potential competitive inhibitor candidates against the MTHFR enzyme from currently existing small-drug molecules by applying the drug repurposing methodology.In order to achieve this, we first identified the druggable pockets in the MTHFR protein using the DoGSiteScorer tool (Volkamer et al., 2012), and figured out the catalytical domain as the most druggable site on MTHFR (Supplementary Figure 1).Next, we followed a pipeline, where 30470 molecules were obtained from the ZINC database (Sterling & Irwin, 2015) and scanned against the catalytical domain of MTHFR in the presence and absence of its cofactors (Figure 2).PyRx molecular docking-screening system, which is among the most cited in silico tools in the literature (Dallakyan & Olson, 2015), was employed.Small molecules determined by this process were subjected to ADMET (absorption, distribution, metabolism, excretion and toxicity) filters, where drug-efficacy and usability potential were further analyzed, and the results were refined (Guan et al., 2019).Finally, molecular dynamics simulations were performed for the small molecules that provided satisfactory results according to binding energies obtained from docking, and the results of the ADMET analysis; revealing potential competitive MTHFR inhibitors for clinical repurposing or further drug development.

Ligand library acquisition
In this study, we used ZINC15 database (https://zinc.docking.org/) (Sterling & Irwin, 2015), with a subset where only 'named' and 'commercially available' drugs were considered for screening against MTHFR.This subset contains FDA approved drugs, and those that are in clinical trials, alongside biogenic compounds.All (30470 molecules) were downloaded in the mol2 file format and then directly utilized for PyRx docking screen.Additionally, the substrate CH 2 -THF (ZINC ID: 4228244) was also downloaded to be used as a control.

Preparation of 3D structures
The crystal structure of MTHFR was obtained from Protein Data Bank (PDB ID: 6FCX) (Froese et al., 2018).Two separate models were used in the docking process, one that contained the cofactors FAD and NAD; and another one that was completely cleaned and did not contain any of these cofactors.Additionally, both models were pre-processed by removing all the water molecules using UCSF Chimera (Pettersen et al., 2004).To validate the pre-processed models, their structural alignment with the original PDB model was checked.For further validation, PROCHECK (Laskowski et al. 1993), and ProSA (Wiederstein & Sippl, 2007) were employed.As the first step of the screening procedure, the energies of downloaded drug molecules were minimized based on default parameters of PyRx (v0.9.9), which is the software used for virtual screening and docking (Dallakyan & Olson, 2015).They were then converted to the pdbqt format using OpenBabel, a program that converts chemical files to a format suitable for PyRx (Dallakyan & Olson, 2015).

PyRx molecular docking screen
For the docking procedure, the site of interest was selected as the binding site of the substrate, CH 2 -THF, as given in a previously published study (Froese et al., 2018).The center of mass for this site was calculated using all the atoms of the residues implicated in the binding site.The grid box utilized here comprised all the binding residues as given by Froese et al. (2018).Coordinates of the center of mass were calculated as À 7.7454, 44.0202 and À 31.1223 on X, Y and Z-axes, respectively.The dimensions were chosen as 23 Å x 23 Å x 23 Å, so that the box would contain all the atoms necessary for binding.Molecular docking was performed with PyRx (Dallakyan & Olson, 2015), which utilizes AutoDock Vina.In this step, docking was performed with the crystal structure containing the cofactors (FAD and NAD) in the binding domain and consecutively with another structure, which did not contain any of the cofactors.Docking was then performed for each downloaded molecule.The threshold for binding affinity was set to À 11.0 kcal/mol.

ADMET analysis
ADMET analysis (absorption, distribution, metabolism, excretion and toxicity) was performed using AdmetSar (Yang et al., 2017, Cheng et al., 2012) and SwissADME (Daina et al., 2017).Each molecule's data were extracted from ZINC15 database using ZINC IDs in the smiles file format.Firstly, smiles lists were presented to SwissAdme profiling tool for investigating the violations of Lipinski (Pfizer) (Lipinski et al., 2001) and Muegge (Bayer) (Muegge et al., 2001) filters.Some of the candidates that have filter violations higher or equal to two (Benet et al., 2016) and bioavailability scores lower than 0.55 were eliminated (Martin, 2005).For further investigation, AdmetSar was employed.Considering human oral bioavailability and human intestinal absorption potentials, candidates having lower acute oral toxicity levels (<2.5) were selected as previously described (Zhu et al., 2009).Finally, drugs that are not available in the market, based on ZINC database's 'available for sale' option, were eliminated.

Molecular dynamics analyses
Drugs and biogenic compounds that were found to be promising candidates for MTHFR inhibition were put in a molecular dynamics (MD) simulation with the MTHFR protein.In all MD simulations, AMBER99SB-ILDN all-atom force field (Lindorff-Larsen et al., 2010) was used and GROMACS 2021.1 software (Abraham et al., 2015) was employed.The docking result and its ligand complexes was used to initiate the MD simulation.To create ligand topologies and force field parameters, Acpype (AnteChamber PYthon Parser interfacE) (Silva & Vranken, 2012) was utilized.The system's topology was modified to include all compound parameters, which were merged into complex topology files.SPC/E water molecules were used to dissolve protein-ligand complexes in a cubic box.Then, to make the simulation system electrically neutral, we replaced solvent molecules with Cl -or Na þ ions.The particle-mesh Ewald (PME) method was used to treat long-range electrostatic interactions.During the simulations, the temperature was kept constant at 300 K. Using the Parrinello-Rahman coupling, the pressure was kept isotopically at 1 bar.The LINCS algorithm was used to constrain all bond lengths with a time step of 2 fs (femtoseconds).Structures were relaxed prior to MD simulations by performing 50000 steps of energy minimization using the steepest descent algorithm, which was followed by 1 ns of equilibration in the NVT (constant number of particles, system volume, and temperature) and NPT (constant number of particles, system pressure, and temperature) ensembles.Finally, unbiased MD simulations were run, where the system atomic coordinates were saved every 10 ps (picoseconds).Resulting trajectories were collected for further analysis.After the trajectories were obtained, GROMACS analysis toolkit was employed to analyze hydrogen bonds, root mean square deviation (RMSD), radius of gyration (Rg), and root mean square fluctuation (RMSF).The binding free energy (DG) was approximated using the Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) followed by decomposition of pre-residue by utilizing gmx_MMPBSA (Vald� es-Tresanco et al., 2021), which was adapted from AmberTools MMPBSA.py(Miller et al., 2012) script to work with GROMACS trajectories.The MMPBSA analysis consisted of 1000 frames and covered the final 50 ns of the simulation as previously described (Balbuena-Rebolledo et al., 2021).Per-residue decomposition is configured to include all residues within the 23 Å proximity of the ligand.

Preparation of the MTHFR structure and its validation
The processed crystal structure (PDB: 6FCX) was put through web-based structure validation tools PROCHECK and ProSA.PROCHECK checks the stereochemistry of a protein complex in detail.It produces a variety of charts in PostScript format and a detailed residue-by-residue list.
These provide an evaluation of the structure's overall quality in comparison to well-defined structures of the same resolution, along with highlighting residues that might also require additional examination (Laskowski et al., 1993).The Ramachandran plot obtained by PROCHECK showed that 87.8% of the residues were found in the most favored regions, whereas 11.2% were in the additional allowed regions with only 0.3% in disallowed regions (Figure 3a).Additionally, the model was structurally aligned with the original 6CFX structure and both models were found to be well-aligned, with an RMSD value of 0.004 Å (Figure 3b).Next, another validation was performed with ProSA, which is a convenient software for checking 3D models of proteins for possible faults.Error detection in experimental structures, hypothetical frameworks, and protein engineering are among its functions (Wiederstein & Sippl, 2007).The z-score by ProSA, which calculates the z-score of a given model and compares it to publicly available structures from PDB, was À 12.12 (Figure 3c).This score is within the range for the MTHFR model when compared to other proteins with similar sizes.Grid box coordinates were decided according to past research (Froese et al., 2018) and it was placed according to the binding site of the substrate, CH 2 -THF (Figure 3d).

Structure based ligand screening and molecular docking
Structure-based drug repurposing generates affinity estimation based on docking and binding simulations of target protein and potential candidates (Vyas et al., 2008).This method depends on cavity structural properties such as polarity, buriedness, hydrophobicity, length, and curvature for processing the 3D target structure (Rognan, 2017).For this purpose, MTHFR was docked using AutoDock Vina via the PyRx platform (Dallakyan & Olson, 2015), with more than 30470 molecules downloaded from the ZINC15 database.The calculated grid box shown in Figure 3d was used for the docking process.Docking was performed on two separate models with and without cofactors in order to only consider the molecules that exhibit high affinity both in the presence and absence of the cofactors.While evaluating the docking results, binding affinity, and RMSD values were taken into consideration.The more negative binding affinities are linked to better interaction between a ligand and a macromolecule (Dallakyan & Olson, 2015).Therefore, at the end of the screening, the binding affinity threshold of the selected molecules was determined to be À 11.00 (kcal/mol) and below, and only those with RMSD values '0' were considered.Additionally, the substrate CH 2 -THF was also docked to find the binding affinity of a molecule known to bind to the catalytic site to compare with the drug molecules.As a result, among the common molecules from the results of the screening with both models, those matching the specified criteria (as also demonstrated in Figure 2) were selected and listed in Table 1.Another supplementary table, showing additional molecules where the threshold is À 10.00 kcal/mol was also provided (Supplementary Table 1).

ADMET profiling
In order to examine the drug-efficacy of the small molecules selected via molecular docking; absorption, distribution, metabolism, elimination and toxicity (ADMET) properties were analyzed for each.At this stage, widely used tools, SwissADME (Daina et al., 2017) and AdmetSar (Yang et al., 2017, Cheng et al., 2012), both of which have been stated to be the gold standard of in silico ADMET profiling, were utilized (Kar & Leszczynski, 2020;Pathania & Singh, 2021).In the first step, the Lipinski (Lipinski et al., 2001) and Muegge filters (Muegge et al., 2001), crucial indicators of pharmacokinetic principles such as molecular weight, lipophilicity, and rotatable bonds that molecules should satisfy in order to be counted as a potential drug candidate (Benet et al., 2016)  were used as selective parameters.ADMET profiles derived from SwissAdme tool revealed that out of the eight molecules obtained from docking, Dabigatran Ethyl Ester (Zinc ID: 2043398), Ramipril Diketopiperazine (Zinc ID: 67665351), Selexipag (Zinc ID: 3990451), Vilanterol (Zinc ID: 3991624) and Rivenprost (Zinc ID: 3940680) showed overall satisfactory results (<1 violations) using the Lipinski filters.These five molecules also displayed convincing Abbott bioavailability scores (¼ 0.55) (Martin, 2005).ADMET results for these molecules were given in Table 2. ADMET profiles of molecules with the binding threshold lowered down to À 10.00 kcal/mol were given in Supplementary Table 2.
As the next step, to investigate the toxicity and further analyze the drug-efficacy features, remaining candidates that showed satisfactory results in the previous phase were presented to the AdmetSar tool.Further pharmacokinetic profile and toxicity assessment via AdmetSar showed that out of the 5 molecules that passed the Lipinski and Muegge filters, only Ramipril Diketopiperazine, Selexipag and Vilanterol showed satisfactory human oral bioavailability and human intestinal absorption characteristics.In order to increase the efficacy of potential inhibitors, these factors have been used as eliminative parameters, as human intestinal absorption and bioavailability are important factors in the efficacy of orally taken drugs (Guan et al., 2019).Further considering the acute oral toxicity scores (<2.5), three candidate molecules were determined among the eight molecules that satisfied the binding threshold, À 11.00 kcal/mol (Table 2).Lastly, to easily enable the application of these findings to possible future clinical research, it was also taken into account whether inhibitor candidates retrieved from ADMET analysis were already available on the market or not.As such, Dabigatran Ethyl Ester was rejected from the list for being unpurchasable.
In conclusion, Vilanterol, Selexipag and Ramipril Diketopiperazine were selected as potential inhibitor candidates (Figure 4).The LigPlotþ (Wallace et al., 1995) representations show the MTHFR amino acids, which interact with drug molecules and the hydrogen bonds between the protein and each of the drugs.All the molecules were found to interact with Tyr285, Leu235, Glu27 and Gln192 residues.Additionally, Ala159, Thr58, Thr191, Ile190, His91, Trp59, His60, Met118, Thr191 and Ile190 were found to be commonly interacting with the substrate and at least one of the drug molecules.

Molecular dynamics simulations
Molecular dynamics (MD) was performed to analyze the dynamics and stability of most promising repurposing candidates (Vilanterol, Selexipag and Ramipril Diketopiperazine) for MTHFR inhibition along with the natural substrate CH 2 -THF (Figure 5).Candidate inhibitor molecules and the The structures were taken from PubChem (Kim et al., 2021).
Manzamine A 2 violations 2 violations 0.17 À 1.026 þ substrate were subjected to an MD simulation of 100 ns together with the MTHFR homodimer.
The Root-Mean-Square deviation (RMSD) of the backbone measures the stability of the structure and atomic change conformation upon ligand binding (Aier et al., 2016).RMSD results of our MD simulations revealed similar characteristics with the natural substrate and inhibitor candidates (Figure 5a).The radius of gyration (Rg) is a measure of the compactness and size of protein structures.Rg value is calculated by the distribution of atoms of a protein around its axis (Arnittali et al., 2019), facilitating the estimation of the pressure exerted on a specific location and showing the strength of a bond between two cross-sections.The radius of gyration decreases as the mass is dispersed closer to the axis of rotation (Sneha & Priya Doss, 2016).Our Rg analysis revealed that the substrate CH 2 -THF reached the lowest Rg values and all the inhibitor candidates had similarly stable results (Figure 5b).The RMSD values of ligands show how much a ligand changes conformation inside the binding pocket after initial binding.While CH 2 -THF fluctuated during the whole simulation time and peaking at around  (Wallace et al., 1995) representations.The yellow lines represent H bonds, whereas red arcs represent Van der Waals interactions between the amino acid residues and drug molecules.
2 nm, Selexipag's RMSD consistently increased, and Vilanterol and Ramipril Diketopiperazine exhibited a similarly stable profile (Figure 5c).The Root Mean Square Fluctuation (RMSF) values, on the other hand, stayed below 0.7 nm for the whole simulation, and the highest peak was observed between residues 120-140 for all ligands (Figure 5d).Additionally, the capabilities of CH 2 -THF and the drug molecules to form hydrogen bonds were also analyzed to assess the conformational changes between the ligand and protein based on time (Figure 6).The H-bonds provide necessary stability to the protein-ligand complex (Kumar et al., 2014).The number of Hbonds between MTHFR and the ligands was calculated by GROMACS at 300 K.For the substrate, the number of hydrogen bonds fluctuated with an average of 0.869.Selexipag and Vilanterol had higher averages at 1.313 and 1.077 respectively, whereas Ramipril Diketopiperazine averaged at 0.140.

Binding free energy analyses with MM-PBSA and per-residue decomposition
Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) calculations were used to examine the binding affinities between protein-ligand complexes that were obtained from MD trajectories (Table 3).The results of MM-PBSA displayed the same order of affinities as molecular docking.The substrate CH 2 -THF exhibited lower affinity towards the MTHFR in comparison to Vilanterol, Selexipag and Ramipril Diket.as was the case with molecular docking.
We also performed decomposition analyses, which revealed the energy contribution of individual residues during the molecular dynamics simulation (Figure 7).According to this calculation, Vilanterol and Selexipag revealed a higher number of interacting residues (jDGj >0.5 kcal/mol) within 23 Å of their respective ligand whereas CH 2 -THF and Ramipril Diket.were only associated with a handful of residues.In line with this, decomposition result for Vilanterol and Selexipag generally agreed well with the corresponding residues in the molecular docking result (Figure 4), even though this relation did not seem as strong for Ramipril Diket.and CH 2 -THF.Moreover, certain residues particularly stood out among the decomposition analyses for all ligands.Specifically; LEU235 appeared for all drugs but not for the substrate, whereas LEU287 was emphasized for the substrate itself as well as Ramipril Diket.and Selexipag.Interestingly, another residue (i.e.TYR285), which is very close to LEU287, was identified to be interacting with Vilanterol; altogether highlighting the potential importance of the amino acids between the positions 285-290.

Discussion
Finding new therapeutic potentials for available drugs to use them for disease treatment other than what they were first described for, is referred to as drug repositioning or repurposing.It is a process that is significant for the development of new drugs, as it can be more cost-effective and time- saving when compared to classical drug discovery (Jarada et al., 2020).This process has been used to find possible drugs for the treatment of cancer (Schein, 2021;Zhang et al., 2020), rare diseases (Roessler et al., 2021) and most recently, Covid-19 (Chakraborty et al., 2021).Drug repurposing trials for Covid-19 illuminated the need for finding rapid treatments during the height of the pandemic.Some of the repurposed antiviral drugs have been approved in several countries for the treatment of Covid-19 (Chakraborty et al., 2021), and several drugs are in clinical trials for cancer treatment (Nowak-Sliwinska et al., 2019).
The contribution of MTHFR in human disease and the need for its inhibitors is well acknowledged.Even though deleterious effects were reported for the complete knockout of MTHFR (Chen et al., 2001), knockdown studies suggest possibly beneficial effects particularly against human cancers in vitro and in vivo (Stankova et al., 2005(Stankova et al., , 2008;;Wang et al., 2022;Wu et al., 2021).While testing MTHFR inhibition preclinically and clinically will address the safety concerns, it is worth mentioning that the complete knockouts of many genes/proteins, whose inhibitors are successfully and widely used in the clinics at the moment, were reported to cause lethal effects whereas the right dosing regimen of inhibition brings clinical benefits.Sorafenib (Escalante & Zalpour, 2011) and Erlotinib (Masuda et al., 2017) are well known examples for this case as the knockouts of their target proteins are reportedly lethal (Haiko et al., 2008;Lee & Threadgill, 2009;Zhang et al., 2010).Essential roles of MTHFR in subcellular metabolic events including carbon metabolism indeed make it an interesting drug target.More recently, MTHFR was linked to efflux transportation with potential implications on drug resistance (Naseem et al., 2022) pointing out a  Newly developed and novel allosteric inhibitors were already recently described for MTHFR (Bezerra et al., 2021;Bhatia et al., 2020), yet no experimentally or clinically used competitive inhibitor was developed.The effectiveness of either type of inhibitor can change depending on the context.Interestingly, competitive and allosteric inhibitors can be used concurrently to overcome drug resistance (R� ea & Hughes, 2022).For instance, both competitive and allosteric inhibitors have been described for AKT with varying side effects and separate effectiveness rates for different disease contexts as extensively reviewed previously (Lazaro et al., 2020).In this study, we performed an in silico drug screen against the catalytical domain using docking and molecular dynamics simulations, and found three candidates for repurposing: Vilanterol, a b2-adrenergic agonist (Wen et al., 2022); Selexipag, a prostacyclin receptor (IP) agonist (Kaufmann et al., 2015), and Ramipril Diketopiperazine, an ACE inhibitor (Kn€ utter et al., 2008).During the evaluation of the binding affinity results, we saw that the substrate had binding affinities of À 7.6 and À 8.6 kcal/mol for models with cofactors and without cofactors, respectively; while the binding affinities were À 12 and À 14.1 kcal/mol for Vilanterol, À 11.8 and À 13.8 kcal/mol for Selexipag, and À 11.4 and À 11.5 kcal/mol for Ramipril Diketopiperazine.(Table 1).As the docking of the substrate was performed as a control, and its binding affinity was found lower than the drug molecules in all analyses described in this article, the molecules pointed out in our study have high potentials to bind to MTHFR even though this findings need to be elaborated further with experimental studies.
Hydrogen bonds are thought to be the main contributors to higher binding affinities (D.Chen et al., 2016), interestingly Ramipril Diket., which shows no H bonds in the LigPlot þ results and an average of 0.14 bonds in the MD results, has higher binding affinity when compared to the substrate.This higher binding affinity might be related to increased hydrophobic and Van der Waals interactions, the latter of which can be seen in LigPlot þ results where Ramipril Diketopiperazine has 14 and the substrate has 10 (Figure 4).The LigPlot þ results also show different numbers of H bonds for the substrate and Vilanterol in comparison to MD results.In LigPlot þ representations, the substrate and Vilanterol had 4 H bonds, Selexipag had 1 and Ramipril Diketopiperazine had 0 (Figure 4).On the other hand in the MD results, Ramipril Diketopiperazine had an average of 0.140 and Selexipag had 1.313 (Figure 6).Conversely, MD simulations showed that the substrate had an average of 0.869 and Vilanterol had 1.077 (Figure 6).The underlying reason might be the fact that in LigPlotþ, the 3 D structures were flattened to a 2 D representation, leaving the side chains with limited movement capacity (Wallace et al., 1995).As this is not a moving simulation unlike MD, the snapshot of the ligand-protein complex might have had the maximum bonds at the time.In fact, for Vilanterol and Selexipag particularly, the formation of Hbonds showed high variability at 50 ns level (Figure 6b and c).In line with this, RMSD ligand results (Figure 5c) also showed that the molecules had higher instability at that time.The RMSF results showed peaks around the 120-140 residue areas for all the molecules with the exception of Vilanterol (Figure 5d).This area contains FAD binding residues, which, as the MD tests have been done with the model that did not contain cofactors, may have led to the decreased stability of the molecules that are bound to the area.
When detailed calculations were performed with MM-PBSA using the MD trajectories we obtained, we observed a significant reduction in binding affinities for the substrate and all drug molecules except Vilanterol, as opposed to docking results.As we saw this reduction for the substrate as well, this result might be attributed to known limitations of molecular dynamics simulations as discussed in the literature (Hollingsworth & Dror, 2018).Still, decomposition analyses revealed some consistent residues, particularly highlighting the amino acids between 285-290 for ligand binding.In future studies, laboratory experiments are necessary to investigate the bindings and subsequent effects of these drug molecules in vitro and validate the in silico results presented in this study possibly via a receiver operating characteristic (ROC) analysis (Al-Sha'er et al., 2022).

Conclusion
In this study, our purpose was to reveal potential competitive inhibitors for MTHFR, which is a crucial protein related to proliferation in cancer cells.We specifically aimed to reveal molecules that are clinically available, bypassing the long period needed to develop de novo drugs.To this end, we performed systematic molecular docking followed by ADMET analysis and molecular dynamics.Based on our analyses, we determined three drug molecules which can possibly be used as competitive inhibitors of MTHFR: Selexipag, Vilanterol and Ramipril Diketopiperazine, particularly highlighting the amino acids between 285-290 for binding.Among the reported candidates, Vilanterol stood out as it showed consistently high levels of binding both in molecular docking and molecular dynamics simulations.

Figure 2 .
Figure 2. Pipeline flowchart of the study.

Figure 3 .
Figure 3. Validation of prepared MTHFR structure model.(a) Ramachandran plot of the processed model.(b) Structural alignment of the model with the original PDB structure.(c) z-score plot of the model within the range of the original structures.(d) The placement of the grid box for molecular docking.

Figure 4 .
Figure 4. Natural substrate and potential competitive MTHFR inhibitors for repurposing.The binding simulations of the (a) CH 2 -THF substrate, (b) Vilanterol, (c) Selexipag and (d) Ramipril Diketopiperazine along with their schematic LigPlotþ(Wallace et al., 1995) representations.The yellow lines represent H bonds, whereas red arcs represent Van der Waals interactions between the amino acid residues and drug molecules.

Figure 5 .
Figure 5. Molecular dynamics analyses of the natural substrate and potential MTHFR inhibitor candidates.(a) RMSD values of the Ca atoms of the backbone of MTHFR upon different ligand binding as a function of time.(b) Radius of gyration of the three drugs and substrate.(c) RMSD values of the ligands.(d) RMSF of all the ligands at different residues.

Table 1 .
Molecules obtained as a result of docking screen and their binding affinities.

Table 2 .
ADMET analysis of the molecules obtained from molecular docking.

Table 3 .
Calculations of the binding free energy of MTHFR complexes using the mmPBSA method.
tantalizing possibility of employing MTHFR inhibition in reducing drug resistance in cancer.