Molecular mechanics parameters for the PLP cofactor were derived in this work, in forms that were covalently bound to lysine (Figure 1). Four protonation states were treated: the phosphate group was singly- or doubly-charged, and the NZ-O₃ atom pair had a proton either on NZ (zwitterionic, keto form “Z”) or O₃ (neutral, enol form “N”). The forms are referred to as Z1, Z2, N1, and N2.

Atomic partial charges were calculated using Gaussian09 (Frisch et al., 2009), with the Hartree-Fock method and the 6-31G* basis set (HF/6-31G* method) (Hehre et al., 1987), consistent with the Amber protein force field (Cornell et al., 1995; Maier et al., 2015). The PLP moiety was covalently bound to a fragment of the Lys side chain, truncated after the CD atom (with the CD methylene replaced by a methyl). Charges were adjusted, as usual, to reproduce the electrostatic potential at nearby points selected according to the Merz-Singh-Kollman scheme (Cornell et al., 1995). Atom types for the PLP moiety were assigned by analogy to standard amino acids or nucleotides in the Amber ff14SB force field. Van der Waals parameters were then available automatically. Geometric parameters were taken from the optimized ab initio geometries, and force constants were inferred by analogy to existing groups, with one exception. For the Schiff base dihedral angle χ = C4-C4′-NZ-CE, a more complex procedure was used, where the distribution sampled in MD simulations was adjusted to match the PDB distribution; see details in Results.

Simulation systems were prepared using Xplor (Brünger, 1992), starting from the ppMGL crystal structure (PDB code 1GC2), which is a homotetramer. The missing loops in 1GC2 were taken from 2O7C, structurally aligned, and minimized by applying harmonic restraints at the junction points, using Xplor. Each model involved a protein tetramer, spherically truncated to eliminate parts far away from the active site. Truncation spheres were centered on the PLP in one specific active site and had radii of 28 Å. Each truncated tetramer was solvated by an octahedral TIP3P water box (Jorgensen et al., 1983). Water molecules having an oxygen less than 2.5 Å from any nonhydrogen protein atom were deleted. A few sodium ions (around 5, depending on the system) were added, to neutralize the total protein charge. The water box had an edge length of 78 Å.

MD simulations were performed with the Amber ff14SB protein force field (Maier et al., 2015) and PLP parameters developed here. Simulations were done with NAMD (Phillips et al., 2005). They used periodic boundary conditions, with a truncated octahedral box and the Particle Mesh Ewald method for long-range electrostatics (Darden et al., 1993). A 12 Å cutoff was applied to van der Waals and real-space electrostatic interactions. Temperature was maintained at 295 K using Langevin dynamics for nonhydrogen atoms with a damping coefficient of 5 ps⁻¹. A pressure of 1 bar was maintained by the Langevin piston Nose-Hoover method, with an oscillation period of 200 fs and a damping time of 100 fs (Martyna et al., 1994; Feller et al., 1995).

All simulations followed the same setup procedure. Systems were first minimized for 1,000 conjugate steps, followed by 100 ps of equilibration in the NVT ensemble, then four segments of 500 ps each in the NPT ensemble, with decreasing harmonic restraints in each equilibration segment. For the MD production, only the residues within 5 Å of the edge of the protein truncation sphere were restrained. The harmonic force constants gradually decreased going from the edge inwards, with four concentric shells having values of 2, 1, 0.5 and 0.25 kcal/mol/Ų.

We analyzed 27 PDB structures that all contained the PLP cofactor, covalently bound to lysine (Lys211 in ppMGL): 11 of ppMGL and 16 of orthologs. Removing symmetry images within each MGL tetramer left 63 independent monomers. The list of structures, including information on ligands and crystallization conditions, is in Supplementary Material, Table S1. Selected structural variables were extracted using Xplor, related to PLP and nearby side chains.

The fragments used for charge calculations are shown in Figures 1, 2. PLP in its lysine-bound form was represented by a fragment where the Lys side chain was truncated (after CD) and methylated. To re-establish the correct total charge in the context of a full lysine side chain, the ab initio charges obtained for atoms close to the junction were adjusted slightly by hand. Atom types were inferred by analogy to groups in the Amber force field. Geometric parameters were taken from the optimized ab initio geometry, with one exception, described next. Final parameters are reported in Supplementary Material, in a format compatible with Xplor (Brünger, 1992) and NAMD (Phillips et al., 2005).

