Taylor & Francis Taylor & Francis Group
In silico drug screen reveals potential competitive MTHFR inhibitors for clinical repurposing
Nazlıgül Keske, Başak Özay, Ezgi Yağmur Tükel, Muratcan Menteş & Cihangir Yandım
To cite this article: Nazlıgül Keske, Başak Özay, Ezgi Yağmur Tükel, Muratcan Menteş & Cihangir Yandım (2023) In silico drug screen reveals potential competitive MTHFR inhibitors for clinical repurposing, Journal of Biomolecular Structure and Dynamics, 41:21, 11818-11831, DOI: 10.1080/07391102.2022.2163697
To link to this article: https://doi.org/10.1080/07391102.2022.2163697
+
View supplementary material ☒
Published online: 04 Jan 2023.
Submit your article to this journal ☒
Article views: 558
View related articles
☒
CrossMark
View Crossmark data
☒
Citing articles: 1 View citing articles ☒
Taylor & Francis Taylor & Francis Group
Check for updates
In silico drug screen reveals potential competitive MTHFR inhibitors for clinical repurposing
Nazlıgül Keskeª*, Başak Özaya*, Ezgi Yağmur Tükela*, Muratcan Menteşª* and Cihangir Yandıma,b
ªFaculty of Engineering, Department of Genetics and Bioengineering, İzmir University of Economics, Balçova, İzmir, Turkey; bİzmir Biomedicine and Genome Center (IBG), Dokuz Eylül University Health Campus, İnciraltı, İzmir, Turkey
Communicated by Ramaswamy H. Sarma
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 car- cinoma, 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 iden- tify 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 stron- gest candidate for MTHFR inhibition. Our results could guide the development of novel MTHFR inhibi- tor compounds, which could be inspired by the drugs brought into the spotlight here. More importantly, these potential candidates could be quickly tested as a repurposing strategy in pre-clin- ical and clinical studies of the cancers mentioned above.
Abbreviations: ACC: Adrenocortical carcinoma; ADMET: Absorption, distribution, metabolism, excretion and toxicity; CH2-THF: 5,10-methylenetetrahydrofolate; CH3-THF: 5-methyltetrahydrofolate; FAD: Flavin adenine dinucleotide; FDA: U.S. Food and drug administration; Fs: Femtoseconds; HCC: Hepatocellular Carcinoma; IP: Prostacyclin receptor; LAML: Acute Myeloid Leukemia; LGG: Low-grade glioma; LIHC: Liver hepatocellular carcinoma; MD: Molecular Dynamics; MM-PBSA: Molecular Mechanics Poisson- Boltzmann Surface Area; MTHFR: Methylenetetrahydrofolate reductase; NAD: Nicotinamide adenine dinucleotide; NPT: Constant number of particles, system pressure, and temperature; NVT: Constant number of particles, system volume, and temperature; PDB: Protein data bank; PME: Particle-mesh Ewald; PRMT4: Arginine methyltransferase 4; Ps: Picoseconds; RG: Radius of gyration; RMSD: Root mean square deviation; RMSF: Root mean square fluctuation; ROC: Receiver operating characteristic; SAM: S- Adenosyl methionine; TCGA: The Cancer Genome Atlas
ARTICLE HISTORY Received 28 July 2022 Accepted 24 December 2022
KEYWORDS Docking; MTHFR; cancer; molecular dynamics; drug screen; PyRx; repurposing; inhibitor
1. 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 metabo- lites (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 conse- quently 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.,
CONTACT Cihangir Yandım Economics, Balçova, İzmir, 35330, Turkey
Faculty of Engineering, Department of Genetics and Bioengineering, İzmir University of
*These authors contributed equally.
+ Supplemental data for this article can be accessed online at https://doi.org/10.1080/07391102.2022.2163697.
2005). MTHFR is an enzyme that connects the folate and methionine cycles in one-carbon metabolism by catalyzing the conversion of 5,10-methylenetetrahydrofolate (CH2-THF) to 5- methyltetrahydrofolate (CH3-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-car- bon 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 prolif- eration, folate and methionine metabolisms, MTHFR has
A
Folic Acid
THF
Methionine
S -adenosylmethionine
5,10 Methylene THF
Thymi- dine
Purines
ONE-CARBON CYCLE
Methylation
FOLATE CYCLE
S -adenosylhomocysteine
5-Methyl THF
Homocysteine
MTHFR
· DNA
· Histones
· DNA Repair
· Epigenetic imprinting
B
ACC
LGG
1.0
Low MTHFR TPM High MTHFR TPM
1.0
Logrank p=0.0061 HR(high)=3.1
Low MTHFR TPM High MTHFR TPM
Logrank p=0.0045
0.8
p(HR)=0.0096
0.8
HR(high)=1.7
Percent survival
n(high)=38 n(low)=38
Percent survival
p(HR)=0.0049 n(high)=257 n(low)=257
0.6
0.6
0.4
0.4
0.2
0.2
0.0
0.0
0
50
100
150
0
50
100
150
200
Months
Months
LIHC
LAML
1.0
Low MTHFR TPM
1.0
High MTHFR TPM
Low MTHFR TPM
High MTHFR TPM
Logrank p=0.0026
Logrank p=0.32
0.8
HR(high)=1.7
HR(high)=1.3
p(HR)=0.0028
0.8
p(HR)=0.33
Percent survival
n(high)=181 n(low)=181
Percent survival
n(high)=53
0.6
0.6
n(low)=53
0.4
0.4
0.2
0.2
0.0
0.0
0
20
40
60
80
100
120
0
20
40
60
80
Months
Months
been studied within the concept of potential cancer treat- ment 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 carcin- oma 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 hep- atocellular 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 inhibi- tor candidates against the MTHFR enzyme from currently existing small-drug molecules by applying the drug repur- posing methodology. In order to achieve this, we first identi- fied 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 data- base (Sterling & Irwin, 2015) and scanned against the catalyt- ical domain of MTHFR in the presence and absence of its cofactors (Figure 2). PyRx molecular docking-screening sys- tem, which is among the most cited in silico tools in the lit- erature (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 fur- ther drug development.
30k+ Molecules From ZINC15 Database
₸
Crystal Structure of MTHFR From PDB (PDB ID: 6FCX)
Pre-Processing of The Crystal Structure
Structure Based In Silico Screening (PyRx 0.9.9, AutoDock Vina)
Binding Energy (≤-11.00 kcal/mol)
Rejected Binding Energy (> -11.00 kcal/mol)
Rejected ADMET Analysis (Violations >1, NO Oral Bioavail.)
ADMET Analysis (Violations <1, Oral toxicity < 2.5)
Molecular Dynamics
2. Materials and method
2.1. 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, along- side biogenic compounds. All (30470 molecules) were down- loaded in the mol2 file format and then directly utilized for PyRx docking screen. Additionally, the substrate CH2-THF (ZINC ID: 4228244) was also downloaded to be used as a control.
2.2. 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 con- tained 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 for- mat suitable for PyRx (Dallakyan & Olson, 2015).
2.3. PyRx molecular docking screen
For the docking procedure, the site of interest was selected as the binding site of the substrate, CH2-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 calcu- lated 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 per- formed for each downloaded molecule. The threshold for binding affinity was set to -11.0 kcal/mol.
2.4. ADMET analysis
ADMET analysis (absorption, distribution, metabolism, excre- tion 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 vio- lations 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 elimi- nated (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.
2.5. Molecular dynamics analyses
Drugs and biogenic compounds that were found to be promis- ing 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 mol- ecules with CI 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 (femtosec- onds). Structures were relaxed prior to MD simulations by per- forming 50000 steps of energy minimization using the steepest descent algorithm, which was followed by 1ns of equilibration in the NVT (constant number of particles, system volume, and temperature) and NPT (constant number of par- ticles, 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 trajec- tories were obtained, GROMACS analysis toolkit was employed to analyze hydrogen bonds, root mean square deviation (RMSD), radius of gyration (Rg), and root mean square fluctu- ation (RMSF). The binding free energy (AG) was approximated using the Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) followed by decomposition of pre-residue by utilizing gmx_MMPBSA (Valdés-Tresanco et al., 2021), which was adapted from AmberTools MMPBSA.py (Miller et al., 2012) script to work with GROMACS trajectories. The MMPBSA ana- lysis 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.
3. Results
3.1. 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 reso- lution, 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 ori- ginal 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
A
180
B
b
PDB:6FCX
Model
135
62
b
-1
:
90
GLU 128 (B)
Psi (degrees)
45
a
HIS 2799 (B)
0
-45ª
-90
TYR 13
(13)
ASP
(A)
-135
ASP 133 (B)
b.
b
-180
-135
-90
-45
0
45
90
135
180
Phi (degrees)
C
D
10
X-ray
NMR
5
0
Z-score
-5
-10
-15
-20
0
200
400
600
800
1000
Number of residues
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, CH2-THF (Figure 3d).
3.2. Structure based ligand screening and molecular docking
Structure-based drug repurposing generates affinity estima- tion based on docking and binding simulations of target pro- tein and potential candidates (Vyas et al., 2008). This method depends on cavity structural properties such as polarity, buri- edness, hydrophobicity, length, and curvature for processing the 3D target structure (Rognan, 2017). For this purpose, MTHFR was docked using AutoDock Vina via the PyRx plat- form (Dallakyan & Olson, 2015), with more than 30470 mole- cules downloaded from the ZINC15 database. The calculated grid box shown in Figure 3d was used for the docking pro- cess. Docking was performed on two separate models with and without cofactors in order to only consider the mole- cules 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 mol- ecules was determined to be -11.00 (kcal/mol) and below, and only those with RMSD values ‘0’ were considered. Additionally, the substrate CH2-THF was also docked to find the binding affinity of a molecule known to bind to the cata- lytic 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 add- itional molecules where the threshold is -10.00 kcal/mol was also provided (Supplementary Table 1).
3.3. 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 uti- lized (Kar & Leszczynski, 2020; Pathania & Singh, 2021). In the first step, the Lipinski (Lipinski et al., 2001) and Muegge fil- ters (Muegge et al., 2001), crucial indicators of pharmacoki- netic 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)
| Zinc ID | Compound name | Structure | Binding energy (kcal/mol) | |
|---|---|---|---|---|
| With cofactors | Without cofactors | |||
| zinc000004726511 zinc000003991624 | Benzyl(trityl)sulfane Vilanterol | -12.9 -12 | -13.2 -14.1 | |
| zinc000003916953 | Mazokalim | -11.9 | -11.3 | |
| zinc000003990451 | Selexipag | -11.8 | -13.8 | |
| 0=$=0 | ||||
| zinc000003940680 zinc000004769727 | Rivenprost | -11.5 | -12 | |
| Manzamine A | -11.5 | -11 | ||
| zinc000002043398 | Dabigatran Ethyl Ester | -11.4 | -11.7 | |
(continued)
| Zinc ID | Compound name | Structure | Binding energy (kcal/mol) | |
|---|---|---|---|---|
| With cofactors | Without cofactors | |||
| zinc000067665351 | Ramipril Diketopiperazine | -11.4 | -11.5 | |
| 0 O H 0 H | ||||
| zinc00004228244 | CH2-THF (Substrate) | N O H 9 HO NHZ HO 0 0 NHI 0 | -7.6 | -8.6 |
The structures were taken from PubChem (Kim et al., 2021).
| ZINC ID | Compound name | ADMET ANALYSIS | |||||
|---|---|---|---|---|---|---|---|
| Lipinski filter (Pfizer) | Muegge filter (Bayer) | Bioavailability score | Human oral bioavailability | Acute oral toxicity | Human intestinal absorption | ||
| zinc000002043398 | Dabigatran Ethyl Ester | 0 violation | 0 violation | 0.55 | - | 2.34 | þ |
| zinc000067665351 | Ramipril Diketopiperazine | 0 violation | 0 violation | 0.55 | þ | 1.952 | þ |
| zinc000003990451 | Selexipag | 0 violation | 0 violation | 0.55 | þ | 2.117 | þ |
| zinc000004726511 | Benzyl(trityl)sulfane | 1 violation | 2 violations | 0.55 | þ | 2.701 | - |
| zinc000003916953 | Mazokalim | 1 violation | 1 violation | 0.55 | - | 3.119 | þ |
| zinc000003991624 | Vilanterol | 0 violation | 1 violation | 0.55 | þ | 2.074 | þ |
| zinc000003940680 | Rivenprost | 0 violation | 0 violation | 0.55 | - | 3.96 | - |
| zinc000004769727 | Manzamine A | 2 violations | 2 violations | 0.17 | - | 1.026 | þ |
were used as selective parameters. ADMET profiles derived from SwissAdme tool revealed that out of the eight mole- cules 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 mole- cules 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 pre- sented 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 mole- cules 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 candi- dates (Figure 4). The LigPlot+ (Wallace et al., 1995) represen- tations show the MTHFR amino acids, which interact with drug molecules and the hydrogen bonds between the pro- tein 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 com- monly interacting with the substrate and at least one of the drug molecules.
3.4. Molecular dynamics simulations
Molecular dynamics (MD) was performed to analyze the dynamics and stability of most promising repurposing candi- dates (Vilanterol, Selexipag and Ramipril Diketopiperazine) for MTHFR inhibition along with the natural substrate CH2- THF (Figure 5). Candidate inhibitor molecules and the
A
Ala159
Tyr285
mmmm
Tle190
His283
Thr191
Met118
CH2-THF
Glu27
TYR 285
LLLLLLLL His91
n
Gin 192
Leu235
Thr58
THR 58
HIS 283
Pro30
PRO 30
His60
Trp59
B
Thr58
Leu287
Ala159
w
Vilanterol
Thr93
Thr191
Glu27
GLU 27
Leu120
Leu235
Gly122
GLN 192
LEU 235
Gin192
Tyr161
TYR 285
Tyr285
Ile190
C
Gin231
Phe29
Leu235
u
Glu27
Leu287
Selexipag
4
His60
m
Tyr285
2
His91
Gin192
3
HIS 91
STALLILLU
Trp59
Thr58
D
Glu27
Ile190
Ramipril Diket.
Tyr285
Thr191
Gin192
A
Leu235
Met118
Leu120
Gly122
Tyr161
Trp59
Asp123
Thr93
His60
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 con- formation 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 distribu- tion 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 CH2-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 CH2-THF fluctu- ated during the whole simulation time and peaking at around
A
Root-Mean-Square Deviation of Backbone
B
Radius of Gyration
4
0.6
3.95
0.5
3.9
RMSD (nm)
3.85
0.4
Rg (nm)
3.8
0.3
3.75
0.2
CH2-THF
3.7
- CH2-THF
Vilanterol
- Vilanterol
0.1
- Selexipag
- Ramipril Diket.
3.65
- Selexipag
- Ramipril Diket.
0
0
25
50
75
3.6
Time (ns)
100
0
20000
40000
Time (ps)
60000
80000
le+05
C
Root-Mean-Square Deviation of Ligands
3
D
Root-Mean-Square Fluctuation
0.8
CH2-THF
CH2-THF
2.5
Vilanterol
0.7
- Vilanterol
- Selexipag
- Selexipag
- Ramipril Diket.
0.6
- Ramipril Diket.
RMSD (nm)
2
0.5
1.5
(nm)
0.4
0.3
1
0.2
0.5
0.1
0
0
25
50
75
0
Time (ns)
100
50
100
150
200
Residue
250
300
2 nm, Selexipag’s RMSD consistently increased, and Vilanterol and Ramipril Diketopiperazine exhibited a similarly stable pro- file (Figure 5c). The Root Mean Square Fluctuation (RMSF) val- ues, on the other hand, stayed below 0.7 nm for the whole simulation, and the highest peak was observed between resi- dues 120-140 for all ligands (Figure 5d). Additionally, the capa- bilities of CH2-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 pro- tein-ligand complex (Kumar et al., 2014). The number of H- bonds 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.
3.5. 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 affin- ities between protein-ligand complexes that were obtained from MD trajectories (Table 3). The results of MM-PBSA dis- played the same order of affinities as molecular docking. The substrate CH2-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 dur- ing the molecular dynamics simulation (Figure 7). According to this calculation, Vilanterol and Selexipag revealed a higher number of interacting residues (|4G| >0.5 kcal/mol) within 23 Å of their respective ligand whereas CH2-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 CH2-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 posi- tions 285-290.
4. 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 repur- posing. It is a process that is significant for the development of new drugs, as it can be more cost-effective and time-
A
CH2-THF
B
Vilanterol
6
6
5
5
4
4
Number
Number
w
w
2
2
1
1
0
0
20000
40000
60000
80000 1e+05
0
Time (ps)
0
20000 40000 60000 80000 1e+05
Time (ps)
C
Selexipag
Ramipril Diket.
6
D
6
5
5
4
4
Number
Number
w
2
2
1
1
0
0
20000
40000
60000
80000
le+05
0
Time (ps)
0
20000
40000
60000
Time (ps)
80000
le+05
| CH2-THF | Vilanterol | Selexipag | Ramipril Diket. | |||||
|---|---|---|---|---|---|---|---|---|
| Avg | SEM | Avg | SEM | Avg | SEM | Avg | SEM | |
| VDWAALS | -28.35 | 0.13 | 40.89 | 0.11 | 47.10 | 0.14 | -38.14 | 0.09 |
| EEL | 223.07 | 0.79 | 305.35 | 0.51 | 35.46 | 0.66 | -2.10 | 0.08 |
| EPB | 231.85 | 0.83 | 315.28 | 0.46 | 50.45 | 0.55 | 17.23 | 0.10 |
| ENPOLAR | -19.00 | 0.08 | 31.94 | 0.06 | -36.65 | 0.09 | -26.60 | 0.05 |
| EDISPER | 35.35 | 0.14 | 51.71 | 0.08 | 62.01 | 0.13 | 44.73 | 0.06 |
| AG gas | -251.42 | 0.82 | 346.24 | 0.52 | -82.56 | 0.66 | -40.23 | 0.12 |
| AG solv | 248.19 | 0.86 | 335.05 | 0.47 | 75.82 | 0.54 | 35.36 | 0.12 |
| AG Total | -3.23 | 0.16 | -11.19 | 0.17 | -6.74 | 0.25 | -4.88 | 0.12 |
*All values given are in kcal/mol.
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 treat- ments 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 treat- ment (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, 2008; Wang et al., 2022; Wu et al., 2021). While testing MTHFR inhibition pre- clinically 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
CH2-THF
Vilanterol
-6
-6
4
-4
Energy (kcal/mol)
Energy (kcal/mol)
-2
-2
1
T
1
1
0
0
1
L
1
2
2
8
8
6
A:PHE:29
A:HIS:60
A:HIS:227
A:LEU:287
6
A:GLU:27
A:ASP:56
A:TYR:161
A:PRO:166
A:GLN:192
A:LEU:193
A:PHE:221
A:GLN:231
A:LEU:232
A:LEU:235
A:SER:236
A:LYS:237
A:TYR:285
Selexipag
Ramipril Diket.
-
6
-6
4
4
Energy (kcal/mol)
Energy (kcal/mol)
-
2
1
2
L
0
Ml
1
1
1
L
0
T
1
2
2
4
8
6
A:GLU:27
A:PHE:29
A:ARG:32
A:ASP:56
A:THR:58
A:HIS:60
A:ASP:123
A:GLN:231
A:LEU:232
A:LYS:234
A:LEU:235
A:LEU:287
6
A:GLU:27
A:PHE:29
A:LEU:235
A:THR:286
A:LEU:287
A:ASN:288
tantalizing possibility of employing MTHFR inhibition in reducing drug resistance in cancer.
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 con- text. Interestingly, competitive and allosteric inhibitors can be used concurrently to overcome drug resistance (Réa & 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 repur- posing: Vilanterol, a ß2-adrenergic agonist (Wen et al., 2022); Selexipag, a prostacyclin receptor (IP) agonist (Kaufmann et al., 2015), and Ramipril Diketopiperazine, an ACE inhibitor (Knütter et al., 2008). During the evaluation of the binding affinity results, we saw that the substrate had binding affin- ities of -7.6 and -8.6 kcal/mol for models with cofactors and without cofactors, respectively; while the binding affin- ities 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 anal- yses 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 dif- ferent numbers of H bonds for the substrate and Vilanterol in comparison to MD results. In LigPlot + representations, the sub- strate and Vilanterol had 4H 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 simu- lations 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 move- ment capacity (Wallace et al., 1995). As this is not a moving simulation unlike MD, the snapshot of the ligand-protein com- plex might have had the maximum bonds at the time. In fact,
for Vilanterol and Selexipag particularly, the formation of H- bonds 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 inves- tigate 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).
5. Conclusion
In this study, our purpose was to reveal potential competi- tive 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 anal- yses, we determined three drug molecules which can pos- sibly be used as competitive inhibitors of MTHFR: Selexipag, Vilanterol and Ramipril Diketopiperazine, particularly high- lighting 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.
Acknowledgement
The high-end computer analyses in this study were performed at TUBITAK ULAKBIM, High Performance and Grid Computing Center (TRUBA Resources). We sincerely thank TUBITAK for providing this platform.
Disclosure statement
No potential conflict of interest was reported by the authors.
Funding
This study was funded by the 2209 A program of The Scientific and Technological Research Council of Turkey (TUBITAK).
Author contributions
Conceptualization (CY); Data curation (CY); Formal analysis (NK, BÖ, EYT, MM, CY); Funding acquisition (NK, BÖ, EYT, CY); Investigation (NK, BÖ, EYT, MM, CY); Methodology (NK, BÖ, EYT, MM, CY); Project administra- tion (CY); Resources (CY); Software (CY); Supervision (CY); Validation (NK, BÖ, EYT, MM, CY); Visualization (NK, BÖ, EYT, MM, CY); Roles/Writing - original draft (NK, BÖ, EYT, MM, CY); Writing - review & editing (CY).
References
Abraham, M. J., Murtola, T., Schulz, R., Páll, S., Smith, J. C., Hess, B., & Lindah, E. (2015). Gromacs: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX, 1-2, 19-25. https://doi.org/10.1016/j.softx.2015.06.001
Aier, I., Varadwaj, P. K., & Raj, U. (2016). Structural insights into conform- ational stability of both wild-type and mutant EZH2 receptor. Scientific Reports, 6, 34984. https://doi.org/10.1038/srep34984
Al-Sha’er, M. A., Basheer, H. A., & Taha, M. O. (2022). Discovery of new PKN2 inhibitory chemotypes via QSAR-guided selection of docking- based pharmacophores. Molecular Diversity, https://doi.org/10.1007/ s11030-022-10434-4
Arnittali, M., Rissanou, A. N., & Harmandaris, V. (2019). Structure of bio- molecules through molecular dynamics simulations. Procedia Computer Science, 156, 69-78. https://doi.org/10.1016/j.procs.2019.08. 181
Balbuena-Rebolledo, I., Padilla-Martínez, I. I., Rosales-Hernández, M. C., & Bello, M. (2021). Repurposing FDA drug compounds against breast cancer by targeting egfr/her2. Pharmaceuticals, 14(8), 791. https://doi. org/10.3390/ph14080791
Benet, L. Z., Hosey, C. M., Ursu, O., & Oprea, T. I. (2016). BDDCS, the Rule of 5 and drugability. In Advanced drug delivery reviews (Vol. 101, pp.89-98). Elsevier B.V. https://doi.org/10.1016/j.addr.2016.05.007
Bezerra, G. A., Holenstein, A., Foster, W. R., Xie, B., Hicks, K. G., Bürer, C., Lutz, S., Mukherjee, A., Sarkar, D., Bhattacharya, D., Rutter, J., Talukdar, A., Brown, P. J., Luo, M., Shi, L., Froese, D. S., & Yue, W. W. (2021). Identification of small molecule allosteric modulators of 5,10-methyle- netetrahydrofolate reductase (MTHFR) by targeting its unique regula- tory domain. Biochimie, 183, 100-107. https://doi.org/10.1016/j.biochi. 2021.01.007
Bhatia, M., Thakur, J., Suyal, S., Oniel, R., Chakraborty, R., Pradhan, S., Sharma, M., Sengupta, S., Laxman, S., Masakapalli, S. K., & Bachhawat, A. K. (2020). Allosteric inhibition of MTHFR prevents futile SAM cycling and maintains nucleotide pools in one-carbon metabolism. The Journal of Biological Chemistry, 295(47), 16037-16057. https://doi.org/ 10.1074/jbc.RA120.015129
Burda, P., Schäfer, A., Suormala, T., Rummel, T., Bürer, C., Heuberger, D., Frapolli, M., Giunta, C., Sokolová, J., Vlášková, H., Kožich, V., Koch, H. G., Fowler, B., Froese, D. S., & Baumgartner, M. R. (2015). Insights into severe 5,10-methylenetetrahydrofolate reductase deficiency: Molecular genetic and enzymatic characterization of 76 patients. Human Mutation, 36(6), 611-621. https://doi.org/10.1002/humu.22779
Chakraborty, C., Sharma, A. R., Bhattacharya, M., Agoramoorthy, G., & Lee, S. S. (2021). The drug repurposing for COVID-19 clinical trials pro- vide very effective therapeutic combinations: Lessons learned from major clinical studies. Frontiers in Pharmacology, 12, 704205. https:// doi.org/10.3389/fphar.2021.704205
Chen, D., Oezguen, N., Urvil, P., Ferguson, C., Dann, S. M., & Savidge, T. C. (2016). Regulation of protein-ligand binding affinity by hydrogen bond pairing. Science Advances, 2(3), 433-443. https://doi.org/10.1126/ sciadv.1501240
Chen, Z., Karaplis, A. C., Ackerman, S. L., Pogribny, I. P., Melnyk, S., Lussier-Cacan, S., Chen, M. F., Pai, A., John, S. W. M., Smith, R. S., Bottiglieri, T., Bagley, P., Selhub, J., Rudnicki, M. A., James, S. J., & Rozen, R. (2001). Mice deficient in methylenetetrahydrofolate reduc- tase exhibit hyperhomocysteinemia and decreased methylation cap- acity, with neuropathology and aortic lipid deposition. Human Molecular Genetics, 10(5), 433-443. https://doi.org/10.1093/hmg/10.5. 433
Chen, J., Gammon, M. D., Chan, W., Palomeque, C., Wetmur, J. G., Kabat, G. C., Teitelbaum, S. L., Britton, J. A., Terry, M. B., Neugut, A. I., & Santella, R. M. (2005). One-carbon metabolism, MTHFR polymor- phisms, and risk of breast cancer. Cancer Research, 65(4), 1606-1614. https://doi.org/10.1158/0008-5472.CAN-04-2630
Cheng, F., Li, W., Zhou, Y., Shen, J., Wu, Z., Liu, G., Lee, P. W., & Tang, Y. (2012). admetSAR: A comprehensive source and free tool for assess- ment of chemical ADMET properties. Journal of Chemical Information and Modeling, 52(11), 3099-3105. https://doi.org/10.1021/ci300367a
Daina, A., Michielin, O., & Zoete, V. (2017). SwissADME: A free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules. Scientific Reports, 7, 42717. https://doi. org/10.1038/srep42717
Dallakyan, S., & Olson, A. J. (2015). Small-molecule library screening by docking with PyRx. Surveys in High Energy Physics, 14(1), 145-173. https://doi.org/10.1080/01422419908228843
Escalante, C. P., & Zalpour, A. (2011). Vascular endothelial growth factor inhibitor-induced hypertension: Basics for primary care providers. Cardiology Research and Practice, 2011, 816897. https://doi.org/10. 4061/2011/816897
Froese, D. S., Kopec, J., Rembeza, E., Bezerra, G. A., Oberholzer, A. E., Suormala, T., Lutz, S., Chalk, R., Borkowska, O., Baumgartner, M. R., & Yue, W. W. (2018). Structural basis for the regulation of human 5,10- methylenetetrahydrofolate reductase by phosphorylation and S- adenosylmethionine inhibition. Nature Communications, 9(1), 2261. https://doi.org/10.1038/s41467-018-04735-2
Goyette, P., Pai, A., Milos, R., Frosst, P., Tran, P., Chen, Z., Chan, M., & Rozen, R. (1998). Gene structure of human and mouse methylenete- trahydrofolate reductase (MTHFR). Mammalian Genome: Official Journal of the International Mammalian Genome Society, 9(8), 652-656. https://doi.org/10.1007/s003359900838
Guan, L., Yang, H., Cai, Y., Sun, L., Di, P., Li, W., Liu, G., & Tang, Y. (2019). ADMET-score-A comprehensive scoring function for evaluation of chemical drug-likeness. MedChemComm, 10(1), 148-157. https://doi. org/10.1039/C8MD00472B
Haiko, P., Makinen, T., Keskitalo, S., Taipale, J., Karkkainen, M. J., Baldwin, M. E., Stacker, S. A., Achen, M. G., & Alitalo, K. (2008). Deletion of vas- cular endothelial growth factor C (VEGF-C) and VEGF-D is not equiva- lent to VEGF receptor 3 deletion in mouse embryos. Molecular and Cellular Biology, 28(15), 4843-4850. https://doi.org/10.1128/mcb. 02214-07
Hollingsworth, S. A., & Dror, R. O. (2018). Molecular dynamics simulation for all. Neuron, 99(6), 1129-1143. https://doi.org/10.1016/j.neuron. 2018.08.011
Jain, M., Nilsson, R., Sharma, S., Madhusudhan, N., Kitami, T., Souza, A. L., Kafri, R., Kirschner, M. W., Clish, C. B., & Mootha, V. K. (2012). Metabolite profiling identifies a key role for glycine in rapid cancer cell proliferation. Science (New York, N.Y.), 336(6084), 1040-1044. https://doi.org/10.1126/science.1218595
Jarada, T. N., Rokne, J. G., & Alhajj, R. (2020). A review of computational drug repositioning: Strategies, approaches, opportunities, challenges, and directions. Journal of Cheminformatics, 12 (1), 1-23. https://doi. org/10.1186/s13321-020-00450-7
Kar, S., & Leszczynski, J. (2020). Open access in silico tools to predict the ADMET profiling of drug candidates. Expert Opinion on Drug Discovery, 15 (12), 1473-1487. https://doi.org/10.1080/17460441.2020. 1798926
Kaufmann, P., Okubo, K., Bruderer, S., Mant, T., Yamada, T., Dingemanse, J., & Mukai, H. (2015). Pharmacokinetics and tolerability of the novel oral prostacyclin IP receptor agonist selexipag. American Journal of Cardiovascular Drugs: Drugs, Devices, and Other Interventions, 15(3), 195-203. https://doi.org/10.1007/s40256-015-0117-4
Kim, S., Chen, J., Cheng, T., Gindulyte, A., He, J., He, S., Li, Q., Shoemaker, B. A., Thiessen, P. A., Yu, B., Zaslavsky, L., Zhang, J., & Bolton, E. E. (2021). PubChem in 2021: New data content and improved web inter- faces. Nucleic Acids Research, 49(D1), D1388-D1395. https://doi.org/10. 1093/nar/gkaa971
Kim, Y. I. (2009). Role of the MTHFR polymorphisms in cancer risk modifi- cation and treatment. Future Oncology (London, England), 5(4), 523- 542. https://doi.org/10.2217/fon.09.26
Knütter, I., Wollesky, C., Kottra, G., Hahn, M. G., Fischer, W., Zebisch, K., Neubert, R. H. H., Daniel, H., & Brandsch, M. (2008). Transport of angiotensin-converting enzyme inhibitors by H +/peptide transport- ers revisited. The Journal of Pharmacology and Experimental Therapeutics, 327(2), 432-441. https://doi.org/10.1124/jpet.108.143339
Kumar, K. M., Anbarasu, A., & Ramaiah, S. (2014). Molecular docking and molecular dynamics studies on ß-lactamases and penicillin binding proteins. Molecular bioSystems, 10(4), 891-900. https://doi.org/10. 1039/c3mb70537d
Laskowski, R. A., MacArthur, M. W., Moss, D. S., & Thornton, J. M. (1993). PROCHECK: a program to check the stereochemical quality of protein structures. Journal of applied crystallography, 26(2), 283-291.
Lazaro, G., Kostaras, E., & Vivanco, I. (2020). Inhibitors in AKTion: ATP- competitive vs allosteric. Biochemical Society Transactions, 48(3), 933- 943. https://doi.org/10.1042/BST20190777
Lee, T. C., & Threadgill, D. W. (2009). Generation and validation of mice carrying a conditional allele of the epidermal growth factor receptor. Genesis (New York, N.Y. : 2000), 47(2), 85-92. https://doi.org/10.1002/ dvg.20464
Lindorff-Larsen, K., Piana, S., Palmo, K., Maragakis, P., Klepeis, J. L., Dror, R. O., & Shaw, D. E. (2010). Improved side-chain torsion potentials for the Amber ff99SB protein force field. Proteins, 78(8), 1950-1958. https://doi.org/10.1002/prot.22711
Lipinski, C. A., Lombardo, F., Dominy, B. W., & Feeney, P. J. (2001). Experimental and computational approaches to estimate solubility and permeability in drug discovery and development q settings. Advanced Drug Delivery Reviews, 46(1-3), 3-26. www.elsevier.com/ locate/drugdeliv https://doi.org/10.1016/S0169-409X(00)00129-0
Liu, G. M., Zeng, H. D., Zhang, C. Y., & Xu, J. W. (2019). Identification of a six-gene signature predicting overall survival for hepatocellular carcin- oma. Cancer Cell International, 19, 138. https://doi.org/10.1186/s12935- 019-0858-2
Martin, Y. C. (2005). A bioavailability score. Journal of Medicinal Chemistry, 48(9), 3164-3170. https://doi.org/10.1021/jm0492002
Masuda, C., Yanagisawa, M., Yorozu, K., Kurasawa, M., Furugaki, K., Ishikura, N., Iwai, T., Sugimoto, M., & Yamamoto, K. (2017). Bevacizumab counteracts VEGF-dependent resistance to erlotinib in an EGFR-mutated NSCLC xenograft model. International Journal of Oncology, 51(2), 425-434. https://doi.org/10.3892/ijo.2017.4036
Mattaini, K. R., Sullivan, M. R., & Vander Heiden, M. G. (2016). The import- ance of serine metabolism in cancer. Journal of Cell Biology, 214 (3), 249-257. https://doi.org/10.1083/jcb.201604085
Mazat, J .- P. (2021). One-carbon metabolism in cancer cells: A critical review based on a core model of central metabolism. Biochemical Society Transactions, 49(1), 1-15. https://doi.org/10.1042/BST20190008
Miller, B. R. I. I. I., McGee, T. D., Jr., Swails, J. M., Homeyer, N., Gohlke, H., & Roitberg, A. E. (2012). MMPBSA.py: An efficient program for end- state free energy calculations. Journal of Chemical Theory and Computation, 8(9), 3314-3321. https://doi.org/10.1021/ct300418h
Muegge, I., Heald, S. L., & Brittelli, D. (2001). Simple selection criteria for drug-like chemical matter. Journal of Medicinal Chemistry, 44(12), 1841-1846. https://doi.org/10.1021/jm015507e
Naseem, A., Pal, A., Gowan, S., Asad, Y., Donovan, A., Temesszentandrási- Ambrus, C., Kis, E., Gaborik, Z., Bhalay, G., & Raynaud, F. (2022). Intracellular metabolomics identifies efflux transporter inhibitors in a routine Caco-2 cell permeability assay-Biological implications. Cells, 11(20), 3286. https://doi.org/10.3390/cells11203286
Newman, A. C., & Maddocks, O. D. K. (2017). One-carbon metabolism in cancer. British Journal of Cancer, 116(12), 1499-1504. https://doi.org/ 10.1038/bjc.2017.118
Nowak-Sliwinska, P., Scapozza, L., & i Altaba, A. R. (2019). Drug repurpos- ing in oncology: Compounds, pathways, phenotypes and computa- tional approaches for colorectal cancer. Biochimica et Biophysica Acta - Reviews on Cancer, 1871 (2), 434-454. https://doi.org/10.1016/j. bbcan.2019.04.005
Pathania, S., & Singh, P. K. (2021). Analyzing FDA-approved drugs for compliance of pharmacokinetic principles: should there be a critical screening parameter in drug designing protocols ?. Expert Opinion on Drug Metabolism and Toxicology, 17 (4), 351-354. https://doi.org/10. 1080/17425255.2021.1865309
Pettersen, E. F., Goddard, T. D., Huang, C. C., Couch, G. S., Greenblatt, D. M., Meng, E. C., & Ferrin, T. E. (2004). UCSF Chimera - A visualiza- tion system for exploratory research and analysis. Journal of Computational Chemistry, 25(13), 1605-1612. https://doi.org/10.1002/ jcc.20084
Réa, D., & Hughes, T. P. (2022). Development of asciminib, a novel allo- steric inhibitor of BCR-ABL1. Critical Reviews in Oncology/Hematology, 171, 103580. https://doi.org/10.1016/j.critrevonc.2022.103580
Roessler, H. I., Knoers, N. V. A. M., van Haelst, M. M., & van Haaften, G. (2021). Drug repurposing for rare diseases. Trends in Pharmacological Sciences, 42 (4), 255-267. https://doi.org/10.1016/j.tips.2021.01.003
Rognan, D. (2017). The impact of in silico screening in the discovery of novel and safer drug candidates. Pharmacology & Therapeutics, 175, 47-66. https://doi.org/10.1016/j.pharmthera.2017.02.034
Schein, C. H. (2021). Repurposing approved drugs for cancer therapy. British Medical Bulletin, 137 (1), 13-27. https://doi.org/10.1093/bmb/ ldaa045
Silva, D., & Vranken, B. F. (2012). ACPYPE-AnteChamber PYthon Parser interfacE. BMC Research Notes, 5(1), 1-8. http://www.biomedcentral. com/1756-0500/5/367
Sneha, P., & Priya Doss, C. G. (2016). Molecular dynamics: New frontier in personalized medicine. Advances in Protein Chemistry and Structural Biology, 102, 181-224. https://doi.org/10.1016/bs.apcsb.2015.09.004
Stankova, J., Lawrance, A., & Rozen, R. (2008). Methylenetetrahydrofolate reductase (MTHFR): A novel target for cancer therapy. Current Pharmaceutical Design, 14(11), 1143-1150. https://doi.org/10.2174/ 138161208784246171
Stankova, J., Shang, J., & Rozen, R. (2005). Antisense inhibition of methyl- enetetrahydrofolate reductase reduces cancer cell survival in vitro and tumor growth in vivo. Clinical Cancer Research: An Official Journal of the American Association for Cancer Research, 11(5), 2047-2052. https://doi.org/10.1158/1078-0432.CCR-04-2047
Sterling, T., & Irwin, J. J. (2015). ZINC 15 - Ligand discovery for everyone. Journal of Chemical Information and Modeling, 55(11), 2324-2337. https://doi.org/10.1021/acs.jcim.5b00559
Su, A., Ling, F., Vaganay, C., Sodaro, G., Benaksas, C., Dal Bello, R., Forget, A., Pardieu, B., Lin, K. H., Rutter, J. C., Bassil, C. F., Fortin, G., Pasanisi, J., Antony-Debré, I., Alexe, G., Benoist, J .- F., Pruvost, A., Pikman, Y., Qi, J., … Puissant, A. (2020). The folate cycle enzyme MTHFR is a critical regulator of cell response to MYC-targeting therapies. Cancer Discovery, 10(12), 1894-1911. https://doi.org/10.1158/2159-8290.CD- 19-0970
Sun, D. F., Weng, Y. R., Chen, Y. X., Lu, R., Wang, X., & Fang, J. Y. (2008). Knock-down of methylenetetrahydrofolate reductase reduces gastric cancer cell survival: An in vitro study. Cell Biology International, 32(8), 879-887. https://doi.org/10.1016/j.cellbi.2008.03.018
Tang, Z., Li, C., Kang, B., Gao, G., Li, C., & Zhang, Z. (2017). GEPIA: A web server for cancer and normal gene expression profiling and inter- active analyses. Nucleic Acids Research, 45(W1), W98-W102. https://doi. org/10.1093/nar/gkx247
Tibbetts, A. S., & Appling, D. R. (2010). Compartmentalization of mamma- lian folate-mediated one-carbon metabolism. Annual Review of Nutrition, 30, 57-81. https://doi.org/10.1146/annurev.nutr.012809. 104810
Valdés-Tresanco, M. S., Valdés-Tresanco, M. E., Valiente, P. A., & Moreno, E. (2021). gmx_MMPBSA: A New Tool to Perform End-State Free Energy Calculations with GROMACS. Journal of Chemical Theory and Computation, 17(10), 6281-6291. https://doi.org/10.1021/acs.jctc. 1c00645
Volkamer, A., Kuhn, D., Rippmann, F., & Rarey, M. (2012). Dogsitescorer: A web server for automatic binding site prediction, analysis and drugg- ability assessment. Bioinformatics (Oxford, England), 28(15), 2074-2075. https://doi.org/10.1093/bioinformatics/bts310
Vyas, V., Jain, A., Jain, A., & Gupta, A. (2008). Virtual screening: A fast tool for drug design. Scientia Pharmaceutica, 76(3), 333-360. https://doi. org/10.3797/scipharm.0803-03
Wallace, A. C., Laskowski, R. A., & Thornton, J. M. (1995). LIGPLOT: A pro- gram to generate schematic diagrams of protein-ligand interactions The LIGPLOT program automatically generates schematic 2-D repre- sentations of protein-ligand complexes from standard Protein Data Bank file input. Protein Engineering, 8 (2), 127-134. http://peds.oxford- journals.org/
Wang, X., Xu, Z., Ren, X., Chen, X., Yi, Q., Zeng, S., Thakur, A., Gong, Z., & Yan, Y. (2022). MTHFR inhibits TRC8-mediated HMOX1 ubiquitination and regulates ferroptosis in ovarian cancer. Clinical and Translational Medicine, 12(9), e1013. https://doi.org/10.1002/ctm2.1013
Wen, Y. C., Chen, W. Y., Tram, V. T. N., Yeh, H. L., Chen, W. H., Jiang, K. C., Abou-Kheir, W., Huang, J., Hsiao, M., & Liu, Y. N. (2022). Pyruvate kin- ase L/R links metabolism dysfunction to neuroendocrine differenti- ation of prostate cancer by ZBTB10 deficiency. Cell Death and Disease, 13(3), 252. https://doi.org/10.1038/s41419-022-04694-z
Wiederstein, M., & Sippl, M. J. (2007). ProSA-web: Interactive web service for the recognition of errors in three-dimensional structures of pro- teins. Nucleic Acids Research, 35(Web Server), W407-W410. https://doi. org/10.1093/nar/gkm290
Wu, M. T., Ye, W. T., Wang, Y. C., Chen, P. M., Liu, J. Y., Tai, C. K., Tang, F. Y., Li, J. R., Liu, C. C., & Chiang, E. P. I .. (2021). MTHFR knockdown assists cell defense against folate depletion induced chromosome segregation and uracil misincorporation in DNA. International Journal of Molecular Sciences, 22(17), 9392. https://doi.org/10.3390/ ijms22179392
Yang, H., Lou, C., Sun, L., L., J., Cai, Y., Wang, Z., Li, W., Liu, G., & Tang, Y. (2017). admetSAR 2.0: Web-service for prediction and optimization of chemical ADMET properties. Bioinformatics, 33(1), 1-7. https://doi.org/ 10.1093/bioinformatics/xxxxxx
Yang, L., Wang, X. W., Zhu, L. P., Wang, H. L., Wang, B., Wu, T., Zhao, Q., Jinsihan, D. L. X. T., & Wang, X. Y. (2016). Relationship between gen- etic polymorphisms of methylenetetrahydrofolate reductase and breast cancer chemotherapy response. Genetics and Molecular Research, 15(3), https://doi.org/10.4238/gmr.15038679
Zhang, Z., Zhou, L., Xie, N., Nice, E. C., Zhang, T., Cui, Y., & Huang, C. (2020). Overcoming cancer therapeutic bottleneck by drug repurpos- ing. Signal Transduction and Targeted Therapy, 5 (1). 1-25. https://doi. org/10.1038/s41392-020-00213-8
Zheng, Y., Ramsamooj, S., Li, Q., Johnson, J. L., Yaron, T. M., Sharra, K., & Cantley, L. C. (2019). Regulation of folate and methionine metabolism by multisite phosphorylation of human methylenetetrahydrofolate reductase. Scientific Reports, 9(1), 4190. https://doi.org/10.1038/ s41598-019-40950-7
Zhang, L., Zhou, F., Han, W., Shen, B., Luo, J., Shibuya, M., & He, Y. (2010). VEGFR-3 ligand-binding and kinase activity are required for lymphangiogenesis but not for angiogenesis. Cell Research, 20(12), 1319-1331. https://doi.org/10.1038/cr.2010.116
Zhu, H., Martin, T. M., Ye, L., Sedykh, A., Young, D. M., & Tropsha, A. (2009). Quantitative structure - Activity relationship modeling of rat acute toxicity by oral exposure. Chemical Research in Toxicology, 22(12), 1913-1921. https://doi.org/10.1021/tx900189p