dUMP/F-dUMP Binding to Thymidylate Synthase: Human Versus Mycobacterium tuberculosis.

Thymidylate synthase is an enzyme that catalyzes deoxythymidine monophosphate (dTMP) synthesis from substrate deoxyuridine monophosphate (dUMP). Thymidylate synthase of Mycobacterium tuberculosis (MtbThyX) is structurally distinct from its human analogue human thymidylate synthase (hThyA), thus drawing attention as an attractive drug target for combating tuberculosis. Fluorodeoxyuridylate (F-dUMP) is a successful inhibitor of both MtbThyX and hThyA, thus limited by poor selectivity. Understanding the dynamics and energetics associated with substrate/inhibitor binding to thymidylate synthase in atomic details remains a fundamental unsolved problem, which is necessary for a new selective inhibitor design. Structural studies of MtbThyX and hThyA bound substrate/inhibitor complexes not only revealed the extensive specific interaction network between protein and ligands but also opened up the possibility of directly computing the energetics of the substrate versus inhibitor recognition. Using experimentally determined structures as a template, we report extensive computer simulations (∼4.5 μs) that allow us to quantitatively estimate ligand selectivity (dUMP vs F-dUMP) by MtbThyX and hThyA. We show that MtbThyX prefers deprotonated dUMP (enolate form) as the substrate, whereas hThyA binds to the keto form of dUMP. Computed energetics clearly show that MtbThyX is less selective between dUMP and F-dUMP, favoring the latter, relative to hThyA. The simulations reveal the role of tyrosine at position 135 (Y135) of hThyA in amplifying the selectivity. The protonation state of the pyrimidine base of the ligand (i.e., keto or enolate) seems to have no role in MtbThyX ligand selectivity. A molecular gate (consists of Y108, K165, H203, and a water molecule) restricts water accessibility and offers a desolvated dry ligand-binding pocket for MtbThyX. The ligand-binding pocket of hThyA is relatively wet and exposed to bulk water.

Introduction

Tuberculosis is a leading infectious disease caused by its etiological agent Mycobacterium tuberculosis (Mtb) inducing ∼1.5 million deaths and ∼10 million reported cases per year.[1] The attempts to eradicate this disease have been hindered by the emergence of multidrug-resistant and extremely drug-resistant strains of Mtb.[2] This resistance to the existing drugs urges the necessity of designing new drugs and discovering novel drug targets. One such target is thymidylate synthase of M. tuberculosis (MtbThyX) which is coded by thyX gene and is structurally distinctive from its counterpart in human (hThyA, coded by thyA gene).[3] Both MtbThyX and hThyA catalyzes the formation of 2′-deoxythymidine 5′-monophosphate (dTMP) from deoxyuridine monophosphate (dUMP) using 5,10-methylenetetrahydrofolate (mTHF) as a carbon-donor cofactor.[4] dTMP is an important precursor of DNA, and its inhibition inflicts a cascade of events leading to “thymine-less” death.[5,6] The mechanism of dUMP → dTMP conversion is also proposed to be different for MtbThyX and hThyA.[7] Unlike hThyA, MtbThyX involves a flavin adenine dinucleotide (FAD) molecule which is sandwiched between dUMP and mTHF during the catalysis.[7,8] Thus, MtbThyX is also called flavin-dependent thymidylate synthase. In Mtb, two types of thymidylate synthase (ThyX and ThyA) are present, but experiments suggest that only ThyX is essential for pathogenic growth, thus making MtbThyX a pharmacologically attractive target.[9,10] Substitution into native substrate dUMP with fluorine is a successful inhibitor (F-dUMP) of thymidylate synthase.[11] F-dUMP is a metabolite form of the anti-cancerous drug 5-fluorouracil (5-FU) and is known to inhibit both ThyX and ThyA.[12,13] Thus, the application of F-dUMP as a drug for M. tuberculosis is limited by poor selectivity.

MtbThyX is a homo-tetramer containing four ligand-binding sites (Figure a). Each ligand-binding site is at the intersection of three monomeric subunits, and the ligand interacts with two subunits (Figure a).[14] The ligand-binding site is deeply buried inside the protein. Note, MtbThyX bound to native dUMP substrate is not known till date. The X-ray structure of F-dUMP bound MtbThyX reveals the intricate interaction network which includes two characteristic features: (A) F-dUMP binding pocket and (B) lid above the Hoogsteen edge of F-dUMP (Figure a).

Figure 1

(a) X-ray structure of MtbThyX (homotetramer; monomeric units are in yellow, cyan, green and purple, PDB 3GWC(16)). Each ligand-binding site (out of four) is at the intersection of three monomeric units. Zoomed in view of (i) binding site and (ii) lid. CW1 and CW2 (crystal waters, red sphere), cofactor (FAD, cyan-line), and ligand (F-dUMP, yellow). (b) X-ray structure of hThyA (Homodimer; monomeric units are in yellow and green). Zoomed in view of substrate (dUMP) binding pocket with antifolate (MTX, yellow-line). Key residues (in sticks) and interaction network shown in black broken lines. Watson–Crick edge (the red surface) is labeled with an arrow; the opposite side of the pyrimidine ring is the Hoogsteen edge.