To optimize parameters for the Schiff base dihedral angle χ = C4-C4′-NZ-CE, torsion energy scans were initially performed for both the neutral and zwitterionic forms of the lysine-bound PLP. A model compound was used, shown in Figure 3, having no phosphate group. We restrained the χ dihedral angle to a series of values from −180° to 180° and minimized each structure with molecular mechanics, using a restraint force constant of 1,000 kcal/mol/Ų and the Xplor program. The minimized structures served as the starting points for quantum mechanical geometry optimization, with a dihedral constraint and the HF/6-31G* method. Care was taken to maintain the orientation of the groups surrounding the Schiff base linkage, to ensure a smooth variation with χ and good consistency between the conformations optimized with molecular mechanics and quantum mechanics. For χ angles from −180° to −120°, we re-optimized the structures and calculated the energies with the MP2 method and the 6–311 (1d, 1p) basis set, chosen for consistency with similar torsion energy scans for dihedral parametrization (MacKerell et al., 2015). The molecular mechanics parameters were fitted by comparing the energies from Xplor and Gaussian. Fitted curves are shown in Figure 3.

The fitted torsion parameters were then employed in MD simulations of the enzyme with PLP covalently bound. Despite the accurate fit, the values of χ sampled in the simulations did not match the values observed in MGL crystal structures, as shown in Figure 3. Therefore, the parameters were adjusted manually and new simulations were performed. After iterating this procedure several times, the MD simulations reproduced approximately the χ distribution inferred from the PDB. The initial and final torsion parameters are reported in Supplementary Material. The parameters for the zwitterionic tautomer (Z1/Z2) reproduced the PDB distribution closely, while the fit was somewhat less good for the N1/N2 forms, despite extensive fitting efforts. The final parameters were then used for all the MD simulations described below.

Before describing MD simulations of MGL-PLP complexes, we present some aspects of the MGL architecture and the PLP environment. MGL belongs to the cystathionine synthase-like SCOP family, within the PLP-dependent transferase superfamily (Andreeva et al., 2008). In the PDB, 47 MGL crystal structures are found, from seven species. The structures are homotetramers. Each monomer contains an N-terminal domain, a PLP-binding domain and a C-terminal domain. The active sites (four per tetramer) are at the monomer interface within a catalytic dimer.

27 of the PDB structures have PLP bound in the active site through a covalent bond to Lys211 (ppMGL numbering). Protonation states of several active site groups, including PLP, are unknown. The C4′-NZ chemical bond linking PLP to Lys211 can be seen as part of an imine Schiff base, and is expected to have a partial double bond character (Jencks, 1986). Nevertheless, in the PDB structures, the C4-C4′-NZ-CE atoms deviate from a planar geometry, as described in the previous section. Eight of the 27 PDB structures (Table S1) have a crystallization anion directly bound to the PLP NZ atom, which is presumably protonated in these cases (keto form). The anion site is close to the position where the carboxylate of the Met substrate will ultimately bind. All 27 structures were crystallized with either sulfate (0.1–1.8 M) or chloride (50 or 250 mM) anions. The NZ atom has two distinct orientations in the 63 monomers, revealed by the distances between NZ and O3 and OG (Ser208) (ppMGL numbering; Figure 4B). The first orientation, “A”, is towards the anion site. In this orientation, NZ and O3 form a hydrogen bond. “A” occurs in 25 structures, solved with either sulfate or 50 mM chloride crystallization buffers. Six have a visible anion, always sulfate. The second orientation, “B”, has NZ turned inwards, and more buried. It occurs in four structures, including one where A and B are seen in different monomers. All four have chloride crystallization conditions; one (5DX5) has a visible Cl⁻. In two of them (1GC0, 1GC2, both ppMGLs), there is negative density in the Fo—Fc electron density map at the NZ position, which suggests A may also be partially occupied, even though it is absent from the PDB coordinate file. Orientation B gives rise to an NZ-OG (Ser208) hydrogen bond, allowed by both NZ protonation states (since OG can be either donor or acceptor). Though NZ is turned towards the phosphate in the four B structures, it is too far to form a hydrogen bond (Figure 4B). Overall, orientation A is predominant, especially with sulfate or low chloride concentrations.

