Molecular Modeling and Molecular Dynamics Simulation Studies on the Selective Binding Mechanism of MTHFD2 Inhibitors

Methylenetetrahydrofolate dehydrogenase 2 (MTHFD2) is a mitochondrial enzyme that plays an important role in purinecarbon metabolism and thymidine biosynthesis. It has attracted broad interest as a novel therapeutic target for cancer. However, a major problem of current MTHFD2 inhibitors is their lack of selectivity and reactivity with its closest isoform, MTHFD1. Recently, the first selective MTHFD2 inhibitor, DS44960156, has been reported and it exhibits a more than 18-fold selectivity for MTHFD2 over MTHFD1. However, mechanism of DS44960156 selective binding to MTHFD2 over MTHFD1 is unknown. In this study, molecular docking, molecular dynamics (MD) simulations, molecular mechanics generalized born/surface area (MM_GBSA) binding free energy calculations, and analysis of the decomposition of binding free energies were used to investigate the selective binding mechanism of DS44960156 to the folate-binding site of MTHFD2 over MTHFD1. The results revealed that contributions from residues Gln100/Gln132, Val55/Asn87, and Gly237/Gly310 in the binding pocket of MTHFD1/MTHFD2 are the key factors responsible for the binding selectivity. These findings explain the selectivity of DS44960156 to MTHFD2 over MTHFD1, and may provide guidance for the future study and design of novel MTHFD2 inhibitors.


Introduction
Anticancer drugs such as alkylating agents, topoisomerase inhibitors, and anti-microtubule agents directly target DNA replication and cell division; however, How to cite this paper: Qu [1]. MTHFD2 is a mitochondrial enzyme that plays an important role in purine carbon metabolism and thymidine biosynthesis. Sex-determining region Y-related high mobility group box 7, SOX7, as a tumor suppressor, can combine with the MTHFD2 promoter to inhibit the expression of MTHFD2.
Therefore, MTHFD2 is not expressed in most healthy adults [2]. However, the induction of KRAS can activate the AKT and ERK1/2 pathways, and recognize and bind to the MTHFD2 promoter through c-Myc to up-regulate the expression of MTHFD2. High expression of MTHFD2 mRNA and protein has been found in patients with various types of cancer [3] [4] [5]. Thus, MTHFD2 inhibitors are considered as potential cancer treatment drugs with minimal side effects. A few of MTHFD2 inhibitors, such as LY345899 [6] [7] and carolacton [8] have been identified. However, an important problem of these inhibitors is the lack of selectivity for MTHFD2 over MTHFD1, which is because of the similar structure of the 2 enzymes. As shown in Figure 1, MTHFD1 shares 80% similarity with MTHFD2 in the folate-binding site. Inhibition of MTHFD1 is considered to be a potential safety risk, because it is broadly expressed in normal tissues [9].
Therefore, selective inhibitors of MTHFD2 are highly attractive for cancer treatment and drug development.
Recently, Kawai et al. reported the first selective inhibitor of MTHFD2, DS44960156 (Figure 2(a)), which exhibits a more than 18-fold selectivity for MTHFD2 over MTHFD1 [9]. X-ray analysis of the crystal structure of the MTHFD2-DS44960156 complex shows that DS44960156 occupies the folate-binding site of MTHFD2 (Figure 2(b)). The key interactions of the MTHFD2-DS44960156 complex are: 1) Four hydrogen bonds: Gln132 and Lys88   in MTHFD2 with C=O of the pyrimidin-4-one in DS44960156; Asn87 in MTHFD2 with C=O of the linker amide in DS44960156; Gly310 in MTHFD2 with S=O of the sultam in DS44960156; 2) Anda π-π interaction between Tyr84 in MTHFD2 and the pyrimidin-4-one in DS44960156. However, because of the highly similar structures between the folate-binding sites of MTHFD1 and MTHFD2, it is unclear why DS44960156 selectively binds to MTHFD2 over MTHFD1.
In this study, molecular modeling and molecular dynamics (MD) simulations were performed to investigate the selective binding mechanism of DS44960156 to MTHFD2 over MTHFD1. We performed MD simulations on the complex structures of DS44960156-MTHFD1 and DS44960156-MTHFD2. The conformation of DS44960156 in the folate-binding site of MTHFD1 and MTHFD2 was compared. Binding free energies were calculated to identify the key residues responsible for the different binding affinity. These results illustrate the selectivity of DS44960156, and may provide guidance for design of novel MTHFD2 inhibitors.

