p-ISSN 1693-5683; e-ISSN 2527-7146 84 Vol. 20, No. 1, May 2023, pp. 84-89 Research Article Molecular Dynamics Study of Human MTHFR Wild-Type and A222V Mutant-Type: A Computational Approach Michael Resta Surya Yanuar, Dita Maria Virginia, Enade Perdana Istyastono* Faculty of Pharmacy, Sanata Dharma University, Campus III Paingan, Maguwoharjo, Depok, Sleman, Yogyakarta, 55282, Indonesia https://doi.org/10.24071/jpsc.005727 J. Pharm. Sci. Community, 2023, 20(1), 84-89 Article Info ABSTRACT Received: 11-01-2023 Revised: 13-01-2023 Accepted: 14-01-2023 *Corresponding author: Enade Perdana Istyastono email: enade@usd.ac.id Keywords: Molecular dynamics simulations; MTHFR; YASARA-Structure Methylenetetrahydrofolate reductase (MTHFR) is an essential enzyme for the homeostasis and metabolism of intracellular folate. However, it is known that one of the commonly found MTHFR gene SNPs (A222V) causes a decrease in activity which decreases folate levels and increases the accumulation of homocysteine in the blood. The purpose of this study was to investigate the molecular dynamics of MTHFR wild- type protein and MTHFR A222V mutant protein in silico. After some preparations, the crystal structure of MTHFR with code 6FCX obtained from the Protein Data Bank was used as the input file for molecular dynamics (MD) simulations. YASARA-Structure was employed for performing 50 ns production run of MD simulations. The deviation of the root-mean-squared deviation value of the backbone atoms (∆RMSDBb) of MTHFR wild-type is consistently less than 2 Å from the beginning of the MD simulations production runs. On the other hand, there is a sudden spike of ∆RMSDBb of MTHFR A222V from 15.01 to 20.00 ns of 2.221 Å. According to the analysis of ∆RMSDBb, MTHFR wild-type protein is considered stable and the MTHFR A222V mutant protein is considered unstable which may lead to various clinical problems for the person possessing this mutation. INTRODUCTION Methylenetetrahydrofolate reductase (MTHFR) is a FAD-dependent flavoenzyme so each subunit of the MTHFR protein binds non- covalently to Flavin Adenine Dinucleotide (FAD). FAD acts as a redox cofactor in several metabolic reactions and NADPH is required as an electron donor. This protein consists N-terminal catalytic region and C-terminal regulatory region which allosteric inhibitor S-adenosylmethionine (SAM) binds and regulates enzyme activity in response to cellular methionine concentration which is very slow and can be reversed by binding to S- adenosylhomocysteine (SAH) (Abhinand et al., 2016; Froese et al., 2018). MTHFR is an important enzyme in the metabolism and homeostasis of intracellular folate. MTHFR catalyzes the irreversible change of 5,10- methylenetetrahydrofolate to 5- methyltetrahydrofolate which is the primary methyl donor for homocysteine remethylation to methionine in 1-C cycle. The dimer MTHFR protein is encoded by the MTHFR gene located at chromosomal locus 1p36.22 of 20373 bp containing 12 exons (NCBI, 2022; Raghubeer et al.., 2021). However, one of the commonly found Single Nucleotide Polymorphism (SNP) of the MTHFR gene is at point 677 C→T which is located in the 4th exon where it causes changes on the 222th amino acid sequence from Alanine to Valine (Castro et al., 2003). The genetic variation of MTHFR A222V results in a thermolabile form of the MTHFR protein it is also known to experience a decrease in its ability to bind to its cofactor, FAD (Abhinand et al., 2016). Genetic variation MTHFR A222V homozygous is known to cause a decrease in activity by 50-60% and can act synergistically if there is a mutation at another point (Abhinand et al., 2016; El Hajj Chehadeh et al., 2016). This decrease in activity is a crucial factor in the http://issn.pdii.lipi.go.id/issn.cgi?daftar&1180428136&1&& http://issn.pdii.lipi.go.id/issn.cgi?daftar&1465346481&1&& https://e-journal.usd.ac.id/index.php/JFSK/index https://creativecommons.org/licenses/by/4.0/ Journal of Pharmaceutical Sciences and Community Molecular Dynamics Study of Human MTHFR... Research Article 85 Yanuar et al. J. Pharm. Sci. Community, 2023, 20(1), 84-89 decrease of folate levels and the accumulation of homocysteine in the blood (El Hajj Chehadeh et al., 2016). The relationship between the MTHFR A222V polymorphism and disease might be associated with two aspects involving homocysteine concentration and altered metabolism of folates and homocysteine (Ueland et al., 2001). The disease that might be involved such as diabetic neuropathy (Wu et al., 2017; Yigit et al., 2013), congenital heart disease (CHD)(Xuan et al., 2014), renal failure (Hishida et al., 2013), cancer (He et al., 2017; Kumar et al., 2018; Li et al., 2021), psychiatric disorders (Zhang et al., 2022), etc. This genetic variation requires further exploration regarding the risk of changes in protein structure and stability. Therefore, in silico studies are important to confirm the structural stability through molecular dynamics (MD) simulations of mutated MTHFR protein. METHODS Materials The crystal structure of MTHFR (6FCX.pdb) obtained from Protein Data Bank (https://www.rcsb.org/structure/6FCX accessed on 25 July 2022) was used as the starting material for MD simulation input files. Instrumentation The molecular dynamics were performed using YASARA-Structure 20.10.4 in Windows 10 Pro personal computer client with 4 Intel(R) Core (TM) i5-7500 CPU @ 3.41 GHz as the processor, RAM 8.00 GB. Input file preparation The input file 6FCX.pdb was loaded to YASARA-Structure 20.10.4 in the client computer with the SeqRes assigned as “Yes” in the loading options. Only amino acids from chains A, FAD as a cofactor, and SAH as a ligand were used in the MD simulations so unnecessary atoms in the MD simulations were removed. The SeqRes module showed missing residues, i.e., Met37, Asp38, Pro39 from the N-terminal loop and Ile161, Gly162, Asp163, Gln164, Trp165, Glu166, Glu167, Glu168, Glu169, Gly170, Gly171, Ser392, Ser393, Ser394, Pro395, Ala396 from the loop. Therefore, the command BuildLoop None, Sequence=MDPE,5-7, Structures=1, Mutate=All, Bumpsum=1.0, SecStr= followed by BuildLoop 949-951, Sequence=PIGDQWEEEEGGF,967-969, Structures=1, Mutate=All, Bumpsum=1.0, SecStr= and BuildLoop 2767-2769, Sequence=NSSSPAF,2780-2782, Structures=1, Mutate=All, Bumpsum=1.0, SecStr= were applied. Before the terminal caps were added, atoms in Gln651 is manually edited to arrange the C-terminal loop with the command DelAtom 9836-9839 and SwapAtom 9834, Oxygen, UpdateBonds=Yes, UpdateHyd=Yes at first. The command pH value=7.4, update=Yes, and CleanAll followed by AddCapObj 1, Type=ACE+NME, Location=all were subsequently performed to add hydrogens correctly at the physiological pH then add terminal caps. The Force Field was set to AMBER14 and the energy minimization experiment was performed with the simulation cell constructed in the cubic shape with a distance of 5 Å around all atoms. The scene was saved as 222A-429A.sce. We included MTHFR E429A in our study as the native protein collected from 6FCX.pdb input file and MTHFRThe scene was then modified by swapping amino acids Ala to Glu at 429th order with the command SwapRes Atom 6268, Glu, Isomer=L to make wild-type protein model of MTHFR as control. The scene was saved as 222A- 429E.sce. The scene was modified again by swapping amino acids Ala to Val at 222th order with the command SwapRes Atom 2890, Val, Isomer=L to make A222V mutated protein model of MTHFR. The scene was saved as 222V- 429E.sce. Molecular dynamics simulations The MD simulations in YASARA-Structure were performed by a macro file that defines the parameter of the simulations using personal computer client. A macro file md_run_50ns_ss10ps.mcr was prepared by modifying the default macro md_run.mcr. It was done to perform 50 ns MD simulations during which snapshots were taken every 10 ps at 310K and water density 0.993 g/mL. The following are the modified simulation parameters by Istyastono et al. (2022), described in YASARA-Structure used in this study: The simulation was performed on YASARA-Structure. The setup included optimization of the hydrogen bond network to increase the solute stability and a pKa prediction to fine-tune the protonation state of protein residues at the chosen pH of 7.4. NaCl ions were added at a physiological concentration of 0.9%, (an excess of either Na or Cl was added to neutralize the cell). After the steepest descent and simulated annealing minimization to eliminate collisions, the simulation was run for 50 ns using the AMBER14 force field for the solute, GAFF2, and AM1BCC for ligands, and TIP3P for water. For Van der Waals https://www.rcsb.org/structure/6FCX Research Article Journal of Pharmaceutical Sciences and Community Molecular Dynamics Study of Human MTHFR ... 86 Yanuar et al. J. Pharm. Sci. Community, 2023, 20(1), 84-89 forces, the cut-off was 8 Å and no cut-off was applied to electrostatic forces (using the Particle Mesh Ewald algorithm). The equations of motion were integrated with multiple timesteps of 1.25 fs for bonded interactions and 2.5 fs for non- bonded interactions at a temperature of 310K and a pressure of 1 atm (NPT ensemble) using algorithms described in detail previously. After examining the solute root-mean-squared deviation (RMSD) as a function of simulation time, the first 7 ns were considered equilibration time and excluded from further analysis. (Istyastono et al., 2022; Krieger et al., 2015). Analysis The first 7 ns were considered equilibration time, while the subsequent 43 ns were considered as the production runs. The snapshots resulting from the MD simulations were analyzed using the default macro “md_analyze.mcr”. Together with all the pdb files resulted from running the “md_convert.mcr” macro, the files “dpp4-kiq-r2md_analysis.tab” and “dpp4-apo-r2md_analysis.tab” resulted from running the “md_analyze.mcr” macro was copied to the computer client for further analysis. The stability of the molecular systems was analyzed by following suggestions by (Liu et al., 2017). RESULTS AND DISCUSSION The RMSD values of the backbone atoms of the MTHFR E429A, MTHFR wild type, and MTHFR A222V were obtained from the files “222A-429A_analysis.tab”, “222A- 429E_analysis.tab”, and “222A- 429E_analysis.tab” respectively. The graphs RMSDBb vs simulation from 0 ns to 50 ns is presented in Figure 1. It shows that these three molecular systems have reached equilibrium at 7 ns simulation time. Therefore, the subsequent 43 ns in the molecular dynamics are considered as the production runs. According to Liu et al. (2017), the protein-ligand complex is considered stable if the deviation of the RMSD value is less than 2 Å in MD simulations (Liu et al., 2017). In this experiment, the stability of the MTHFR backbone atoms was analyzed by examining the deviation of the RMSD of the backbone atoms (∆RMSDBb) values in every 5 ns of the total 50 ns simulation time as recommended by Istyastono et al. (Istyastono et al., 2021). Figure 2 shows the average value of ∆RMSDBb in every 5 ns during the MD simulations of three molecular systems. We can see that the backbone atoms of the MTHFR wild type were reported stable consistently from the beginning of the MD simulations production runs. Even if there is a spike in 15.01 – 20 .00 ns but it is still in the acceptable range of ∆RMSD value a protein- ligand complex is considered stable and subsequently stabilizes itself until it reached 50,00 ns. On the other hand, the ∆RMSDBb value resulting from mutated MTHFR E429A protein- ligand complex shows stabilizing its complex by decreased value of ∆RMSDBb from 0 to 20 ns but there are spikes occur from 25.01 – 30.00 ns and 35.01 – 40.00 ns of 2.453 Å and 2.092 Å respectively that considered unstable. ∆RMSDBb value resulted from mutated MTHFR A222V also indicating that this protein-ligand complex is unstable because there is a sudden spike of ∆RMSDBb from 15.01 to 20.00 ns of 2.221 Å. Figure 1. The RMSD values of the MTHFR backbone atoms of the MD simulations Journal of Pharmaceutical Sciences and Community Molecular Dynamics Study of Human MTHFR... Research Article 87 Yanuar et al. J. Pharm. Sci. Community, 2023, 20(1), 84-89 Figure 2. The average value of the RMSDBb deviation in every 5 ns (∆RMSDBb) during the MD simulations Instrument to Assess Misconception Levels This result supports the previous studies that MTHFR A222V might be involved in several diseases as mentioned in various studies described below. There are significant differences between genotype and allele frequencies of MTHFR A222V with diabetic peripheral neuropathy (p=0.003) and retinopathy (p=0.039) (Yigit et al., 2013). To support that, a meta-analysis by Wu et al (2016) reported that there is an association between MTHFR A222V and diabetic neuropathy in Turkey population (OR=1.43; 95% CI 1.08-1.90; p=0.014)(Wu et al., 2017). Xuan et al (2014) conducted a meta-analysis based on 9329 cases and 15076 controls to observe the association between MTHFR A222V gene polymorphism and the risk for CHD. They detected that there was a significant association in all genetic models of overall children populations for T vs C (OR 51.248, 95% CI: 1.093–1.426; P = 0.001), TT vs CC (OR 51.485, 95% CI: 1.140–1.935; P = 0.003), and TT vs CT (OR 51.312, 95% CI: 1.100–1.565; P = 0.003), the dominant model (OR 5 1.240, 95% CI: 1.053–1.461; P = 0.010), and the recessive model (OR 5 1.410, 95% CI: 1.139–1.724; P = 0.001. In addition, significant associations were observed in the overall maternal population in all genetic models for T vs C (OR 51.215, 95% CI: 1.085–1.361; P = 0.001), TT vs CC (OR 51.488, 95% CI: 1.169–1.859; P = 0.001), TT vs CT (OR 5 1.315, 95% CI: 1.042–1.659; P = 0.021), the dominant model (OR 5 1.258, 95% CI: 1.144– 1.383; P = 2.14e-6), and the recessive model (OR 5 1.408, 95% CI: 1.128–1.757; P = 0.002) (Xuan et al., 2014). Hishida et al. (2013) found a significant association between the subject with the T/T genotype of MTHFR A222V polymorphism and the elevated risk of chronic kidney disease (CKD) using cross-sectional data. When those with MTHFR A222V C/C were defined as references, those with MTHFR A222V C/T and T/T demonstrated the aORs for CKD of 1.14 (95 % CI 0.93–1.40) and 1.39 (1.06–1.82), respectively (Hishida et al., 2013). A meta-analysis of 50 studies (19,260 cases and 26,364 controls) presented by He and Shen (2017) showed that there was a significant association between MTHFR A222V polymorphism with breast cancer risk for homozygous model (OR =1.20, CI: 1.12– 1.28, P<0.05), recessive model (OR =1.19, CI: 1.11– 1.27, P<0.05), dominant model (OR =1.19, CI: 1.79–1.95, P<0.05), and allele model (OR =1.19, CI: 1.12–1.28, P<0.05) from overall population analysis and with ovarian cancer risk for homozygous model (OR =1.43, CI: 1.30–1.59, P<0.05), recessive model (OR =1.35, CI: 1.23– 1.48, P<0.05), dominant model (OR =1.20, CI: 1.12–1.28, P<0.05), and allele (OR =1.19, CI: 1.13–1.25, P<0.05) only in Asians (He and Shen, 2017). In addition, Li et al (2021) suggested that the MTHFR A222V polymorphism is significantly associated with susceptibility to triple-negative breast cancer from 980 southwestern China subjects (490 breast cancer cases and 490 cancer-free controls) (OR = 2.83, 95% CI: 1.12– 9.51, P = 0.0401)(Li et al., 2021). On the other hand, a meta-analysis by Kumar et al. (2018), of 29 studies with 6520 cases and 9192 controls included resulting that there were significant Research Article Journal of Pharmaceutical Sciences and Community Molecular Dynamics Study of Human MTHFR ... 88 Yanuar et al. J. Pharm. Sci. Community, 2023, 20(1), 84-89 associations between A222V polymorphism and esophageal cancer risk using overall comparisons in five genetic models (T vs. C: OR = 1.20, 95% CI = 1.1–1.27, p = 0.05) (Kumar and Rai, 2018). In their meta-analysis, Zhang et al (2022) found that MTHFR A222V polymorphism is significantly related to schizophrenia and major depression in the overall population. MTHFR A222V has been linked to an increased risk of bipolar disorder in the recessive model (TT vs. CT + CC). Ethnic subgroup analysis shows that schizophrenia for T vs. C: OR = 1.19, 95% CI = 1.11–1.29, P < 0.001; for TT + CT vs. CC: OR = 1.22, 95% CI = 1.10–1.35, P < 0.001; for TT vs. CT + CC: OR = 1.31, 95% CI = 1.16–1.48, P < 0.001; for TT vs. CC: OR = 1.46, 95% CI = 1.24–1.72, P < 0.001 and major depression for T vs. C: OR = 1.46, 95% CI = 1.21–1.77, P < 0.001; for TT + CT vs. CC: OR = 1.52, 95% CI = 1.11–2.08, P = 0.009; for TT vs. CT + CC: OR = 1.75, 95% CI = 1.34– 2.28, P < 0.001; for TT vs. CC: OR = 1.89, 95% CI = 1.40– 2.57, P < 0.001 significantly correlate with MTHFR A222V in Asian populations (Zhang et al., 2022). This concludes our finding that MTHFR A222V mutant protein is considered unstable by molecular dynamics study in silico which may lead to various clinical problems for the person possessing this mutation as described above from other studies. This study had several limitations. We did not measure the average value of the RMSD ligand movement deviation (∆RMSDligmove) in every 5 ns and MM/PBSA analysis. Therefore, we suggest future studies to examine the protein stability via ∆RMSDligmove and MM/PBSA analysis to give more results to be analyzed so it’s described the real MTHFR wildtype and MTHFR A222V molecular dynamics. CONCLUSIONS The present study revealed that MTHFR wild-type protein is stable and MTHFR A222V mutant protein is unstable during MD simulations in silico. Further study is needed to examine ∆RMSDligmove and MM/PBSA analysis for more precise results. ACKNOWLEDGEMENTS This research was supported and funded by Ministry of Education, Culture, Research, and Technology (1989.9/LL5-INT/PG.02.00/2022, 028 Penel.LPPM-USD/V/2022) and Institute for Research and Community Services, Universitas Sanata Dharma (Research Grant Contract No. 007 Penel./LPPM-USD/II/2022). CONFLICT OF INTEREST The authors declare no conflict of interest in this study. REFERENCES Abhinand, P.A., Shaikh, F., Bhakat, S., Radadiya, A., Bhaskar, L.V.K.S., Shah, A., Ragunath, P.K., 2016. Insights on the Structural Perturbations in Human MTHFR Ala222Val Mutant by Protein Modeling and Molecular Dynamics. Journal of Biomolecular Structure and Dynamics, 34(4): 892–905. Castro, R., Rivera, I., Ravasco, P., Jakobs, C., Blom, H.J., Camilo, M.E., de Almeida, I.T., 2003. 5,10-Methylenetetrahydrofolate Reductase 677C→T and 1298A→ C Mutations are Genetic Determinants of Elevated Homocysteine. QJM - Monthly Journal of the Association of Physicians, 96(4): 297–303. Chehadeh, E.H., S.W., Jelinek, H.F., Al Mahmeed, W.A., Tay, G.K., Odama, U.O., Elghazali, G.E.B., Al Safar, H.S., 2016. Relationship Between MTHFR C677T and A1298C Gene Polymorphisms and Complications of Type 2 Diabetes Mellitus in an Emirati Population. Meta Gene, 9: 70–75. Froese, D.S., Kopec, J., Rembeza, E., Bezerra, G.A., Oberholzer, A.E., Suormala, T., et al., 2018. Structural Basis for the Regulation of Human 5,10-Methylenetetrahydrofolate Reductase by Phosphorylation and S- adenosylmethionine Inhibition. Nature Communications, 9(1): 1–13. He, L., Shen, Y., 2017. MTHFR C677T Polymorphism and Breast, Ovarian Cancer Risk: a meta-analysis of 19,260 patients. OncoTargets and Therapy, 10: 227–238. Hishida, A., Okada, R., Guang, Y., Naito, M., Wakai, K., Hosono, S., Nakamura, K., Turin, T.C., Suzuki, S., Niimura, H., Mikami, H., et al., 2013. MTHFR, MTR and MTRR Polymorphisms and Risk of Chronic Kidney Disease in Japanese: Cross- Sectional Data from the J-MICC Study. International Urology and Nephrology, 45(6): 1613–1620. Istyastono, E.P., Prasasty, V.D., 2021. Computer- Aided Discovery of Pentapeptide AEYTR as a Potent Acetylcholinesterase Inhibitor. Indonesian Journal of Chemistry, 21(1): 243–250. Istyastono, E.P., Riswanto, F.D.O., 2022. Molecular Dynamics Simulations of the Caffeic Acid Interactions to Dipeptidyl Journal of Pharmaceutical Sciences and Community Molecular Dynamics Study of Human MTHFR... Research Article 89 Yanuar et al. J. Pharm. Sci. Community, 2023, 20(1), 84-89 Peptidase IV. International Journal of Applied Pharmaceutics, 14(4): 274–278. Krieger, E., Vriend, G., 2015. New Ways to Boost Molecular Dynamics Simulations. Journal of Computational Chemistry, 36(13): 996– 1007. Kumar, P., Rai, V., 2018. MTHFR C677T Polymorphism and Risk of Esophageal Cancer: An Updated Meta-Analysis. Egyptian Journal of Medical Human Genetics, 19(4): 273–284. Li, Z., Zhang, J., Zou, W., Xu, Q., Li, S., Wu, J., Zhu, L., Zhang, Yunjiao, et al., 2021. The Methylenetetrahydrofolate Reductase (MTHFR) C677T Gene Polymorphism is Associated with Breast Cancer Subtype Susceptibility in Southwestern China. PLoS ONE, 16(7 July): 1–12. Liu, K., Watanabe, E., Kokubo, H., 2017. Exploring the Stability of Ligand Binding Modes to Proteins by Molecular Dynamics Simulations. Journal of Computer-Aided Molecular Design, 31(2): 201–211. NCBI, 2022. MTHFR methylenetetrahydrofolate reductase [Homo sapiens (human)]. URL https://www.ncbi.nlm.nih.gov/gene/452 4 (accessed 11.29.22). Raghubeer, S., Matsha E., T., 2021. Methylenetetrahydrofolate (MTHFR), the One-Carbon Cycle, and Cardiovascular Risks. Nutrients, 13: 4562. Ueland, P.M., Hustad, S., Schneede, J., Refsum, H., Vollset, S.E., 2001. Biological and Clinical Implications of the MTHFR C677T Polymorphism. Trends Pharmacol. Sci., 22(4): 195–201. Wu, S., Han, Y., Hu, Q., Zhang, X., Cui, G., Li, Z., Guan, Y., 2017. Effects of Common Polymorphisms in the MTHFR and ACE Genes on Diabetic Peripheral Neuropathy Progression: a Meta-Analysis. Molecular Neurobiology, 54(4): 2435–2444. Xuan, C., Li, H., Zhao, J.X., Wang, H.W., Wang, Y., Ning, C.P., Liu, Z., Zhang, B.B., et al., 2014. Association between MTHFR polymorphisms and congenital heart disease: A meta-analysis based on 9,329 cases and 15,076 controls. Scientific Reports, 4: 1–13. Yigit, S., Karakus, N., Inanir, A., 2013. Association of MTHFR Gene C677T Mutation with Diabetic Peripheral Neuropathy and Diabetic Retinopathy. Molecular Vision, 19(July): 1626–1630. Zhang, Y.X., Yang, L.P., Gai, C., Cheng, C.C., Guo, Z.Y., Sun, H.M., Hu, D., 2022. Association Between Variants of MTHFR Genes and Psychiatric Disorders: A Meta-Analysis. Frontiers in Psychiatry, 13(August): 1–19.