(A) F-dUMP binding pocket: Watson–Crick edge of the pyrimidine base forms direct interaction with R199, R107, and crystal-water (CW1).[15] The proximity of R199 and the Watson–crick edge of pyrimidine (distance between N3 of F-dUMP and guanidinium nitrogen of R199 ∼ 2.8 Å)[15] and the conserved crystal-water (CW1) were isolated in the structures of ThyX from various species.[16−19] The phosphate moiety of the F-dUMP is surrounded by side chains of R87, R172, S105, and E92.[16] The nature of amino-acid at 92nd position is also highly conserved, either as glutamic acid (E92) or aspartic acid (D92).20F-dUMP in the MtbThyX binding pocket is highly compact, where the phosphate moiety is placed very close to the Hoogsteen edge of its pyrimidine base (Figure a). The presence of phosphate close to the Hoogsteen edge of the ligand has been proposed to facilitate methyl-transfer reaction.[17] (B) Lid over F-dUMP: Interactions between residues Y108, K165, and the backbone of H203 with water (CW2, crystal water) at the center result in the formation of a lid above the Hoogsteen edge of the F-dUMP ligand. The lid could be thought of as the lining of the ligand-binding pocket.

ThyA is a well-studied classical thymidylate synthase, which is a homo-dimeric protein.[21] Each monomer of ThyA binds to one dUMP ligand (two dUMP per ThyA protein) resulting in the formation of the binary complex, followed by cofactor (mTHF) binding resulting in the formation of a ternary complex.[8] The ternary complex of hThyA bound to dUMP and antifolate (methotrexate: MTX) is shown in Figure b. Watson–Crick edge of dUMP is in interaction with the side-chain of N226 and backbone NH of D218 (Figure b).[8,22] Residue Y135 was proposed to act as a base and plays a key role in the methyl-transfer reaction.[23] Structural studies of thymidylate synthase (ThyA and ThyX) with various ligands (e.g., dUMP, F-dUMP, Br-dUMP, F-U, etc.) revealed the intricate interaction network involved in the structural integrity of those complexes.[16,20,24] Biochemical studies showed that catalytic activity (kcat) of MtbThyX is surprisingly low compared to hThyA, whereas, the substrate (dUMP) binding affinity (Km) is more or less similar.[12]

Despite structural and biochemical studies, several key issues related to ligand-binding to thymidylate synthase require attention. Direct link between 3D structure and energetics (of ligand selectivity: dUMP versus FdUMP) of thymidylate synthase is unknown. Further, the proximity of R199 to the pyrimidine base of F-dUMP in the X-ray structure of MtbThyX (Figure a) gives rise to several possible models of the ligand in the MtbThyX binding pocket [viz., Neutral(R199): keto(F-dUMP), Charged (R199): enol (F-dUMP), and charged (R199): deprotonated/enolate (F-dUMP–1) pairs]. Thus tautomerization (keto, enol) and protonation state (neutral, protonated, deprotonated) of the ligand is certainly of considerable importance for MtbThyX and has not been addressed in previous studies. The intrinsic stabilities of the nature of the ligand (keto, enol, enolate) in the protein binding pocket is not known. Does the nature of the ligand alter the MtbThyX selectivity (dUMP vs F-dUMP)? What is the function of the lid (Figure a) in MtbThyX?

The near-atomic resolution (resolution ∼ 2 Å) X-ray structures of MtbThyX and hThyA are sufficiently accurate for performing a structure-based computational evaluation of the dynamics and energetics associated with ligand-binding. Structure-based molecular dynamics (MD) simulations were performed previously to study the dynamics and stability of human ThyA in complexes with various ligands.[25] However, MtbThyX complexes have not been investigated using MD simulations. Using the available X-ray structures as a template, we report computer simulations of MtbThyX and hThyA to decipher the energetics of ligand recognition (dUMP vs F-dUMP), thereby comparing the dynamics and energetics. This paper is focused on addressing the above-mentioned issues that are relevant to ligand-binding to MtbThyX and hThyA.

Methods

MD Setup

Initial structures of MtbThyX (bound to F-dUMP and FAD) and hThyA [bound to dUMP and MTX (antifolate)] were obtained from the protein data bank entry 3GWC (resolution 1.9 Å) and 5X66 (resolution 1.99 Å), respectively.[16,22] Models of dUMP bound MtbThyX and F-dUMP bound hThyA were obtained by atomic replacement (“F” by “H” vice versa) at C-5 carbon. Because MtbThyX is a large protein (88 × 86 × 65 Å3) with four ligand-binding sites, we considered a spherically truncated protein model centered at the ligand. Our truncated simulation model included protein residues having at least one non-hydrogen atom within a 25 Å radius sphere, centered at the N1-atom of the nitrogenous base of the ligand (dUMP or F-dUMP). The objective of this work is to probe the local structural fluctuations of the ligand-binding site in response to different forms (keto, enol, and enolate) of the ligands (dUMP/F-dUMP), rather than sampling distant larger scale conformational changes that certainly demand longer simulation time scales for convergence. Harmonic restraints were imposed on the protein-heavy atoms between 22 and 25 Å region (buffer region). Harmonic restraints were gradually incremented toward the solute boundary (22–23 Å with force constant = 3 kcal mol–1 Å–2, 23–24 Å with force constant = 4 kcal mol–1 Å–2, 24–25 Å with force constant = 5 kcal mol–1 Å–2, >25 Å region with force constant = 10 kcal mol–1 Å–2). Solute atoms within 22 Å were fully flexible in the MD trajectory, except FAD molecule which was harmonically retrained with force constant = 10 kcal mol–1 Å–2 to its experimentally determined structure.

We have considered full human ThyA dimer without truncation for MD simulations. A model of the binary complex (without MTX) was obtained from PDB 5X5D.[22] Both hThyA and truncated MtbThyX models were overlaid with a cubic water box of 90 Å edge. The total charge of the simulation box was neutralized by adding counter ions Na+/Cl–. Two Cl– ions were added in MtbThyX, and ten Na+ ions were added in hThyA for overall charge neutralization. MD simulation models consist of 70,519 and 70,208 atoms for MtbThyX and hThyA complexes, respectively.