System Preparations
The crystal structure of MTHFD2-DS44960156 (PDB code 6JIB) [9] was downloaded from the RCSB Protein Data Bank, and the structure of MTHFD1-DS44960156 was predicted by molecular docking. Docking of DS44960156 to the folate-binding site of MTHFD1 was performed using the AutoDock Vina package [10], and the crystal structure of MTHFD1-LY345899 (PDB code 6ECQ) [11] with LY345899 deleted was used for the docking. All non-bonded hetero-atoms and water molecules were removed fromMTHFD1, polar hydrogen atoms and subsequently the Kollman united atom partial charges were added and assigned to MTHFD1. All hydrogen atoms, including polar and nonpolar atoms, and the Gasteiger-Hückel atomic charge were added and assigned to DS44960156. The MTHFD1 receptor was treated as rigid, and DS44960156 was treated as flexible during docking. We created a 40 × 40 × 40 Å 3 rectangular parallelepiped around the cavity, ensuring that it is the right size to allow the ligand to rotate in it. The best pose of the complex, according to the criteria of lowest energy, was selected for MD simulation steps.

MD Simulations
MD simulations of the complex structures of DS44960156 with MTHFD2 and MTHFD1 were performed. The general AMBER force field (GAFF) [12] was used to parameterize the inhibitor DS44960156 and NADP. The electrostatic potential of DS44960156 was calculated at the HF/6-31G* level using the program Gaussian 09. Afterward, the missing force field parameters were generated using the antechamber module in AMBER 20 [13]. The AMBER ff14SB force field was used to parameterize MTHFD1 and MTHFD2. The complexes were solvated in TIP3P [14] water boxes with an internal offset distance of 10.0 Å. Calcium and chlorine ions were added to neutralize the system.
All of the MD simulations were performed using GROMACS-5.0.7 software [15]. Firstly, the systems were subjected to energy minimizations in 2 steps [16].
The solvent water molecules were optimized first, while the coordinates of other molecules were fixed, followed by an unrestrained minimization of all the atoms in the system. Secondly, the system was gradually heated to 300 K under constant volume. Thirdly, a 5 ns simulation under constant volume and temperature was performed with the non-hydrogen atoms of the protein fixed, followed by a 5 ns simulation under constant pressure and temperature. Finally, 300 ns production simulation was performed. The temperature was maintained at 300 K with a coupling coefficient of 1.0 ps through the V-rescale method. The particle-mesh Ewald method [17] was used to treat long-range electrostatic interactions with a real space cutoff of 1 nm and the cutoff of the van der Waals interactions was set to 1 nm. Bonds involving hydrogen atoms were constrained using the LINCS algorithm [18]. The time step was set to 2 fs and the periodic boundary conditions were applied in all directions.

MM/GBSA Free Energy Calculations and Decomposition
The binding free energy ( bind G ∆ ) between DS44960156 and MTHFD1/2 was calculated by the MM/GBSA method: [19] [20] [21] [22] ( ) where complex G , receptor G , and ligand G are the free energies of the complex, receptor, and ligand, respectively; MM E ∆ is the gas-phase interaction energy between the receptor and the ligand; ∆ is calculated using the GB model [24]. The binding free energies were then decomposed into contributions of protein-ligand interaction pairs using where , represents the fluctuation of protein-ligand interaction energy around the average energy, and … stands for a statistical average.