The PLP aromatic ring stacks upon Y114 (Figure 4A). Its nitrogen atom N1 interacts closely with D186, which indicates (Jencks, 1986) that N1 is protonated. The PLP phosphate forms hydrogen bonds with the S208 and T210 side chains and the G89 and M90 backbone NH groups. In most of the structures, the phosphate interacts with the side chains of Y59* and R61*; the asterisk indicates they are contributed by the other monomer of the catalytic dimer. Out of 63 subunits, 10 have the 52–62 loop disordered (ppMGL numbers); in the other 53, the interactions are present in all but 7 (Y59*) and 5 (R61*), for an occupancy of around 75%. R61* also hydrogen bonds to Y114 in 45 of 63 of the PDB subunits (71%). All these residues are conserved across orthologs. All the interactions are compatible with either a singly-protonated or a fully-deprotonated phosphate. We also expect that either atom O3 or NZ should carry a proton (Figure 1).

With force field parameters in hand, we proceeded to simulate ppMGL with bound PLP in four distinct protonation states: −1 or −2 phosphate and zwitterionic (keto) or neutral (enol) form of the Schiff base linkage. We refer to these forms as Z1, Z2, N1, N2 (Figure 1). For each form, we ran two MD simulations of 600 ns each, for a total of 4.8 microseconds of simulation. Rms deviations from the starting structures were stable on these timescales, indicating that all the models were reasonably stable, with deviations of 0.5–1 Å for backbone atoms and 1.5–2 Å for side chains. We first describe the interactions of the PLP phosphate group, and compare results for the singly- and doubly-charged forms. Next, we consider other PLP interactions, and compare the results for the neutral and zwitterionic forms.

Phosphate deprotonation was first simulated for the model compound in solution (upper leg, Figure 2). The protonated compound was prepared through 20 ns of MD. The transformation was simulated in five steps, with five distinct λ values and 20 ns of MD for each. Each simulation is referred to as a window. The endpoint (or midpoint) of one window served as the starting point for the next. The last frame of the λ = 1 window was used to initiate the reverse transformation. Postprocessing was done with BAR. Agreement between the forward and backward runs was excellent, with ΔG _S values of −80.4 and −80.5 kcal/mol, indicating that sampling was adequate. Postprocessing with thermodynamic integration (TI) gave the same results, within 0.1 kcal/mol (Table 1).

In proteins, ionic transformations like deprotonation are challenging and require careful sampling. The different FEP runs for the protein complex are schematized in Supplementary Material, Supplementary Figure S1. The complex was first simulated in its N1 state (Figure 5). The MD frame after 50 ns was used to initiate a “forward” FEP run (run A). Each λ value was simulated for 40 ns. Each subsequent window was started from the midpoint of the previous window (MD frame after 20 ns). For the final, λ = 1 window, the last frame was used to initiate a backward run, done in the same way (run A′). Postprocessing with BAR, the two estimated free energy changes for deprotonation ΔG _P were −79.1 and −87.0 kcal/mol. The large difference indicates insufficient conformational sampling, despite 400 ns of total MD. Per-window contributions from BAR are shown in Table 1, along with free energy derivatives from TI. The A/A’ difference is seen to arise mainly from the λ = 0 and 0.25 windows.

To identify the structures involved, we collected the same structural variables as above. Results are in Figure 6. The main A/A′ difference involves the interaction between the PLP phosphate and Tyr59* from the second subunit. Above (Figure 5), we described a total of 1,200 ns of MD each for the singly- (N1, Z1) and doubly- (N2, Z2) protonated states. The occupancy of the phosphate–Tyr59* interaction was 100% for N2+Z2, presumably due to the −2 phosphate charge. It was around 8% for N1+Z1. In run A, the interaction was present about 1/4 of the time near the N1 endpoint (λ = 0, 0.25), then maintained throughout. In the reverse run A′, the interaction was maintained during the entire run. The stronger interaction in the λ = 0, 0.25 windows of run A’ led to a more negative ΔG _p value, as expected.