MD simulations were performed under periodic boundary condition using the particle mesh Ewald method with tinfoil boundary condition for long-range electrostatics.[26−28] A cutoff of 16 Å was considered for calculating the Van Der Waals interactions. Temperature and pressure were maintained at 310 K and 1 bar for performing simulations considering NPT ensemble. The temperature was controlled using Langevin dynamics (coupling coefficient of 5 ps–1), whereas pressure was controlled by a Langevin piston.[29,30] Standard CHARMM36 force-field has been used for the protein, water, and ions.[31,32] TIP3P model was used to model waters. Non-standard ligand (FAD, MTX, F-dUMP(k), F-dUMP(d), F-dUMP(e), dUMP(d)) parameters were obtained from CGenFF.[33−36] CHARMM and NAMD software were used for system preparation and calculations.[37,38] For adequate sampling, every system was run five times with different initial velocities and minimization steps. Minimization was done using the algorithm steepest descent for 100 steps followed by adopted basis Newton–Raphson algorithm for 100 steps. The simulation involves 300 ps of equilibration in which the system is heated up to 310 K, and the harmonic restraints on the heavy atoms were gradually released for the 22 Å inner region of MtbThyX. No restraint was applied for human ThyA. After equilibration, each replica was subjected to 50 ns production dynamics for structural analysis and energy calculations.

The size of the simulation system is crucial, and the distant region of the biomolecule may contribute significantly to the ligand-binding. As a check, we also considered the simulation model by including full tetrameric MtbThyX protein with bound ligand (dUMP) or inhibitor (F-dUMP) and subjected to MD simulation in a different software GROMACS.[39] The full MtbThyX protein was overlaid in a 115 Å edge water-box. The whole system (153,442 atoms including twelve Na+ ions) was subjected to minimization for 1800–2000 steps (steepest descent algorithm) followed by equilibration of 1 ns in NPT ensemble, where temperature (310 K) and pressure (1 bar) were controlled by Berendsen thermostat and Parrinello-Rahman algorithm.[40,41] No restraint was applied for full MtbThyX, and five independent trajectories were generated by varying initial velocities. Last 50 ns of each post-equilibrated trajectory (out of 5) of each simulation model was considered for structural analysis and binding free energy calculations. The MD structures of the ligand-binding pocket, obtained from full MtbThyX, are more or less identical to the truncated 25 Å model (Tables S1 and S2). The same MD protocol described for full MtbThyX simulation was used to study the dynamics of hThyA bound to nolatrexate (antifolate) (PDB 5X67) in GROMACS. This nolatrexate-bound hThyA system overlaid with a smaller water-box of edge 92 Å, and the MD model consists of total 77,803 atoms. Structure and energetics obtained from simulations of MTX bound hThyA are more or less identical to nolatrexate bound ThyA (Tables S3 and S4). For enolate dUMP/F-dUMP, we have also adopted an ad hoc parameterization (distributing −1 charge over the Watson crick edge of uridine) approach to create the force-field parameters. MD structures obtained with ad hoc enolate parameters also reproduced the X-ray structures (results not shown). Results obtained from different MD setup using different simulation protocols (including two different software) confirmed the robustness of the simulated results. The uncertainty in the averaged ΔGbind (∼1.0 kcal/mol) is well within the acceptable statistical uncertainty. Good structural and energetic convergence (Tables S1–S4 and 1) was obtained by performing a total of ∼4.5 μs of simulations using two different software: NAMD (2.4 and 0.6 μs for MtbThyX and hThyA, respectively) and GROMACS (1 and 0.5 μs for MtbThyX and hThyA, respectively).

Table 1

MD Trajectory Averaged MM/PBSA Ligand-Binding Free Energy (ΔGbind) and Its Components: ΔGMM (Molecular Mechanics), ΔGPB (Polar Solvation), and ΔGSA (Non-Polar Solvation)a

	system	model	ΔGMM	ΔGPB	ΔGSA	ΔGbind	ΔΔGbind
CHARMM	MtbThyX	dUMP(d)	–370.73 ± 1.10	236.56 ± 0.40	–4.00 ± 0.005	–138.18 ± 1.08	2.69±1.65
		F-dUMP(d)	–376.33 ± 1.21	239.46 ± 0.44	–4.00 ± 0.005	–140.87 ± 1.26
		dUMP(k)	–300.51 ± 0.85	221.00 ± 0.30	–3.88 ± 0.004	–83.95 ± 0.97	1.8±1.55
		F-dUMP(k)	–305.83 ± 0.96	225.51 ± 0.62	–4.09 ± 0.004	–85.75 ± 1.21
	hThyA	dUMP(k)	–294.22 ± 1.02	221.03 ± 0.63	–3.84 ± 0.005	–77.10 ± 1.25	6.75±1.51
		F-dUMP(k)	–306.07 ± 0.66	226.05 ± 0.46	–3.84 ± 0.005	–83.85 ± 0.86
GROMACS	MtbThyX	dUMP(d)	–489.25 ± 1.01	344.63 ± 0.51	–3.06 ± 0.005	–147.68 ± 0.59	2.38±0.75
		F-dUMP(d)	–493.04 ± 0.85	346.11 ± 0.44	–3.14 ± 0.005	–150.06 ± 0.47
		dUMP(k)	–292.41 ± 0.67	243.31 ± 0.39	–3.22 ± 0.006	–52.34 ± 0.38	3.84±0.55
		F-dUMP(k)	–295.60 ± 0.83	242.64 ± 0.53	–3.22 ± 0.007	–56.18 ± 0.40
	hThyA	dUMP(k)	–394.68 ± 2.23	302.81 ± 0.85	–3.84 ± 0.005	–95.34 ± 1.12	6.21±1.82
		F-dUMP(k)	–414.84 ± 2.01	316.90 ± 1.10	–3.84 ± 0.005	–101.55 ± 1.44