Molecular Docking
In this study, the complex structure of MTHFD1-DS44960156 was predicted by molecular docking. To examine the prediction accuracy, we first performed re-docking of DS44960156 into MTHFD2. The re-docked conformation was then superimposed onto the crystal conformation (Figure 3

Structural Stability
MD simulations can accurately identify correct binding conformations, and provide valuable information of the dynamic properties of receptor-ligand interactions [27]. We then performed 300 ns MD simulations on the crystal structure of MTHFD2-DS44960156 and the docked structure of MTHFD1-DS44960156. The RMSD for all the C α atoms of the 2 proteins relative to the initial structure was  Figure 4(a), both systems reached stability from the beginning of the simulations. To further understand the stability of DS44960156 in the folate-binding sites of MTHFD2 and MTHFD1, the RMSD of heavy atoms of DS44960156 was calculated. As shown in Figure 4(b), DS44960156 exhibited significantly larger conformational dynamics in MTHFD1 than in MTHFD2.
These results indicate that DS44960156 is more stable in the folate-binding site of MTHFD2 which may be due to its high binding affinity. The difference of the RMSDs between DS44960156 in the folate-binding sites of MTHFD1 and MTHFD2 indicates that DS44960156 may have a different binding mode in MTHFD1 and MTHFD2. To further understand this finding, we superimposed the frames of the MD simulations of MTHFD1-DS44960156 and MTHFD2-DS44960156. As shown in Figure 5, DS44960156 occupied the folate-binding site ofMTHFD1 as well as MTHFD2, with the terminal benzene ring points out of the pocket, and carboxyl group points into the pocket. However, DS44960156 in MTHFD1 is farther from the folate-binding site than in MTHFD2, which may cause the different binding affinity, thus the difference of the RMSDs.

Binding Free Energy
To further understand the binding affinity between DS44960156 and MTHFD1/2, MM_GBSA binding free energy calculations were performed on the last 50 ns MD trajectories. As shown in Table 1

Residue-Specific Binding Free Energies
In order to explore the reasons for the stronger binding affinity to MTHFD2, the decomposition of bind G ∆ for the MTHFD1-DS44960156 and MTHFD2-DS44960156 systems was performed using the MM_GBSA method. Residues that contributed more than 5 kJ/mol were defined as hot-spot residues. As shown in Table 2, all of the residues in the binding pocket of MTHFD2 were hot-spot residues. However, there were only 2 hot-spot residues (Tyr 52 and Lys 56) in MTHFD1. These results indicate that the contributions from residues Gln100/Gln132, Val55/Asn87, and Gly237/Gly310 in the binding pockets of MTHFD1/2 are the key factors responsible for the binding selectivity of DS44960156.
The positions of these residues in the 2 complexes are shown in Figure 6. It  can be seen that MTHFD2 forms 4 hydrogen bonds (Gln132/Lys88 with C=O of the pyrimidin-4-one; Asn87 with C=O of the linker amide; Gly310 with S=O of the sultam), and a strong π-π interaction (Tyr84 and pyrimidin-4-one) with DS44960156. However, in the MTHFD1-DS44960156 complex, DS44960156 only forms a hydrogen bond and a weak π-π interaction with Lys56 and Tyr52 in MTHFD1. The potential reasons for these interaction differences may be that residue Asn 87 is replaced by a non-polar residue (Val 55) in MTHFD1 which cannot interact with the C=O of the linker amide in DS44960156. This causes the position of DS44960156 in MTHFD1 to be farther from the folate-binding site than in MTHFD2. Therefore, the C=O of pyrimidin-4-one in DS44960156 cannot interact with MTHFD1, and the π-π interaction between Tyr84 and DS44960156 is weak. These interaction differences may be responsible for the difference in selectivity.

Conclusion
In this study, molecular docking, MD simulations, MM_GBSA binding free energy calculations, and calculations of the decomposition of binding free energies were performed to explore the selective binding mechanism of DS44960156 to the folate-binding site of MTHFD2 over that of MTHFD1. MM_GBSA binding free energy calculations showed that bind G ∆ was 52.7 kJ/mol lower in the MTHFD1-DS44960156 complex compared with the MTHFD2-DS44960156 complex. The per-residue free energy decomposition shows that all the residues in the binding pocket of MTHFD2 have strong interactions with DS44960156. Nonetheless, the interactions between residues Gln100, Val55 and Gly237 in the binding pocket of MTHFD1 and DS44960156 are quite weak. These results indicate that the contributions from residues Gln100/Gln132, Val55/Asn87, and Gly237/ Gly310 in the folate binding pockets of MTHFD1 and MTHFD2 are the key factors responsible for the binding selectivity of DS44960156. These results might be useful for designing novel selective MTHFD2 inhibitors.

Conflicts of Interest
The author declares no conflicts of interest regarding the publication of this paper.