Two additional forward/backward pairs of runs were done in the same way (runs B/B′ and C/C’; Supplementary Figure S1), using state N1 as the starting point, but choosing starting MD frames where the Tyr59*–phosphate was broken (Figure 5). In these runs, each subsequent λ window was initiated from the endpoint, not the midpoint of the previous window, thus benefiting from a more thorough equilibration. The structural variables from all these runs are included in Figure 6. In all four runs, the Tyr59*–phosphate interaction never formed, even at the λ = 1 endpoints (N2 state). BAR postprocessing gave deprotonation free energies ΔG _P of -74.5/-74.6 kcal/mol for B/B′ and -70.1/-74.3 kcal/mol for C/C’. TI values were very similar, with a mean unsigned TI/BAR difference of 0.3 kcal/mol. Thus, B/B′ and C/C′ were in reasonable agreement with each other and in fair agreement with run A. The deprotonation energies for the first half of the three runs were a bit less favorable than for run A, which is expected, since the Tyr59*–phosphate interaction was present about 1/4 of the time in run A. Run C’ gave somewhat less negative values for all windows.

Finally, another pair of runs was done, which started from the midpoint λ = 0.5 value then ran until λ = 1 and back (runs d, d’; lowercase names recall that these runs spanned only part of the λ range). Window lengths were doubled to 80 ns. The initial frame for run d at λ = 0.5 was chosen so that the Tyr59*–phosphate interaction was formed. The interaction was then maintained throughout the two runs. BAR and TI results are included in Table 1.

Overall, based on the 2,400 ns of equilibrium of MD and 1,680 ns of FEP simulations, it appears that for the λ ≤ 0.25 range, runs A/B/B’/C/C′ are representative of the most important conformations, with the Tyr59*–phosphate interaction formed 10% of the time (1/4 of run A′). For the range λ ≥ 0.75, runs A/A’/d/d’ are representative. For the midpoint window, we assumed we can average over B/B’/C/C’/d/d’. To implement this averaging, we focussed on the TI data, and averaged the free energy derivative values over runs appropriately (A/B/B’/C/C’ for λ = 0 and 0.25, and so on). The mean derivatives (Table 1) were then integrated, giving a total deprotonation free energy of −78.3 kcal/mol. This total does not include box size corrections.

Box size corrections were detailed in Methods. Briefly, PME imposes a shift of the electrostatic potential, such that the mean potential throughout the box is zero. For the protein solution, the mean potential far from the protein was obtained by averaging over 100 snaphots from a 10 ns MD run. Around 97,000 grid points out of 372000 were over 12 Å from the protein and considered to represent the bulk solvent region. Averaging over these grid points and MD snapshots, the mean potential was found to be Φ = −2.9 kcal/mol/e. The mean potential in the model compound solution was zero. Upshifting the protein box potential by − Φ brings the potential in the bulk-solvent region to the same level as in the model compound simulation. This shift added a contribution to the deprotonation free energy that was −1 ×Φ = 2.9 kcal/mol. A second free energy correction (Methods) was considered to be negligible.

Thus, we obtained a final value of ΔG _p = −81.2 kcal/mol, and ΔG _p − ΔG _s = −0.8 kcal/mol. The phosphate pK _a in the protein is given by Eq. 1. The model compound has a solution pK _a of 6.3 (Limbach et al., 2011), so that the pK _a in the protein is 5.7. Thus, the PLP phosphate group is predominantly dinegative in MGL. At a physiological pH of 7.5, the dinegative state is predicted to have a population of over 98%.

Table 2 reports the N2 → Z2 free energy change in the protein, using a TI method. The free energy derivatives from two separate 600 ns simulations of each state are reported. They were integrated with a trapezoidal method, yielding the free energy change Δ G N 2 → Z 2 prot (MM), where MM stands for “molecular mechanics”. In the gas phase, we computed the energy for a single N2 structure and a single Z2 structuree. The gas phase free energy change Δ G N 2 → Z 2 gas (MM) is simply the difference. The gas phase value from quantum mechanics Δ G N 2 → Z 2 gas (QM) = -0.4 kcal/mol was obtained by Lin and Gao (Lin and Gao, 2010) (B3LYP/6-311G (d,p)//B3LYP/6-311G (d,p) level of theory). Finally, the free difference in the protein was computed as Δ G N 2 → Z 2 prot = Δ G N 2 → Z 2 prot M M − Δ G N 2 → Z 2 gas M M + Δ G N 2 → Z 2 gas Q M (3) The final result was Δ G N 2 → Z 2 prot = −1.3 ± 1.0 kcal/mol. The uncertainty was estimated by taking results for the two halves of the protein MD segment for each endpoint. The FEP calculations thus indicate a moderate preference for the zwitterionic, keto state Z2. This is consistent with our MD/PDB comparison, above. The preference is, coincidentally, the same as the value obtained by Lin and Gao for another PLP enzyme (Lin and Gao, 2010).