The standard error of mean (s.e.m) was reported as error associated with free energy differences (ΔGbind, ΔGMM, ΔGPB, and ΔGSA). Free energy estimation was done by performing MD simulations of two different MD models using two different software (CHARMM and GROMACS; see Methods). Relative binding free energy was estimated as ΔΔGbind = ΔGbind(dUMP) – ΔGbind(F-dUMP), and error was calculated by propagating the standard error of the mean associated with the averaged ΔGbind. One hundred evenly spaced snapshots of last 25 ns MD trajectory were considered for free energy estimate (in kcal/mol).

Free Energy Calculation by a MMPBSA Method

Protein and ligand-binding can be described with an equation

Binding affinity should in principle be calculated from the absolute free energy difference of the complex (GC), free protein (GP), and free ligand (GL) using the formulawhere ⟨ ⟩ represent averaging over the three independent MD trajectories (C, P, and L in water). However, a common practice is to compute the approximate binding affinity by simulating only “C” and generate “P” and “L” ensemble by simply deleting appropriate atoms. Thus, binding free energy could be written as

A popular choice of above mentioned single trajectory averaging is MM/PBSA. The binding free energy (ΔGbinding) is composed of three parts, molecular mechanics energy (ΔGMM), solvation free energy (ΔGsol), and entropic contribution (TΔS).

ΔGMM includes the electrostatic (ΔGELE) and Van der Waals (ΔGVDW) interaction energy terms. TΔS term is the entropic contribution to the binding free energy at temperature T. Because our objective was to compute binding free energy differences (ΔΔGbinding) of F-dUMP versus dUMP binding to protein (MtbThyX or hThyA), the entropic contribution of the protein part will certainly cancel. The relative entropic contribution of ligands (F-dUMP and dUMP) is expected to be very small (as the ligands are also very similar differing by only a single atom, H or F); thus, ΔGbinding was calculated by excluding the entropy part. ΔGsol was obtained by adding the electrostatic (ΔGPB) and nonpolar (ΔGSA) solvation binding free energies.

The electrostatic part of solvation free energy is obtained by numerically solving Poisson–Boltzmann equation[42−44]where T, ε, q, ρf, Φ(r), I, κ2, and e are referring to absolute temperature, dielectric constant, atomic number, fixed charge density, dimensionless electrostatic potential, ionic strength of the solution, Debye–Hückle parameter, and solute dielectric constant, respectively. The fixed dielectric constant was used for calculation. A grid spacing of 0.4, an ionic strength of 0.15 M, a temperature of 300 K, and a solute dielectric constant of 2 for MtbThyX and 1 for hThyA were used for calculating ΔGPB. Nonpolar solvation free energy was evaluated by calculating solvent accessible surface area using the formula[45]where γ and β are solvation parameters with a constant value of 0.00542 kcal/mol·Å2 and 0.92 kcal/mol, respectively.[45] CHARMM and GROMACS software was used to calculate the binding energy for 100 MD snapshots (evenly spaced) extracted from last 25 ns of all the trajectories generated from NAMD and GROMACS, respectively.[37,39,46−49] Default parameters of g_mmpbsa module (except solute dielectric of 1 was used for hThyA) implemented in GROMACS software was adopted to compute MM-PBSA free energy calculations.[43,44]

Electronic Structure Calculation

Ab initio quantum calculations were performed for keto, enol, and enolate forms of ligands (dUMP and F-dUMP) in the active site of MtbThyX. Representative MD structures of the ligand-binding site were selected at first. Then, the ligand, R199, R107, and CW1 were extracted. The sugar–phosphate moiety of the ligand and backbone of the amino acids (R199, R107) was replaced by a methyl group to reduce the computational cost. The resulting simplified ligand-binding site model was then subjected to gas-phase optimization. Gas-phase optimized coordinates were further optimized in the aqueous dielectric and subjected to normal mode analysis. M06-2X functional in combination with 6–311++G** basis set with an implicit SMD solvent model was used for performing ab initio calculations.[50,51] Calculations were performed with the Gaussian16 program.[52] The overall charge and the number of atoms were the same for both keto and enolate ligand-bound binding pocket model. The keto form of the ligand was modeled with deprotonated/neutral R199, whereas the deprotonated ligand was modeled with positively charged R199. Free energies were calculated by including the normal mode vibrational frequencies. Interatomic distances in the ab initio structures are more or less identical to the MD structures. The focus was to compute a quantity that approximates the relative ligand preference (keto vs enolate) at the active site. We computed the free energy difference of keto and enolate ligands in the simplified binding pocket model. Note, the entropy term was computed by including only the vibrational frequencies. Thus, the rotation and translation entropy were ignored.

Results

Comparison of MD and X-ray Structures of Mtb-ThyX

The X-ray structure of F-dUMP bound MtbThyX (Figure a) confirms the precise positioning of R199 (with a low average B factor of ∼14) close to the Watson–Crick edge of F-dUMP (distance < 2.8 Å).[16] Note, the relative position and/orientation of R199 and the Watson–Crick edge of substrate analogue (F-dUMP and Br-dUMP) are conserved and confirmed by experimentally determined structures of ThyX.[16,20] It has been hypothesized that hydrogen bonding between R199 and dUMP is required for the correct positioning of the substrate (dUMP) relative to cofactor (FAD) and facile methylation.[19] The proximity of R199 and F-dUMP Watson–Crick edge prompted us to examine three different models of ligands [keto (k), enol (e) and deprotonated/enolate (d)] in the MtbThyX binding pocket (abbreviated them as F-dUMP(k), F-dUMP(e), and F-dUMP(d), the same terminology is used for dUMP). To satisfy the hydrogen bonding requirement, models were considered either with charged R199 or neutral R199 (deprotonated) depending on the nature of the ligand viz., F-dUMP(d), F-dUMP(k) (Figure a,b).

Figure 2

MD (colored, after 50 ns post equilibrated trajectory) and X-ray (light grey) structures of the binding pocket: (a) enolate ligand (F-dUMP(d)) bound (b) keto (F-dUMP(k)) with deprotonated R199. (c) Average rmsd of protein heavy atoms during 50 ns of the trajectory for different models of MtbThyX (indicated by different colours). Trajectory averaged rmsd with standard deviation is given in Å. (d) Disruption of lid architecture in the absence of CW2. MD structures in the presence of CW2 (intact lid or gate-closed conformation) and in absence of CW2 (lid disruption or gate-open conformation). Absence of CW2 results in gate opening by the alteration of the position and orientation of the key residues (R107, Y108 and K165 shown by arrows), leading to water exposure (red mesh) of the ligand-binding pocket.

Structural similarly between MD and template X-ray structures of the ligand (dUMP, F-dUMP)-bound MtbThyX active site was analyzed by calculating the root mean square deviation (rmsd) of protein-heavy atoms relative to its template X-ray structure (PDB 3GWC).[16] MD trajectory average rmsd was found to be small (less than 1 Å) for different models, suggesting that MD structures are more or less similar to its template X-ray structure. rmsd versus time plots (Figure c) confirmed structural convergence during 50 ns MD trajectory. Enol tautomer of the ligand results in large rmsd (∼0.97 ± 0.07 Å) relative to its keto (∼0.82 ± 0.02 Å) and enolate (∼0.82 ± 0.02 Å) analogue (Figure c). MD structures of the MtbThyX binding pocket with bound F-dUMP(d) and F-dUMP(k) ligands are more or less identical to its template X-ray structure (Figure a,b). The ligand–protein hydrogen bonds observed in the X-ray structure are intact throughout the MD trajectory (Figure a,b, Tables S1, and S2). The enol tautomer, on the other hand, disrupts the local environment by losing two H-bonding interactions, F-dUMP(e): R107 (Figure S2, Table S1) and CW1: E92 (Figure S2). Binding pocket distortion rejects the possibility of enol tautomer of the ligand in the MtbThyX binding pocket. The crystal water (CW1) present in the binding pocket is stable in the MD trajectory and seems to act as a bridge between ligand (F-dUMP) and protein (E92). CW1 forms stable hydrogen bonds with O4 of dUMP and side chain of E92 (Figure a,b). The proximity of E92 side-chain and the phosphate moiety of F-dUMP observed in the X-ray structure (Figure a, Table S1) suggest two possibilities, either (1) E92 side-chain is protonated (Neutral E92) or (2) phosphate of F-dUMP is protonated. MD simulations suggest that the structural integrity of the binding site observed in the X-ray structure was reproduced more appropriately in the case of the protonated E92 model than protonated phosphate moiety (Figure S3), even though, the pKa of E92 side chain (pKa = 4.3) is less than the pKa of the phosphate of UMP (pKa = 7.2) at neutral pH = 7.[53,54] In response to protonation, the phosphate group of the ligand alters the orientations, resulting in the loss of salt-bridge interactions with R172 (Figure S3). The phosphate of F-dUMP forms direct interactions (Figure a) with the side-chains of two arginines and a serine (R87, R172, S105) and the backbone of Q106 (not shown). It can be argued that salt bridge interaction between the phosphate of F-dUMP and two arginines (R87 and R172) might be responsible for lowering the pKa, resulting in preferential deprotonation of phosphate of F-dUMP in the binding pocket.

Y108, K165, the backbone of H203, and a water molecule (CW2) seem to form a lid above the Hoogsteen edge of F-dUMP, possibly restricting water accessibility (Figure a). Simulation suggests that the presence of the water molecule (CW2) is crucial for the structural integrity of the binding pocket (Figure d). Absence of this crystal water results in severe disruption of the binding pocket, which include altered protein residues (see residues K165, Y108, and R107) and solvent exposer of the dUMP binding pocket (Figure d). MD trajectories indicate that the absence of CW2 disrupts the lid and allows water accessibility to the Hoogsteen edge of the F-dUMP binding pocket. Experimental studies showed that K165A mutation results in enzyme inactivity; thus, K165 is crucial for MtbThyX catalysis.[20] K165A mutation most likely disrupts the lid and exposes the binding pocket to water. Ligand–protein interactions observed in the X-ray structure were intact in the MD trajectory (Figure S4a,b, Table S1) for both dUMP(k)/F-dUMP(k) and dUMP(d)/F-dUMP(d) ligands. On the other hand, the binding pocket was observed to be severely disintegrated in dUMP(e) (Figure S2b, Table S1) because of the loss of a key hydrogen bond between dUMP(e) and protein side-chain [R107 with an average heavy atom bond distance of 5.54 ± 0.54 Å (Table S1)].

Comparison of MD and X-ray Structures of hThyA

Using MTX and dUMP-bound human ThyA structure (PDB 5X66) as a template (Figure b), we performed simulations with bound dUMP(k) and F-dUMP(k) ligands.[22] Ligand-binding pocket of ThyA (Figure b) is distinctly different from MtbThyX (Figure a). Polar residues (side-chain of N226, backbone −NH of D218) form H-bonds with the Watson–Crick edge of the ligand in hThyA. Thus, enol and enolate forms of the ligand is unlikely in hThyA (clearly disrupt ligand-N226 interaction). Therefore, we have considered the keto form of the ligand in hThyA. To understand the role of MTX toward the structural integrity in hThyA, we performed simulations in the presence (ternary complex) and absence (binary complex) of MTX (Figure a,b). rmsd calculations suggest that F-dUMP bound ThyA has higher average rmsd (1.53 ± 0.06 Å) than its dUMP bound analogue (rmsd 1.19 ± 0.06 Å) (Figure S1). MD simulations suggest that overall flexibility of hThyA protein increases in the absence of MTX (binary complex), compared to the MTX bound ternary complex (Figure c). A drastic increase in the RMSF of C-terminal end (residue 309–313) was observed in case of binary complex (RMSF of 5.12 ± 1.87 Å) relative to the ternary complex of hThyA (RMSF of 1.31 ± 0.37 Å) (Figure c). It can be argued that the C-terminal end is crucial for folate/antifolate binding in hThyA, as reported earlier.[22] MD trajectory of the F-dUMP bound ternary complex suggests that the side chain of Y135 is involved in h-bond interaction with the fluorine (63.2% of the 50 ns trajectory with an average heavy atom bond-distance of 3.42 ± 0.36 Å, Figure d). Y135 seems to satisfy the h-bonding requirement of the Hoogsteen edge of F-dUMP by forming a H-bond with the F atom of F-dUMP during MD trajectory, as suggested previously.[25] The X-ray structure of ThyA in Escherichia coli species has shown Y135 and F-dUMP interaction, which is in line with our MD results of hThyA.[55] Hoogsteen edge of the ligand in the hThyA binding pocket is wet (accessible to bulk water).

Figure 3

MD (colored, after 50 ns post equilibrated trajectory) and X-ray (light grey) structures of hThyA: (a) ternary complex, that is, with antifolate (MTX, line), (b) binary complex, that is, without antifolate (MTX). Position of C-terminal (residue 309–313, in cartoon) was altered relative to its X-ray structure in the case of a binary complex, shown in black arrow. (c) Residue-wise RMSF of protein-heavy atoms in the binary and ternary complexes. C-terminal region is shown in the rectangular box. (d) F-dUMP bound complex shows direct interaction between Y135 with fluorine of the ligand (F-dUMP). The same interaction is absent in the dUMP-bound X-ray structure of the complex, and Y135 was found to be distant from the Hoogsteen edge of the substrate (grey).

Keto Versus Enolate Preference in MtbThyX—Quantum Chemical Calculations

Classical MD simulations suggest that both the keto (dUMP(k)/F-dUMP(k)) and enolate (dUMP(d)/F-dUMP(d)) ligands in the MtbThyX binding pocket can preserve the interaction network observed in the X-ray structure (PDB 3GWC, Figure a). On the other hand, consideration of enol tautomer of the ligand results in binding pocket disruption during MD production runs (Figure S2). Comparison of structural integrity between X-ray and MD structures suggests that both keto and enolate ligands are the reasonable models in MtbThyX binding pocket. We first addressed the question of keto versus enolate ligand preference in the MtbThyX binding pocket using quantum chemical calculations. Keto/enolate ligand, R199, R107 and CW1 were selected from a representative MD snapshot (Figure ) and subjected to ab initio calculation by methyl capping (see Methods). Ab initio calculations confirmed that the simplified models of the ligand-binding pocket of MtbThyX were true minima on the potential energy hyper-surface, supported by positive normal mode frequencies (Table S5).

Figure 4

Ligand-binding pocket after ab initio optimization (a) keto ligand (R-dUMP(k), R = H, F) with neutral R199. (b) Enolate ligand (R-dUMP(d), R = H, F) with charged R199. CW1 is crystal water. Ligand, R107 and R199 were subjected to methyl capping. N3–H3 bond distance (in Å) of the ligand in the binding pocket (complex) and free form is given within the box. Ab initio Free energy difference between model “b” and “a” is given as ΔG (in kcal/mol).

Free energy calculations in the continuum water dielectric suggest that enolate F-U(d)/U(d) is energetically more favored over keto F-U(k)/U(k) by 9.5/8.6 kcal/mol (Figure ). The stabilization of the F-U(d)/U(d) over its keto analogue, however, has to be corrected with the free energy cost required for deprotonating the uridine or arginine (R199) in solution at pH 7. With a pKa of 9.2 (for N1 of uridine)[56] and 12.1 (for arginine side-chain),[57] this correction thus becomes +3 kcal/mol and +7 kcal/mol for model “b” and model “a”, respectively (Figure ), resulting in a net effect on the relative free energy difference of 13.5/12.6 kcal/mol in favor of enolate ligand F-U(d)/U(d) compared to the neutral keto F-U(k)/U(k) in the MtbThyX binding pocket. Noticeable N3–H3 bond-length elongation (by 0.05–0.06 Å) in the binding pocket was observed for the keto ligand relative to its unbound free form (Figure ). Direct interaction between R199 and R107 with the ligand results in N3–H3 bond-length elongation. Note, N3–H3 bond-length elongation is reduced (by 0.02 Å) if optimization is done by placing R107 away from the ligand (Figure S5). Contrary to the earlier proposal[17] that R199 is solely responsible for the ionizing dUMP in the ThyX pocket, our results suggest that R107 is likely to play a key role in stabilizing the enolate form of the ligand. Ab initio result suggests that the enolate ligand is preferred in the MtbThyX binding pocket, relative to the keto form of the ligand.

Ligand Binding Affinity—MM/PBSA

Absolute binding free energy was estimated for ligand (dUMP and F-dUMP) binding to MtbThyX and hThyA using the MM/PBSA method (Table ). The computed binding free energies reveal several remarkable features. First, the enolate ligand is favored over keto form in the MtbThyX binding pocket, corroborating with quantum chemical predictions. Binding affinity is predicted to be the lowest for enol tautomer (relative to keto/enolate form of the ligand) to the MtbThyX binding pocket (Table S6). Second, F-dUMP is always favored over dUMP (Figure ). Third, the discriminatory power of hThyA in favor of F-dUMP binding relative to dUMP is larger (ΔΔG ∼ 6–7 kcal/mol) than MtbThyX (ΔΔG ∼ 2–3 kcal/mol). Thus, we may conclude MtbThyX is weakly selective and hThyA is highly selective for favorable F-dUMP binding relative to substrate dUMP. Fourth, the discrimination strength of MtbThyX in favoring F-dUMP binding over dUMP is more or less independent of the protonation state (keto or enolate) of the ligand (Table , Figure ). We performed a statistical significance test (z-test) for the estimated ΔΔG’s. Estimated Z-statistic (Z0 > 2.5 and p-values < 0.01) indicate that the sign of (ΔΔG) is correct and statistically significant (greater than 99% confidence level) (Figure S6). Molecular mechanics term (ΔGMM) primarily drives ligand-binding to the protein, whereas polar solvation term (ΔGPB) disfavors ligand-binding (Table ). ΔGMM includes electrostatic (ΔGELE) and Van der Waals (ΔGVDW) part with the former as a major contributor for favorable ligand-binding (Table S6). Non-polar solvation (ΔGSA) was favorable for ligand-binding, but it is insignificant relative to the ΔGMM.

Figure 5

Computed MM/PBSA binding free energy differences (ΔΔGbind, error bars in 1 s.e.m) of F-dUMP and dUMP binding to MtbThyX and hThyA. Positive sign of ΔΔGbind indicates favorable F-dUMP binding relative to dUMP. The magnitude of ΔΔGbind indicates the strength of selectivity. Results obtained by using two different software [CHARMM (grey) and GROMACS (white)] with different MD models (see Methods) were compared with experiment[12] (green).

Discussion

Understanding the mechanism in terms of atomic structure is the key for designing therapeutic molecules for combating disease. In this paper, we focused on understanding the dynamics and energetics of ligand-binding to MtbThyX and compared the same with its human analogue hThyA. The proximity of R199 and Watson–Crick edge of pyridine base of the ligand in the X-ray structure of the MtbThyX binding pocket (Figure a) indicates different possible models for the pair, (1) F-dUMP(k): neutral R199, (2) F-dUMP(d): positively charged R199, and (3) F-dUMP(e): neutral R199. We propose, model (2), that is, enolate dUMP(d)/F-dUMP(d) as the most appropriate ligand for the MtbThyX ligand-binding pocket (Figure a). Consideration of model (3), that is, enol dUMP(e)/F-dUMP(e) as a possible ligand in the MtbThyX binding pocket disrupts the structural integrity (Figure S2), thus being unlikely to be a correct model. MD simulations of MtbThyX with either model (1) or model (2) result in very similar structures. Quantum chemical calculations suggest that the keto form of the ligand is energetically less favorable in the binding site by ∼9 kcal/mol relative to its deprotonated form. The relative preference is boosted further (by 3 kcal/mol) if deprotonation free energy of R199 and dUMP is considered at neutral pH. NMR studies based on C14 chemical shift on ThyX isolated from Thermotoga maritima proposed the presence of deprotonated dUMP in the ligand-binding site, in line with our results.[17] R199- and R107-induced N3–H3 bond length elongation (Figure ) suggests that arginines (R199, R107) are most likely play a key role in facilitating deprotonation of keto dUMP in the MtbThyX binding pocket. Note, displacement of R107 away from the ligand (keto-form) results in shortening of N3–H3 bond-length (Figure S5). The X-ray structure clearly shows that E92 side-chain is located very close to the phosphate of the substrate dUMP, which suggests the protonated form of either E92 (side-chain) or dUMP (phosphate) in the substrate-bound MtbThyX. It is energetically more favorable to protonate the phosphate of UMP (pKa = 7.2) at neutral pH in solution than to protonate the side-chain of E92 (pKa = 4.3).[53,54,58] However, MD simulations support protonated E92 as a possible model. Consideration of phosphate protonation in the dUMP ligand disrupts the binding pocket during MD simulation (Figure S3). Direct interaction between phosphate of dUMP and two arginine’ (R87, R172) favored deprotonated dUMP phosphate in the MtbThyX (Figure a,b). The lid composed of Y108, K165, H203, and CW2 (Figure a) acts as a gate and seems to provide hydrophobic cover to the ligand-binding pocket. Absence of CW2 in the MD setup clearly shows binding pocket disruption and water accessibility through the gate (Figure d). Experiments also suggested that mutation of the gate residue (K165A) results in loss of enzyme activity, in lines with our findings.[20] On the contrary, the Hoogsteen edge of the ligand in the hThyA binding pocket is exposed to the bulk water. Human ThyA prefers the keto form of the substrate in its binding pocket. Simulation of ligand (dUMP or F-dUMP) bound ThyA in presence and absence of MTX indicate that MTX plays a key role in stabilizing the C-terminal end of hThyA (Figure a,b). The overall flexibility of the ThyA also increases substantially in the absence of MTX (Figure c). Previously, it has been proposed that methyl-tetrahydrofolate (mTHF), which occupies the same position as MTX in the ThyA catalytic pocket, increases the stability of ThyA.[22] The same argument holds in case of antifolate MTX. Folate (mTHF)/antifolate (MTX) seems to induce conformation transition (highly flexible (binary) → less flexible (ternary) conformation) by stabilizing the C-terminal part of hThyA.

Both hThyA and MtbThyX prefer F-dUMP binding over dUMP (Figure ). Note, the strength of discrimination of hThyA (ΔΔG ∼ 6–7 kcal/mol) favoring F-dUMP over dUMP is stronger than MtbThyX (ΔΔG ∼ 2–3 kcal/mol). These results complement the experimental studies, which proposed that F-dUMP is a more potent inhibitor of human ThyA with a Ki value of 1.7 nM compared to MtbThyX (Ki = 100 nM).[12] Similar Michaelis constant (Km) associated with dUMP binding to MtbThyX (Km = 3 μM) and hThyA (Km = 2.5 μM) suggests that binding affinity of dUMP is more or less the same for MtbThyX and hThyA. Experimental relative binding free energy (F-dUMP vs dUMP, favouring former) is estimated to be ΔΔGexp ∼ 2 kcal/mol (for MtbThyX) and 4.5 kcal/mol (for hThyA) at 37 °C (assuming Km ≈ Ki).[12] The magnitude and sign of relative binding free energies calculated from MM/PBSA method (ΔΔGbind) are in excellent agreement with the experiments (Figure ). Human ThyA discriminates between F-dUMP(k) and dUMP(k), strongly favoring the former (by ∼6–7 kcal/mol) by establishing new protein–ligand interaction involving Y135 side-chain and fluorine of the F-dUMP ligand (Figure d). Note, residue Y135 which is located away from the Hoogsteen edge of dUMP reorients and forms a new H-bond with the fluorine in case the of F-dUMP (Figure d). The simulations reveal that Y135 plays a key role in amplifying the strength of selectivity. On the other hand, in the MtbThyX binding pocket, R107 is located close to the Hoogsteen edge, irrespective of the nature of the ligand, dUMP or F-dUMP. R107 is solely responsible to satisfy the hydrogen-bonding requirement of the Hoogsteen edge of the ligand. Thus, R107 forms either single H-bonding interaction with O4 (in case of dUMP) or bifurcated H-bonds with O4 and F (in case of F-dUMP), leading to weak discrimination between dUMP(d) and F-dUMP(d) (Figure a,b). A free energy decomposition analysis indicates that R107 (of MtbThyX) and Y135 (of hThyA) residues primarily contribute to the estimated ligand selectivity (ΔΔG) by its electrostatic term (ΔΔGELE ∼4 and 1.5 kcal/mol for R107 and Y135, respectively). A detailed free energy decomposition analysis and its structural basis for ligand recognition (dUMP vs F-dUMP) are outside the scope of this paper and will be reported elsewhere.

Conclusions

In the present study, a structure-based MD study was performed for ligand/inhibitor (dUMP/F-dUMP) bound MtbThyX and hThyA protein complexes. The architecture of the ligand-binding pocket (viz., the proximity of R199 and Watson–Crick edge of inhibitor F-dUMP in the X-ray structure MtbThyX, Figure a) suggests different possible forms of the ligand (keto, enol, and enolate), thus raising ambiguity in the atomic model of the binding pocket. Classical MD suggest the consideration that both enolate and keto forms of the substrate preserve the structural integrity of the MtbThyX binding pocket (Figure a,b), whereas enol form causes severe distortion of the interaction network in the binding pocket. Gate formed by Y108, K165, and H203 is in a closed conformation in the presence of CW2 (Figure d). Absence of CW2 opens the gate and allows water entry to the ligand-binding pocket, thus being crucial for structural stability of the binding pocket (Figure d). The gate above the Hoogsteen edge of the ligand in MtbThyX provides dry-desolvated pocket, contrary to hThyA, which is wet and water-exposed. Residue E92 of the MtbThyX pocket is proposed to be neutral. The protonation state of the ligand and E92 is shown to influence the structural stability of the MtbThyX protein. Electronic structure calculations suggest that the enolate form of the ligand is energetically preferred relative to its keto analogue in MtbThyX (Figure ). Deprotonation of the ligand seems to be facilitated by two arginines (R199 and R107), indicated by bond-length elongation in the keto form of the ligand (Figure ). The binding pocket of hThyA includes the keto form of the ligand. The presence of MTX lowers the overall flexibility of hThyA protein, with the maximum reduction of RMSF in the C-terminal end (Figure c). Both MtbThyX and hThyA prefer inhibitor (F-dUMP) binding over its substrate (dUMP). Substrate/inhibitor binding to protein is stabilized primarily by molecular mechanics terms (electrostatic and van der Waals interactions) (Table ). Calculations of binding free energy differences (using MM/PBSA) indicate that hThyA is more selective (ΔΔG = 6–7 kcal/mol) between the substrate (dUMP) and inhibitor (F-dUMP), favoring latter, relative to MtbThyX (ΔΔG = 3–4 kcal/mol) (Figure ). The protonation state of the pyrimidine base of the ligand (keto or enolate) is found to play no significant role in tuning the discriminatory power of MtbThyX (favouring F-dUMP over dUMP binding) (Table , Figure ). Y135 plays a crucial role in boosting the discriminatory power of hThyA by forming new protein–ligand interaction, involving the fluorine atom of F-dUMP and side-chain of Y135 (Figure d).

We propose (1) protonation state of the substrate is very different between MtbThyX and hThyA. MtbThyX prefers the deprotonated/enolate dUMP ligand, whereas hThyA prefers keto dUMP as its substrate. (2) MtbThyX ligand-binding pocket is relatively dry to its human analogue (hThyA). (3) hThyA is more selective than MtbThyX in favor of inhibitor (F-dUMP) binding relative to its substrate (dUMP). To the best of our knowledge, this is the first MD study of MtbThyX, which has attempted to establish the link between 3D structures and free energy. The study could be considered as a stepping stone for new antitubercular drugs. We propose that a molecule with an appropriate substituent that can replace the crystal water (CW1) and mimic the deprotonated form of dUMP may be a new inhibitor for MtbThyX with enhanced selectivity. Moreover, a small molecule that can bind to an allosteric site and disrupt the molecular gate in the MtbThyX may also inhibit dUMP → F-dUMP catalysis in M. tuberculosis.

Despite the limitations (viz., sampling, convergence, force fields accuracy), we believe that the present study is valuable for (1) characterizing the key protonation states of the binding pocket residues, (2) establishing a bridge between atomic structures and free energies, and (3) encouraging further experimental investigations.