Article Cite This: J. Phys. Chem. B 2018, 122, 5038−5048
pubs.acs.org/JPCB
Accurate Methyl Group Dynamics in Protein Simulations with AMBER Force Fields Falk Hoffmann,† Frans A. A. Mulder,‡ and Lars V. Schaf̈ er*,† †
Theoretical Chemistry, Ruhr-University Bochum, D-44780 Bochum, Germany Interdisciplinary Nanoscience Center (iNANO) and Department of Chemistry, University of Aarhus, DK-8000 Aarhus, Denmark
‡
Downloaded via WESTERN UNIV on October 26, 2018 at 05:10:54 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.
S Supporting Information *
ABSTRACT: An approach is presented to directly simulate the dynamics of methyl groups in protein side-chains, as accessible via NMR spin relaxation measurements, by all-atom MD simulations. The method, which does not rely on NMR information or any system-specific adjustable parameters, is based on calculating the time-correlation functions (TCFs) of the C−H bonds in methyl groups and explicitly takes the truncation of the TCFs due to overall tumbling of the molecule into account. Using ubiquitin as a model protein, we show (i) that an accurate description of the methyl dynamics requires reparametrization of the potential energy barriers of methyl group rotation in the AMBER ff99SB*-ILDN force field (and related parameter sets), which was done with CCSD(T) coupled cluster calculations of isolated dipeptides as reference, and (ii) that the TIP4P/2005 solvation model yields overall tumbling correlation times that are in close agreement with experimental data. The methyl axis squared order parameters S2axis and associated correlation times τf, obtained within the Lipari−Szabo formalism, are in good agreement with the values derived from NMR deuterium relaxation experiments. Importantly, the relaxation rates and spectral densities derived from MD and NMR agree as well, enabling a direct comparison without assumptions inherent to simplified motional models.
■
protein backbone dihedral angles.25,26 However, while backbone motions can now be reliably predicted by MD on a routine basis, accurate simulation of amino acid side-chain dynamics, as for example probed by deuterium relaxation of methyl groups, is more challenging.27−31 NMR relaxation rates of 13CH22H methyl groups in 13C-, fractionally 2H-labeled proteins report on reorientation motions of 13C−2H bond vectors on time scales faster than the overall tumbling of the molecule, i.e., from picoseconds up to a few nanoseconds for a typical small protein. Measurement of up to five different relaxation rates at a given magnetic field yields spectral densities at discrete frequencies J(0), J(ωD), and J(2ωD), where ωD is the Larmor frequency of deuterium.32 In addition to the dynamic parameters derived from the data with the LS formalism, such as S2 and τf, it is beneficial to also compare the experimentally directly accessible relaxation rates as well as spectral densities to MD simulations, which is the topic of the present work. Here, we explore how well the dynamics of methyl groups in protein side-chains, as accessible via NMR relaxation spectroscopy, can be predicted by MD simulations. Ubiquitin is used as a test system for this purpose, because it is one of the bestcharacterized small proteins, including the availability of NMR data. In addition to generalized order parameters S2 and associated correlation times τf of methyl groups, we focus on the comparison of relaxation rates and spectral densities from MD and NMR. As detailed in the Methods section, we calculate
INTRODUCTION Molecular dynamics (MD) simulations have become an established technique to study structural, dynamic, and thermodynamic properties of biomolecular systems, thanks to steady improvements of algorithms, potential energy functions, and hardware.1 A particularly powerful combination is to use MD simulations in conjunction with NMR spectroscopy,2,3 as NMR spin relaxation techniques can probe dynamic motions of biomolecules on a broad range of time scales in a site-specific (i.e., atom- or residue-resolved) manner.4−8 Interpreting the motional parameters derived from NMR, such as generalized order parameters S2 and associated correlation times τf, requires a model for the motions, which in turn can be provided by MD simulations. However, for a sound evaluation of such motional models, e.g., for the motions of backbone or side-chain bond vectors, it is important to quantitatively link MD simulations to NMR experiments. This aspect is especially relevant also for establishing relations between order parameters and conformational entropy.9−15 In general, the predictive power of MD simulations crucially depends on the conformational space explored within the given simulation time (sampling) and the accuracy of the underlying potential energy function (force field). Advances in MD force fields16−22 have enabled a nearly quantitative interpretation of protein backbone dynamics measured by NMR relaxation in terms of squared generalized order parameters of amide NH groups, S2NH, as derived from the NMR data using the Lipari− Szabo (LS) formalism.23,24 NMR relaxation data has been used not only to validate MD force fields a posteriori, but also as target data for optimizing the potential energy functions for © 2018 American Chemical Society
Received: March 22, 2018 Revised: April 24, 2018 Published: April 25, 2018 5038
DOI: 10.1021/acs.jpcb.8b02769 J. Phys. Chem. B 2018, 122, 5038−5048
Article
The Journal of Physical Chemistry B
protein side-chain methyl groups from MD simulations. We found that an accurate description of the methyl dynamics in the simulations with the AMBER ff99SB*-ILDN force field18−20 (and related AMBER parameter sets) required adjusting the energy barriers of CH3 spinning in all methylcontaining amino acids except THR. This reparametrization was done on the basis of CCSD(T) coupled cluster quantum chemical calculations of isolated dipeptides. Our MD simulation approach enables the prediction of relaxation rates and spectral densities, as well as of LS motional parameters S2 and τf, without resorting to any information from NMR experiments. Our procedure is inspired by the previous work of Bremi an co-workers27 and Skrynnikov and co-workers,28 with one key difference being that no experimental data are required for the overall rotational correlation time of the protein τR. Instead, extended MD simulations with the TIP4P/2005 water model38 directly yield rotational correlation times that agree with experiments.
the time-correlation functions (TCFs) of the methyl C−H bond vectors from the MD simulations C(t ) = ⟨P2[μ ⃗ (0) ·μ ⃗ (t )]⟩
(1)
with the second order Legendre polynomial P2[x] = (3x − 1)/ 2, the unit vector of the bond μ⃗(t), and ⟨···⟩ indicating the average over all time step differences t. The total TCF C(t) describes the reorientation motion of the bond in the laboratory frame via the angle θ(t) formed by a bond vector during the time interval t (the scalar product of the unit vectors is cos θ, i.e., P2[x] = P2[cos θ(t)]). The spectral density is the Fourier transform of the TCF 2
J(ω) =
∫0
∞
C(t )e−iωt dt
(2)
The total TCF can be factorized into internal and overall TCFs, C(t) = Cint(t)CO(t). The assumption behind this separation is that the two types of motions are independent. This is always the case if overall tumbling is sufficiently slower than the internal dynamics.33 However, the separation assumption can be valid even for residues that undergo internal motions on the same time scale as tumbling.34 For an isotropically tumbling molecule with a global rotational tumbling correlation time τR, the overall motion can be described by a single-exponential decay, CO(t) = e−t/τR. For small proteins in aqueous solution close to room temperature, τR is on the order of a few nanoseconds. The internal dynamics of backbone vectors in structurally rigid elements of proteins are typically much faster, simplifying their interpretation in terms of LS models. However, for methyl groups, slow internal dynamics due to jumps between side-chain rotamer wells might occur on similar (or even slower) time scales as overall tumbling, which can complicate matters. The square of the generalized order parameter is defined as the long-time limit of Cint(t) S2 = lim C int(t ) t →∞
■
METHODS Quantum Chemical Calculations. To calculate the adiabatic energy barriers of methyl group rotation (spinning around its 3-fold axis), relaxed potential energy surface (PES) scans were performed for isolated blocked dipeptides, ACE-XNME, of the methyl-containing amino acids (X = ALA, MET, THR, LEU, ILE, VAL). The M06-2X density functional39 was used in conjunction with the correlation-consistent double-ζ ccpVDZ basis set.40,41 For selected structures (in the vicinity of the minima and maxima), single-point energies were calculated at the CCSD(T)/cc-pVTZ coupled cluster level. The geometry optimizations were carried out by constraining the dihedral angle describing the methyl group spinning around the C− Cmethyl bond at 5° increments and minimizing the potential energy in the orthogonal degrees of freedom. Due to the 3-fold symmetry of the methyl group, the total scan range was restricted to 120°. The dihedral angle scans were initiated from structures that were geometry-optimized without constraints in a previous step, i.e., from a minimum on the PES. The blocked dipeptide structures used in these quantum chemical calculations were created with PyMOL, setting the backbone ϕ and ψ dihedral angles to initial values of −135° and +135°, respectively, thus representing an extended backbone conformation. The χ1 side-chain dihedral angle was initially set to −60° (i.e., gauche−) for all amino acids except for VAL and THR, where χ1 was set to 180° (trans) and +60° (gauche+), respectively. These values for χ1 were chosen because they correspond to the most populated rotamer state for each of the investigated side-chains.42 The other χ angles in the side-chains (χ2 and χ3, if present) were set to 180°, i.e., the trans conformation. To obtain the PES for methyl rotation with the force field, energy minimizations of the blocked dipeptides in vacuo were carried out using a conjugate gradient algorithm implemented in GROMACS version 5.0.6.43 A GROMACS version compiled in double precision was used for these energy minimizations. We used the AMBER ff99SB*-ILDN18−20 force field in this study, including the modified backbone parameters of Best and Hummer.19 The methyl rotation dihedral angle terms were introduced in the ff99 parameter set44 and carried over to subsequent parameter sets, such as ff99SB18 and ff99SBILDN,20 ff03,45 and ff14SB.46 Hence, our modifications should be transferable to all these AMBER force fields. In the force field energy minimizations in vacuo, no cut-offs were used for the nonbonded interactions. To increment the dihedral angle in
(3)
2
S describes the amplitude of internal reorientation motions and varies between 0 for a completely unrestricted motion of the bond to 1 for a bond vector that is rigid in the molecular frame. Under the assumption that the motion of the C−Cmethyl bond vector that connects the CH3 group to the rest of the side-chain is independent from the spinning of the methyl group around its 3-fold C3v symmetry axis, the internal TCF Cint(t) can be further separated, Cint(t) = Caxis(t)·CCH3(t); the corresponding order parameter of the methyl symmetry axis is coined S2axis. Realistically capturing the correlation times of methyl group spinning by MD simulation proved to be challenging in the past,35−37 and efforts were made to overcome this issue by reducing the methyl rotation barriers for alanine residues in the CHARMM22 force field.35,36 A different approach was taken by Showalter and co-workers,30 who parametrized the methyl group rotation TCFs according to CCH3(t) = 1/9 + 8/9e−t/τCH3; i.e., CCH3(t) decays with a characteristic correlation time τCH3 from 1 to 1/9, which is the theoretical value due to complete motional averaging around the 3-fold axis, assuming ideal tetrahedrality. Showalter and coworkers treated τCH3 as a fit parameter to minimize the difference between the spectral densities from the MD simulations and the NMR experiments. Here, instead of using such a parametrization that involves fitting against available NMR relaxation data, we directly calculate the reorientational TCFs of the C−H bond vectors in 5039
DOI: 10.1021/acs.jpcb.8b02769 J. Phys. Chem. B 2018, 122, 5038−5048
Article
The Journal of Physical Chemistry B
Furthermore, we repeated all NPT simulations without using bond-length constraints on the protein; in these simulations, the equations of motion were integrated with 1 fs (instead of 2 fs) time steps. The results are highly similar to the ones obtained with constraints. Calculation of Order Parameters, Tumbling Times, and Relaxation Rates. Time-correlation functions (TCFs) for the three Cmethyl−Himethyl (i = 1, 2, 3) vector orientations for all methyl groups of all 10 trajectories were calculated up to a maximum lag-time of 50 ns. Internal TCFs were calculated after aligning all trajectory frames to the initial structure. The individual TCFs of the three C−H vectors of the same methyl group were averaged, as were the TCFs from the 10 individual simulations. We used a spectral density approach for calculating the methyl axis order parameter S2axis. In this approach, the internal TCFs of the Cmethyl−Hmethyl vectors (i.e., after removing overall tumbling) were fitted with six exponentials plus an offset27
5° steps, as was done in the quantum chemical calculations described above, a stiff harmonic dihedral angle restraining potential with a force constant of 20 000 kJ/(mol rad2) was used. The contribution of this dihedral angle restraint to the potential energy was less than 0.1 kJ/mol for any of the final energy-minimized structures but was nevertheless excluded. As described in the Results and Discussion section, to reproduce the CCSD(T) potential energy barriers as close as possible with the force field, the force constants (amplitudes) of the dihedral angle potential energy terms describing rotation of the methyl groups, kdih, were adjusted. MD Simulations. All MD simulations were carried out with GROMACS version 5.0.6.43 The X-ray structure of ubiquitin (PDB code 1UBQ)47 was used as the starting structure for the simulations. Crystallographic water molecules were kept. The protein was centered in a cubic box with a minimum distance between protein and the box boundary of 1.2 nm. The system was solvated with 16 235 TIP4P/200538 water molecules. Sodium and chloride ions were added at a concentration of 0.15 mol/L to yield an overall neutral simulation system. Prior to MD simulation, the system was energy-minimized (steepest descent) and equilibrated in the NPT ensemble for 200 ps with harmonic position restraints on all protein heavy atoms, with force constants of 1000 kJ/(mol nm2). The AMBER ff99SB*ILDN protein force field18−20 was used. To keep the temperature constant at 300 K, the velocity rescaling thermostat of Bussi and co-workers48 was applied, with coupling time constants of τT = 0.1 ps. Isotropic Parrinello− Rahman pressure coupling was used to maintain constant 1 bar pressure, with a 2 ps coupling time constant and a compressibility of 4.5 × 10−5 bar−1. The SETTLE49 and LINCS50 constraint algorithms were applied to constrain the internal degrees of freedom of water molecules and the bonds in the protein, respectively, allowing for integrating the equations of motion with 2 fs time steps using the velocity Verlet integrator.51 Lennard-Jones (6, 12) interactions were smoothly shifted to zero at a cut-off distance of 1.0 nm; this cut-off distance was also used for the short-range Coulomb interactions. Analytical dispersion corrections were added to energy and pressure to correct for the truncation of the Lennard-Jones interactions. Long-range Coulomb interactions were treated with the particle mesh Ewald (PME) method52 with a 0.12 nm grid spacing and cubic spline interpolation. Finally, ten 100 ns production MD simulations were carried out in the NPT ensemble by choosing different random seeds for generating initial velocities from a Maxwell−Boltzmann distribution at 300 K. Coordinates were saved to disk every 1 ps. Additionally, we performed simulations in the NVE ensemble. Prior to the production runs, we equilibrated the system in the NVT ensemble for 100 ps with the same simulation parameters as described above. We started this equilibration from the final structure of one of the previous NPT simulations. From this 100 ps NVT equilibration, we initiated ten 100 ns NVE simulations using snapshots taken every 10 ps. To avoid temperature drift, the NVE simulations were carried out with a GROMACS version that was compiled in double precision; the number of LINCS iterations to correct the rotational bond lengthening was increased to 8, and the maximum allowed error for pair interactions per particle caused by the Verlet buffer was reduced to 5 × 10−8 kJ/(mol ps). This ensured total temperature drifts of less than 0.2 K over the course of the 100 ns NVE simulations.
6
C int(t ) =
∑ A i e −t / τ
i
2 + S long
(4)
i=1
where S2long would correspond to the square of the generalized order parameter in the long-time limit. The different decay times could represent different kinds of motions, although these parameters do not necessarily need to be associated with a physical meaning. In the fit to eq 4, we used ∑i6= 1Ai + S2long = 1, 0 ≤ Ai ≤ 1, 0 ≤ S2long ≤ 1, 0 ≤ τ1 ≤ 10 ps, 0 ≤ τ2 ≤ 50 ps, 0 ≤ τ3 ≤ 200 ps, 0 ≤ τ4 ≤ 500 ps, 0 ≤ τ5 ≤ 1000 ps, and τ6 ≥ 0. To obtain the total TCF that also includes overall tumbling, the resulting internal TCF was multiplied with a singleexponential function28 C(t ) = C int(t )CO(t ) = C int(t )e−t / τR
(5)
using the isotropic tumbling time τR = 4.37 ns obtained from the diffusion tensor of the backbone N−H vectors. To that end, we determined the residue-specific N−H tumbling times from our MD simulations using the approach described by Maragakis and co-workers,53 which includes a direct fit of the total TCFs (i.e., including overall tumbling) to an LS model eff
2 2 C(t ) = S NH e−t / τR,i + (1 − S NH )e−t / τc,i
(6)
with effective correlation time τeff c = (τRτf)/(τR + τf). The fit parameters are the correlation times τR,i and τf as well as S2NH. These residue-specific global tumbling times τR,i were used together with the energy-minimized X-ray structure of ubiquitin to determine the isotropic rotational tumbling time with the program QUADRIC, downloaded from the Web site of the Palmer group (www.palmer.hs.columbia.edu). The value of τR = 4.37 ns obtained from our MD simulations agrees with available NMR data on ubiquitin, which are 3.5 ns,54 4.03 ns,55 4.09 ns,56 and 4.16 ns.57 The total TCF (eq 5) was converted into a spectral density via 6
J(ω) =
∑ i=1
Ai τi i 2 1 + (ωτeff )
+
2 S long τR
1 + (ωτR )2
(7)
with τeff i = (τiτR)/(τi + τR). Then, the values of J(0), J(ωD), and J(2ωD) were fitted to the two-parameter Lipari−Szabo (LS2) model23,24 5040
DOI: 10.1021/acs.jpcb.8b02769 J. Phys. Chem. B 2018, 122, 5038−5048
Article
The Journal of Physical Chemistry B Table 1. Potential Energy Barriers (in kJ/mol) of Methyl Group Spinning in Blocked Dipeptidesa original FF reparametrized FF CCSD(T)
ALA
MET
THR
VAL
LEU
ILE
15.5 14.2 14.2
9.0 7.2 7.1
11.0 11.0 11.4
18.4/17.3 13.1/12.1 14.0/11.5
16.8/16.2 13.9/13.3 14.1/12.9
17.4/13.5 12.4/10.7 12.2/10.7
a
For VAL and LEU, the values before and after the slash are for the trans/gauche- methyl groups, respectively. For ILE, the values before and after the slash are for the Cγ/Cδ methyl groups, respectively. THR was not reparametrized.
J(ω) =
1 2 ⎞ ⎛ 1 2 1 − 9 Saxis τceff ⎟ 2 ⎜ 9 SaxisτR + 5 ⎜⎜ 1 + (ωτR )2 1 + (ωτceff )2 ⎟⎟ ⎠ ⎝
(
■
RESULTS AND DISCUSSION Methyl Group Rotation Barriers in Dipeptides. Table 1 summarizes the results from our quantum chemical calculations of the potential energy barriers of methyl group spinning in amino acid side-chains, as obtained for isolated blocked dipeptides. We found that, for all side-chains except THR, the barriers in the AMBER ff99SB*-ILDN force field, which originate from the ff99 parameter set,44 were too high compared to the reference values obtained from CCSD(T)/ cc-pVTZ quantum chemical calculations. As described in detail below, this overestimation sensitively affects the deuterium relaxation rates and spectral densities obtained from the force field MD simulations of ubiquitin. Interestingly, the barriers are overestimated by the force field to different extents for the different methyl groups. For example, in ALA, the barrier is too high by 1.3 kJ/mol in the original force field, whereas for ILECγ, the barrier is overestimated by 5.2 kJ/mol (Table 1 and Figure 1; the potential energy curves for all side-chains are
)
(8)
as this model was also used to interpret the experimental NMR relaxation data of ubiquitin.58 The fit was carried out according to the following procedure: 1. Grid search within S2axis ∈ [0, 0.01, ..., 2S2long] and τf ∈ [0 2 eff 6 ps, 1 ps, ..., 2τeff fit ]. Here, we used Slong and τfit = ∑i = 1Aiτi as starting points for the fitting parameters. 2. Direct (gradient-based) fit of S2axis and τf with the results from the previous grid search as starting points. This yielded the final parameter pair (S2axis; τf) of the fit. At every step of the grid search or gradient-based minimization, we determined Jmodel(0), Jmodel(ωD), and Jmodel(2ωD) with eq 8, and calculated the longitudinal and transverse relaxation rates, R(Dz) and R(Dy), respectively, as well as χ2 via R(Dz ) =
2 2 3 ⎛ qQ e ⎞ ⎜ ⎟ [J(ωD) + 4J(2ωD)] 16 ⎝ ℏ ⎠
R(Dy) =
2 2 1 ⎛ qQ e ⎞ ⎜ ⎟ [9J(0) + 15J(ωD) + 6J(2ωD)] 32 ⎝ ℏ ⎠
(9)
(10)
χ2 =
⎛ R model(i) − R analytical(i) ⎞2 ⎟⎟ ∑ ⎜⎜ R error(i) ⎠ i ∈ {Dz , Dy} ⎝
Figure 1. Potential energies along the methyl group spinning dihedral angle in alanine blocked dipeptide. Relaxed potential energy surface scans were carried out at the M06-2X/cc-pVDZ level (magenta curve). Single-point energies of these optimized geometries were calculated with CCSD(T)/cc-pVTZ (black). Energies from force field calculations with the original and reparametrized AMBER ff99SB*-ILDN parameter sets are plotted in red and blue, respectively.
(11)
The fitting parameter pair (S2axis; τf) for eq 8 was determined by minimizing χ2. The errors Rerror(i) are the standard errors of the mean of the 10 relaxation rates Rmodel(i) obtained from the 10 individual MD simulations. Ranalytical(i) are the relaxation rates directly derived from J(0), J(ωD), J(2ωD) via eq 7. The Ranalytical values are reported below, because invoking the LS2 model is not required here. Our above approach is based on the separation of the total TCF into internal and overall motions, C(t) = Cint(t)CO(t), and assumes a certain analytical form of both CO (singleexponential) and Cint (6-exponential plus offset, eq 4). This allows one to gauge the contributions of internal and overall motions, and avoids issues linked to noise in the TCFs at long lag-times. To test the validity of this approach, we calculated the raw TCF directly from our MD simulations, i.e., without any model for the internal or overall motion. We compare this total TCF to C(t) = Cint(t)CO(t), where CO(t) is the overall TCF described by the single-exponential fit, and Cint(t) is the raw internal TCF directly obtained from the MD simulations (i.e., without assuming a certain functional form), or from the fit to eq 4. Figure S1 in the Supporting Information shows that the TCFs are correctly captured up to lag-times of 10 ns, demonstrating the validity of the applied assumptions.
shown in Figure S2 in Supporting Information). To resolve this discrepancy, we reparametrized the force field in a methyl group-specific manner. The force constants of the dihedral angle potential energy function, kdih, were lowered such that the reparametrized force field matched the potential energy barriers from the CCSD(T) calculations (Figure 1 and Table 2). Reassuringly, we found that the corresponding methyl groups in different side-chains required the same lowering of kdih (compare Cγ in VAL and ILE, and Cδ in LEU and ILE, respectively, in Table 2). In the quantum chemical calculations, the geometry optimizations were carried out with the M06-2X density functional. For all investigated side-chains apart from MET, we found the M06-2X/cc-pVDZ potential energy barriers to be slightly higher than in the CCSD(T) coupled cluster calculations (Figure S2 in Supporting Information). To check the influence of the basis set size on the energies, we used 5041
DOI: 10.1021/acs.jpcb.8b02769 J. Phys. Chem. B 2018, 122, 5038−5048
Article
The Journal of Physical Chemistry B
internal motions, τf, and then turn to deuterium relaxation rates and spectral densities. We calculated the LS parameters S2axis and τf from our MD simulations with the spectral density mapping approach described in the Methods section, and compare our MD results to experimental values from Liao and co-workers.58 Figure 2A and Table 3 compare the NMR order parameters to
Table 2. Difference of the Force Constants kdih of the Dihedral Energy Terms for Methyl Group Spinning between the Reparametrized and the Original Force Fielda methyl group β
ALA C MET Cϵ VAL Cγ LEU Cδ ILE Cγ ILE Cδ
Δkdih [kJ/mol] −0.06964 −0.31380 −0.30220 −0.16270 −0.30220 −0.16270
In the original force field, the kdih values of the H−C−C−Hmethyl and C−C−C−Hmethyl dihedral angles are 0.62760 and 0.66944 kJ/mol, respectively; for the C−S−C−Hmethyl dihedral angle in MET, kdih = 1.39467 kJ/mol in the original force field. THR was not reparametrized.
a
ethane, as this small molecule renders a systematic comparison computationally readily feasible. We calculated the adiabatic potential energy barriers for methyl rotation in ethane and compared the energies obtained from CCSD(T) calculations with the cc-pVTZ triple-ζ basis to the cc-pVQZ and aug-ccpVQZ quadruple-ζ basis sets. For ethane in vacuo, M06-2X/ccpVDZ yields a potential energy barrier of 12.3 kJ/mol, compared to 11.7, 11.4, and 11.4 kJ/mol from the CCSD(T) calculations with the cc-pVTZ, cc-pVQZ, and aug-cc-pVQZ basis sets, respectively. Assuming that these results are transferable to the amino acid hydrocarbon side-chains, we conclude that the CCSD(T) calculations with the cc-pVTZ basis set yield highly accurate energetics. The required lowering of kdih can be rationalized a priori based on the following simple considerations.36 In the AMBER force field, each set of four covalently connected atoms contributes one dihedral energy term to the force field potential energy. Hence, for all methyl-containing side-chains except MET, nine such terms contribute to the overall methyl group rotation barrier. The potential energy term for the dihedral angle ϕ describing the spinning of the methyl group around its 3-fold axis has the general form Vdih = kdih [1 + cos(3ϕ)]. Thus, each term contributes ca. 2kdih to the methyl rotation barrier.36 For example, for ALA, reducing kdih by 0.069 64 kJ/mol (Table 2) should reduce the barrier by ca. 2 × 9 × 0.069 64 kJ/mol = 1.25 kJ/mol, in close agreement with the actual reduction of 1.3 kJ/mol (from 15.5 to 14.2 kJ/mol, Table 1). Similar considerations for the methyl groups in the other side-chains, and taking into account that for MET there are only three (instead of nine) dihedral terms, readily yields the required adjustments of kdih. In addition to the dihedral terms, the rotation barrier in the force field also depends on the nonbonded (electrostatic and van der Waals) interactions. However, because changes to the nonbonded interactions might have undesired effects beyond the dihedral barriers, we decided to adjust the methyl rotation barriers only through kdih, hence leaving the rest of the force field unaffected. Due to the isolated nature of the rotation of a methyl group at the terminus of a hydrocarbon chain, no effects related to couplings to other energy terms need to be taken into account. Methyl Group Dynamics in Ubiquitin. To investigate the influence of the force field reparametrization on the methyl group dynamics in proteins we used ubiquitin as a model system, an experimentally very well-characterized small protein that is frequently used in simulation studies. We first focus on methyl axis order parameters, S2axis, and correlation times of fast
Figure 2. Methyl axis order parameter S2axis (A) and corresponding correlation time τf (B) from MD simulations with the AMBER ff99SB*-ILDN force field without (red dots) and with (blue triangles) adjusted methyl rotation barriers are correlated to the experimental values from Liao and co-workers.58 The statistical uncertainties are the standard errors of the mean from the 10 individual MD simulations.
Table 3. Pearson and Spearman Correlation Coefficients (RP and RS, Respectively) and Root Mean Square Deviation (RMSD) of Methyl Order Parameter S2axis between NMR and MD simulation
RP
RS
RMSD
original FF (NPT) reparametrized FF (NPT) reparametrized FF (NVE)
0.91 0.85 0.85
0.91 0.88 0.86
0.10 0.13 0.13
the MD values for both force fields, the original AMBER ff99SB*-ILDN and the one with reparametrized barriers for methyl spinning. Correlation coefficients and root mean square deviations (RMSDs) are summarized in Table 3. The results obtained from our simulations in the NVE ensemble are very similar and consistent with the NPT simulations, showing that the thermostat used in the latter does not have a large influence on the results (see Figure S3 in Supporting Information). Our results show that, concerning S2axis, both force fields agree equally well with the NMR experiment, with an RMSD in S2axis of about 0.1. This agreement is considered encouraging, given that obtaining accurate methyl order parameters from MD simulations is considered much harder than for backbone amides, for which RMSDs in S2NH of about 0.05 between MD and NMR are typically achieved. For ubiquitin, previous MD simulation studies reported RMSDs from the NMR-derived54 5042
DOI: 10.1021/acs.jpcb.8b02769 J. Phys. Chem. B 2018, 122, 5038−5048
Article
The Journal of Physical Chemistry B S2axis of 0.125 (Long and co-workers, obtained from 1 μs simulations with the AMBER ff99SBnmr1-ILDN force field59), ca. 0.15 (Bowman,60 obtained from 1 μs simulations with the AMBER ff03, ff99SB-ILDN, and CHARMM27 force fields), and 0.104 (Kasinath and co-workers,15 obtained from a 260 ns simulation with the CHARMM27 force field). To assess the level of agreement that could ideally be expected from the MD simulations, we correlated the methyl axis order parameters obtained from two different deuterium relaxation experiments of ubiquitin, carried out by Lee and coworkers54 and Liao and co-workers.58 Interestingly, the RMSD in S2axis between these two experimental NMR data sets is 0.11 and thus comparable to the deviation between MD and NMR, although the Pearson correlation coefficient of 0.98 is very high.58 The differences between these two NMR data sets seem to be systematic, with the S2axis values derived from Lee and coworkers being on average 14% higher than the values of Liao and co-workers.58 Nevertheless, we consider the RMSD of about 0.1 between our MD simulations and the NMR experiments as fairly good, although there is clearly still room for improvements. Restricting the comparison of MD to NMR to correlation coefficients and RMSDs might not completely unveil the uncertainty of the predictions, and further analysis might be warranted. For example, as typical S2axis values vary substantially between different types of methyl groups, it could be interesting to compare how well the MD simulations do in predicting the site-to-site variation, as compared to a null model that simply assigns the average value of each methyl type to all methyls of that type. However, as the main focus of this work is not on S2axis, but on relaxation rates and spectral densities (see below), we leave such an analysis to future work. Possible reasons for the remaining discrepancies between NMR and MD are shortcomings of the force field and incomplete sampling. Sampling the rotameric transitions around all χ angles of a side-chain is expected to have a strong influence on S2axis, with low S2axis values indicating population of more than a single rotamer well.61−63 Insufficient sampling in our MD simulations would thus be expected to lead to a systematic overestimation of S2axis, which indeed seems to be the case for the low-order parameters (S2axis < 0.4), which are all above the diagonal in Figure 2A. Fully addressing the sampling issue is beyond the scope of the present work. However, it is conceivable that the populations of and/or interconversion between rotamer wells in the ensemble present in the NMR experiments differs from that in the MD simulations, as in the latter case we are limited to a few repeats of limited time (here, ten 100 ns simulations) and with a single molecule. The noticeable deviations in S2axis between MD and NMR that are observed for VAL70-Cγ1,γ2 and for ILE13-Cγ2 hint in this direction. For VAL70, which is part of a hydrophobic patch on the ubiquitin surface, our MD simulations predict S2axis for Cγ1/ Cγ2 of 0.59/0.68, respectively, which are considerably higher than the NMR deuterium relaxation values of 0.38/0.34.58 Likewise, for ILE13-Cγ2, MD predicts S2axis = 0.75, whereas NMR relaxation yields 0.52. VAL70 and ILE13-Cγ2 have been identified by NMR as residues in which the C−Cmethyl vector undergoes slow motions linked to rotameric jumps.58 Indeed, as shown in Figure S4 in Supporting Information, VAL70 undergoes rather slow rotameric transitions around the χ1 sidechain dihedral angle in the MD simulations, with dwell times in the rotamer wells that are often longer than the global tumbling time. The most populated state of the VAL70 side-chain is
gauche− whereas the gauche+ conformer is hardly populated, in qualitative agreement with rotamer populations determined by NMR from three-bond scalar couplings and dipolar couplings.64 However, the rate of rotameric transitions in VAL70 might be underestimated in the MD simulations, and rotameric averaging is thus incomplete within τR. This interpretation is further supported by analyzing Cint(t) from MD in terms of its long-time limit (eq 3). Averaging Cint(t) between 25 and 50 ns yields S2axis = 0.39 for VAL70-Cγ2, in close agreement with the value of 0.34 from NMR relaxation. Consistent with these findings, the order parameters of ILE13Cγ2 and both VAL70 methyl groups were also overestimated in the MD simulations of Liao and co-workers.58 Figure 2B shows that, compared to S2axis, the improvements due to force field reparametrization are much more evident in the correlation time of fast internal dynamics τf, which is severely overestimated by the original force field due to too slow methyl spinning. The RMSD with respect to the NMR values drops from 77.9 ps for the original force field to 20.3 ps for the reparametrized force field. However, instead of discussing this in much detail, we would like to focus on relaxation rates and spectral densities. Methyl order parameters and associated τf’s are determined with the formalism of Lipari and Szabo,23,24 which assumes that the internal dynamics and the overall tumbling of the protein are independent. This assumption is fulfilled if these two processes occur on different time scales; i.e, the fast internal dynamics are separated from the slow overall tumbling. This is usually the case for the dynamics of backbone N−H vectors, where the fast internal motions happen on the picosecond time scale while the typical tumbling time of a small protein is a few nanoseconds. In the case of methyl dynamics, fast (ps) librational motions are accompanied by jumps of the methyl group between rotamer states. These jumps around side-chain χ dihedral angles can happen on a broad range of time scales that might even be slower than global rotational tumbling, hence complicating interpretations in terms of the (simplified) Lipari−Szabo models. Indeed, for a fraction of the methyl-bearing sidechains in HIV protease, ubiquitin, and TNfn3, the generalized methyl order parameters derived from experimentally measured 3 J coupling constants, which include motions over the entire range up to milliseconds, were systematically lower than the corresponding S2 values derived from deuterium relaxation.64,65 This indicates rotameric averaging on time scales that are too slow to be observable in conventional relaxation experiments. It will be interesting to study the nature of such motions and how they affect S2 in future work. In any case, it is thus beneficial not to restrict the comparison of MD and NMR only to S2 and τf. The relaxation rates are directly accessible to NMR relaxation experiments and do not depend on model assumptions. We thus calculated the relaxation rates R(Dy) and R(Dz) from our MD simulations and compare them with experiment (Figure 3, Table 4, and Table S1 in Supporting Information). Reparametrization of the methyl rotation barriers in the force field yields relaxation rates that are in much better agreement with the experiment than the original force field, with a reduction of the RMSD by about a factor 1.3 and 2.2 for R(Dy) and R(Dz), respectively. The stronger improvement for R(Dz) is understandable, because it reflects more on fast (high frequency) motions than R(Dy) does, so adjusting the very fast (sub-nanosecond) methyl spinning rates in the force field is expected to have a more pronounced effect. As observed for the order parameters, the NPT and NVE simulations yield highly 5043
DOI: 10.1021/acs.jpcb.8b02769 J. Phys. Chem. B 2018, 122, 5038−5048
Article
The Journal of Physical Chemistry B
Figure 3. Relaxation rates RMD(Dy) (A) and RMD(Dz) (B) from the simulations with the original AMBER ff99SB*-ILDN force field (red dots) and the one with reparametrized methyl spinning barriers (blue triangles) are correlated to the experimental relaxation rates from Liao and co-workers.58 The statistical uncertainties are the standard errors of the mean from the 10 individual MD simulations.
Table 4. Root Mean Square Deviation (RMSD) of Relaxation Rates R(Dy) and R(Dz) between MD and NMR simulation
RMSD in R(Dy) [s−1]
RMSD in R(Dz) [s−1]
original FF (NPT) reparametrized FF (NPT) reparametrized FF (NVE)
21.3 15.8 16.6
21.0 9.3 9.4
similar results also for the relaxation rates (Table 4 and Table S2 in Supporting Information). Likewise, the relaxation rates obtained from MD simulations carried out without bond-length constraints on the protein (see Methods section) yield very similar results, see Table S3 in Supporting Information. The agreement between MD and NMR critically depends on the spectral densities. We calculated the spectral density of every methyl group from our MD simulations and fitted the values at J(0), J(ωD), and J(2ωD) with the LS2 model (see Methods section), which was also used in the experimental NMR study.58 Figure 4 compares these MD-based spectral densities to the experimental curves, which were constructed from the reported values for S2axis and the correlation times of fast internal motions τf, and of global tumbling τR = 5.0 ns.58 We chose one representative methyl group for every amino acid type; the spectral densities of all methyls are shown in Figure S5 in Supporting Information. Figure 4 shows that adjusting the methyl rotation barrier in the force field improves the spectral densities in comparison to the original force field, especially at high frequencies. The improvement is particularly pronounced for ILE, LEU, and VAL methyl groups, for which the potential energy barriers were reduced most drastically (see Tables 1 and 2). For MET, conclusions are difficult, because there is only one MET in ubiquitin, which in addition is the first residue of the polypeptide chain (MET1).
Figure 4. Spectral densities of selected methyl groups from the simulations with the original AMBER ff99SB*-ILDN force field (left panels) and the reparametrized force field (right panels). The spectral density from MD is shown in red, and the values at 0, ωD, and 2ωD are indicated by red triangles. The LS2 fit to these three data points is shown as a dashed cyan line, and the fitting parameters S2axis and τf are reported. The spectral densities from NMR, obtained with a global tumbling time of 5.0 ns in D2O,58 are shown in black. The black curves in the insets show the corresponding NMR spectral densities for a global tumbling correlation time of 4.16 ns57 from 15N relaxation in 90% H2O/10% D2O.
Figure 4 compares the spectral densities from NMR experiments in 99.9% D2O58 with those from MD simulations in 100% H2O. The TIP4P/2005 water model used in our MD simulations yields a global tumbling correlation time of τR = 4.37 ± 0.27 ns. The statistical uncertainty of ca. 6%, estimated from the standard error of the mean of the 10 individual MD simulations, shows that the sampling suffices to compute τR with reasonable precision. However, D2O has a higher viscosity than H2O, and hence τR is somewhat longer. In their NMR deuterium relaxation study, Liao and co-workers58 used the ratio of the viscosities of D2O and H2O to obtain the global 5044
DOI: 10.1021/acs.jpcb.8b02769 J. Phys. Chem. B 2018, 122, 5038−5048
Article
The Journal of Physical Chemistry B tumbling correlation time of 5.0 ns from the value of 4.16 ns measured in previous 15N relaxation experiments in 90% H2O/ 10% D2O.57 As the value of τR influences the spectral density, especially in the low-frequency range, the insets in Figure 4 show the NMR spectral densities that would be obtained with a global tumbling time of 4.16 ns in H2O (assuming that S2axis and τf are unchanged). As expected from the similar global tumbling times of 4.37 ns (MD) and 4.16 ns (NMR), the agreement between the MD and these hypothetical NMR spectral densities is even better, especially close to ω = 0. A similar result would have been obtained if we would have scaled our MD-derived τR with the D2O/H2O viscosity ratio to mimic the slower tumbling in D2O in an a posteriori manner. We do not apply such scaling here, because our aim is to make predictions based on MD simulations only, by taking the truncation of the internal TCFs due to overall tumbling explicitly into account. Figure 4 shows that this successfully works. Finally, to focus on fast (sub-nanosecond) motions on the time scale of methyl group spinning, which are not easily extractable from Figure 4, we generated “synthetic” NMR data by calculating the internal TCFs of the C−Hmethyl bond vector reorientation from the NMR-derived LS parameters58 and compare them with the internal TCFs obtained from our MD simulations (eq 4). The results for the five representative methyl groups are shown in Figure 5, and for all methyl groups in Figure S6 in Supporting Information. For all amino acid types, the adjusted force field yields TCFs that are in good agreement with NMR. We close this discussion by comparing our present approach to alternative ones in the literature that also take overall tumbling explicitly into account. In our all-atom MD simulations of ubiquitin, the overall rotational tumbling of the protein is correctly captured with the TIP4P/2005 water model; a combined MD sampling of 1 μs (i.e., about 200 times the rotational correlation time τR of ubiquitin, every individual 100 ns simulation exceeds τR by about a factor 20) suffices to obtain statistically precise τR values. Alternative water models that yield realistic diffusional behavior, such as SPC/Eb, were devised and successfully used to reproduce site-specific spectral densities of backbone amide bonds in ubiquitin and protein G.66 An alternative route was followed by Chen and coworkers,67,68 who (i) also focused exclusively on backbone amide bond motions (and not on side-chains), and (ii) described overall rotational tumbling with computationally very efficient coarse-grained MD simulations with the MARTINI force field. However, this approach is not generally applicable, because the unavoidable elastic network potentials required to restrain the protein structure in these coarse-grained simulations impede the extension to somewhat more flexible proteins. In addition, the coarse-grained water force field cannot be expected to yield realistic rotational diffusion times at different conditions, such as higher temperature or pressure. Recently, Anderson and co-workers69 described an MD simulation approach to predict time-correlation functions of backbone amide and nonmethyl C−H bond vectors that employs rotational velocity rescaling.70 As for the other approaches mentioned above, it would be interesting to extend these methods also to the dynamics of side-chain methyl groups in the future.
Figure 5. TCFs of selected methyl groups. The internal TCFs from MD for the original and the reparametrized AMBER ff99SB*-ILDN force field are shown in red and blue, respectively. The internal TCFs from NMR, calculated from S2axis and τf reported by Liao and coworkers,58 are shown in green.
densities, as well as the associated Lipari−Szabo dynamic parameters S2axis and τf, from the time-correlation functions of C−H bond vector reorientation in protein side-chain methyl groups. To achieve agreement with NMR deuterium spin relaxation experiments, we reparametrized the potential energy barriers of methyl group spinning in the used AMBER ff99SB*ILDN force field, which was done using energies from coupled cluster quantum chemical calculations of isolated dipeptides as reference. The validity of these force field modifications, which also apply to AMBER parameter sets related to ff99SB*-ILDN, was demonstrated by comparing the methyl group dynamics in ubiquitin in MD simulations with NMR relaxation data. Our MD simulation approach does not require NMR data or system-specific adjustable parameters. Hence, MD simulations can be used to predict NMR relaxation data of methyl groups in protein side-chains, also for unknown or experimentally difficult to access proteins. Vice versa, experimental NMR relaxation data
■
CONCLUSIONS In the present work, we describe an MD simulation approach to directly calculate deuterium relaxation rates, spectral 5045
DOI: 10.1021/acs.jpcb.8b02769 J. Phys. Chem. B 2018, 122, 5038−5048
Article
The Journal of Physical Chemistry B can be used to gauge protein force fields and foster continued efforts for their improvement, especially of side-chains. The good agreement between MD and NMR reported in this work increases the confidence into the dynamics seen in the MD simulations, and make it possible to investigate the microscopic origins and physical models that are appropriate for describing the experimental observations in close detail. As discussed, the Lipari−Szabo models suffer from limitations, such as the assumed statistical independence of internal and overall motions, and that only simple models for relaxation can be invoked, which in many cases defies reality given the complexity of side-chain motions. Being able to compare simulations and experiments at the level of relaxation rates has the advantage that these simplifying assumptions are not needed any more.
■
(5) Sapienza, P. J.; Lee, A. L. Using NMR to Study Fast Dynamics in Proteins: Methods and Applications. Curr. Opin. Pharmacol. 2010, 10, 723−730. (6) Mittermaier, A.; Kay, L. E. New Tools Provide New Insights in NMR Studies of Protein Dynamics. Science 2006, 312, 224−228. (7) Palmer, A. G. NMR Characterization of the Dynamics of Biomacromolecules. Chem. Rev. 2004, 104, 3623−3640. (8) Brüschweiler, R. New Approaches to the Dynamic Interpretation and Prediction of NMR Relaxation Data from Proteins. Curr. Opin. Struct. Biol. 2003, 13, 175−183. (9) Akke, M.; Brüschweiler, R.; Palmer, A. G. NMR Order Parameters and Free Energy: An Analytical Approach and its Application to Cooperative Calcium(2+) Binding by Calbindin D9k. J. Am. Chem. Soc. 1993, 115, 9832−9833. (10) Yang, D.; Kay, L. E. Contributions to Conformational Entropy Arising from Bond Vector Fluctuations Measured from NMR-Derived Order Parameters: Application to Protein Folding. J. Mol. Biol. 1996, 263, 369−382. (11) Li, Z.; Raychaudhuri, S.; Wand, A. J. Insights Into the Local Residual Entropy of Proteins Provided by NMR Relaxation. Protein Sci. 1996, 5, 2647−2650. (12) LeMaster, D. M. NMR Relaxation Order Parameter Analysis of the Dynamics of Protein Side Chains. J. Am. Chem. Soc. 1999, 121, 1726−1742. (13) Frederick, K. K.; Marlow, M. S.; Valentine, K. G.; Wand, A. J. Conformational Entropy in Molecular Recognition by Proteins. Nature 2007, 448, 325−329. (14) Petit, C. M.; Zhang, J.; Sapienza, P. J.; Fuentes, E. J.; Lee, A. L. Hidden Dynamic Allostery in a PDZ Domain. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 18249−18254. (15) Kasinath, V.; Sharp, K. A.; Wand, A. J. Microscopic Insights into the NMR Relaxation-Based Protein Conformational Entropy Meter. J. Am. Chem. Soc. 2013, 135, 15092−15100. (16) Best, R. B.; Zhu, X.; Shim, J.; Lopes, P. E. M.; Mittal, J.; Feig, M.; MacKerell, A. D. Optimization of the Additive CHARMM AllAtom Protein Force Field Targeting Improved Sampling of the Backbone ϕ, ψ and Side-Chain χ1 and χ2 Dihedral Angles. J. Chem. Theory Comput. 2012, 8, 3257−3273. (17) Huang, J.; MacKerell, A. D. CHARMM36 All-atom Additive Protein Force Field: Validation Based on Comparison to NMR Data. J. Comput. Chem. 2013, 34, 2135−2145. (18) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of Multiple Amber Force Fields and Development of Improved Protein Backbone Parameters. Proteins: Struct., Funct., Genet. 2006, 65, 712−725. (19) Best, R. B.; Hummer, G. Optimized Molecular Dynamics Force Fields Applied to the Helix−Coil Transition of Polypeptides. J. Phys. Chem. B 2009, 113, 9004−9015. (20) Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved Side-Chain Torsion Potentials for the Amber ff99SB Protein Force Field. Proteins: Struct., Funct., Genet. 2010, 78, 1950−1958. (21) Debiec, K. T.; Cerutti, D. S.; Baker, L. R.; Gronenborn, A. M.; Case, D. A.; Chong, L. T. Further Along the Road Less Traveled: AMBER ff15ipq, an Original Protein Force Field Built on a SelfConsistent Physical Model. J. Chem. Theory Comput. 2016, 12, 3926− 3947. (22) Wang, L.-P.; McKiernan, K. A.; Gomes, J.; Beauchamp, K. A.; Head-Gordon, T.; Rice, J. E.; Swope, W. C.; Martinez, T. J.; Pande, V. S. Building a More Predictive Protein Force Field: A Systematic and Reproducible Route to AMBER-FB15. J. Phys. Chem. B 2017, 121, 4023−4039. (23) Lipari, G.; Szabo, A. Model-Free Approach to the Interpretation of Nuclear Magnetic Resonance Relaxation in Macromolecules. 1. Theory and Range of Validity. J. Am. Chem. Soc. 1982, 104, 4546− 4559. (24) Lipari, G.; Szabo, A. Model-Free Approach to the Interpretation of Nuclear Magnetic Resonance Relaxation in Macromolecules. 2.
ASSOCIATED CONTENT
S Supporting Information *
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.8b02769 A python script to modify the methyl rotation barriers in GROMACS force field topology files (.itp or .top files) has been uploaded to our Web site www.molecular-simulation.org. Tables with relaxation rates, validation of factorization of total TCF into internal and overall TCFs and of their functional forms, potential energy curves for methyl rotation in all side-chains, correlation plot of S2axis from NPT and NVE simulations, plot of the rotameric transitions of the VAL70 side-chain, plots of spectral densities from MD and NMR, and plots of TCFs from MD and NMR (PDF)
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. ORCID
Frans A. A. Mulder: 0000-0001-9427-8406 Lars V. Schäfer: 0000-0002-8498-3061 Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS We thank Vitali Tugarinov for providing the relaxation rates from the NMR experiments described in ref 58, and Kresten Lindorff-Larsen for helpful discussions. This work was supported by Deutsche Forschungsgemeinschaft (DFG) through Cluster of Excellence RESOLV (EXC 1069) and Emmy-Noether grant SCHA 1574/3-1 to L.V.S. The Steinbuch Centre for Computing (SCC), Karlsruhe/Germany, is acknowledged for providing computational resources.
■
REFERENCES
(1) Dror, R. O.; Dirks, R. M.; Grossman, J.; Xu, H.; Shaw, D. E. Biomolecular Simulation: A Computational Microscope for Molecular Biology. Annu. Rev. Biophys. 2012, 41, 429−452. (2) Case, D. A. Molecular Dynamics and NMR Spin Relaxation in Proteins. Acc. Chem. Res. 2002, 35, 325−331. (3) Lindorff-Larsen, K.; Best, R. B.; DePristo, M. A.; Dobson, C. M.; Vendruscolo, M. Simultaneous Determination of Protein Structure and Dynamics. Nature 2005, 433, 128−132. (4) Peng, J. W. Exposing the Moving Parts of Proteins with NMR Spectroscopy. J. Phys. Chem. Lett. 2012, 3, 1039−1051. 5046
DOI: 10.1021/acs.jpcb.8b02769 J. Phys. Chem. B 2018, 122, 5038−5048
Article
The Journal of Physical Chemistry B Analysis of Experimental Results. J. Am. Chem. Soc. 1982, 104, 4559− 4570. (25) Li, D.-W.; Brüschweiler, R. NMR-Based Protein Potentials. Angew. Chem., Int. Ed. 2010, 49, 6778−6780. (26) Li, D.-W.; Brüschweiler, R. Iterative Optimization of Molecular Mechanics Force Fields from NMR Data of Full-Length Proteins. J. Chem. Theory Comput. 2011, 7, 1773−1782. (27) Bremi, T.; Brüschweiler, R.; Ernst, R. R. A Protocol for the Interpretation of Side-Chain Dynamics Based on NMR Relaxation: Application to Phenylalanines in Antamanide. J. Am. Chem. Soc. 1997, 119, 4272−4284. (28) Skrynnikov, N. R.; Millet, O.; Kay, L. E. Deuterium Spin Probes of Side-Chain Dynamics in Proteins. 2. Spectral Density Mapping and Identification of Nanosecond Time-Scale Side-Chain Motions. J. Am. Chem. Soc. 2002, 124, 6449−6460. (29) Best, R. B.; Clarke, J.; Karplus, M. What Contributions to Protein Side-Chain Dynamics are Probed by NMR Experiments? A Molecular Dynamics Simulation Analysis. J. Mol. Biol. 2005, 349, 185− 203. (30) Showalter, S. A.; Johnson, E.; Rance, M.; Brüschweiler, R. Toward Quantitative Interpretation of Methyl Side-Chain Dynamics from NMR by Molecular Dynamics Simulations. J. Am. Chem. Soc. 2007, 129, 14146−14147. (31) O’Brien, E. S.; Wand, A. J.; Sharp, K. A. On the Ability of Molecular Dynamics Force Fields to Recapitulate NMR Derived Protein Side Chain NMR Order Parameters. Protein Sci. 2016, 25, 1156−1160. (32) Millet, O.; Muhandiram, D. R.; Skrynnikov, N. R.; Kay, L. E. Deuterium Spin Probes of Side-Chain Dynamics in Proteins. 1. Measurement of Five Relaxation Rates per Deuteron in 13C-Labeled and Fractionally 2H-Enriched Proteins in Solution. J. Am. Chem. Soc. 2002, 124, 6439−6448. (33) Halle, B.; Wennerströ m, H. Interpretation of Magnetic Resonance Data from Water Nuclei in Heterogeneous Systems. J. Chem. Phys. 1981, 75, 1928−1943. (34) Wong, V.; Case, D. A. Evaluating Rotational Diffusion from Protein MD Simulations. J. Phys. Chem. B 2008, 112, 6013−6024. (35) Chatfield, D. C.; Szabo, A.; Brooks, B. R. Molecular Dynamics of Staphylococcal Nuclease: Comparison of Simulation with 15N and 13 C NMR Relaxation Data. J. Am. Chem. Soc. 1998, 120, 5301−5311. (36) Chatfield, D. C.; Wong, S. E. Methyl Motional Parameters in Crystalline L-Alanine: Molecular Dynamics Simulation and NMR. J. Phys. Chem. B 2000, 104, 11342−11348. (37) Xue, Y.; Pavlova, M. S.; Ryabov, Y. E.; Reif, B.; Skrynnikov, N. R. Methyl Rotation Barriers in Proteins from 2H Relaxation Data. Implications for Protein Structure. J. Am. Chem. Soc. 2007, 129, 6827− 6838. (38) Abascal, J. L. F.; Vega, C. A General Purpose Model for the Condensed Phases of Water: TIP4P/2005. J. Chem. Phys. 2005, 123, 234505. (39) Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2008, 120, 215−241. (40) Dunning, T. H., Jr. Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007−1023. (41) Kendall, R. A.; Dunning, T. H., Jr.; Harrison, R. J. Electron Affinities of the First-row Atoms Revisited. Systematic Basis Sets and Wave Functions. J. Chem. Phys. 1992, 96, 6796−6806. (42) Lovell, S. C.; Word, J. M.; Richardson, J. S.; Richardson, D. C. The Penultimate Rotamer Library. Proteins: Struct., Funct., Genet. 2000, 40, 389−408. (43) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High Performance Molecular Simulations through Multi-Level Parallelism from Laptops to Supercomputers. SoftwareX 2015, 1−2, 19−25.
(44) Wang, J.; Cieplak, P.; Kollman, P. A. How Well Does a Restrained Electrostatic Potential (RESP) Model Perform in Calculating Conformational Energies of Organic and Biological Molecules? J. Comput. Chem. 2000, 21, 1049−1074. (45) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; et al. A Point-Charge Force Field for Molecular Mechanics Simulations of Proteins Based on Condensed-Phase Quantum Mechanical Calculations. J. Comput. Chem. 2003, 24, 1999−2012. (46) Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E.; Simmerling, C. ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB. J. Chem. Theory Comput. 2015, 11, 3696−3713. (47) Vijay-Kumar, S.; Bugg, C. E.; Cook, W. J. Structure of Ubiquitin Refined at 1.8 Å Resolution. J. Mol. Biol. 1987, 194, 531−544. (48) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101. (49) Miyamoto, S.; Kollman, P. A. Settle: An Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models. J. Comput. Chem. 1992, 13, 952−962. (50) Hess, B. P-LINCS: A Parallel Linear Constraint Solver for Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 116−122. (51) Swope, W. C.; Andersen, H. C.; Berens, P. H.; Wilson, K. R. A Computer Simulation Method for the Calculation of Equilibrium Constants for the Formation of Physical Clusters of Molecules: Application to Small Water Clusters. J. Chem. Phys. 1982, 76, 637− 649. (52) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (53) Maragakis, P.; Lindorff-Larsen, K.; Eastwood, M. P.; Dror, R. O.; Klepeis, J. L.; Arkin, I. T.; Jensen, M. Ø.; Xu, H.; Trbovic, N.; Friesner, R. A.; et al. Microsecond Molecular Dynamics Simulation Shows Effect of Slow Loop Dynamics on Backbone Amide Order Parameters of Proteins. J. Phys. Chem. B 2008, 112, 6155−6158. (54) Lee, A. L.; Flynn, P. F.; Wand, A. J. Comparison of 2H and 13C NMR Relaxation Techniques for the Study of Protein Methyl Group Dynamics in Solution. J. Am. Chem. Soc. 1999, 121, 2891−2902. (55) Lienin, S. F.; Bremi, T.; Brutscher, B.; Brüschweiler, R.; Ernst, R. R. Anisotropic Intramolecular Backbone Dynamics of Ubiquitin Characterized by NMR Relaxation and MD Computer Simulation. J. Am. Chem. Soc. 1998, 120, 9870−9879. (56) Tjandra, N.; Feller, S. E.; Pastor, R. W.; Bax, A. Rotational Diffusion Anisotropy of Human Ubiquitin from 15N NMR Relaxation. J. Am. Chem. Soc. 1995, 117, 12562−12566. (57) Sheppard, D.; Li, D.-W.; Brüschweiler, R.; Tugarinov, V. Deuterium Spin Probes of Backbone Order in Proteins: 2H NMR Relaxation Study of Deuterated Carbon α Sites. J. Am. Chem. Soc. 2009, 131, 15853−15865. (58) Liao, X.; Long, D.; Li, D.-W.; Brüschweiler, R.; Tugarinov, V. Probing Side-Chain Dynamics in Proteins by the Measurement of Nine Deuterium Relaxation Rates per Methyl Group. J. Phys. Chem. B 2012, 116, 606−620. (59) Long, D.; Li, D.-W.; Walter, K. F. A.; Griesinger, C.; Brüschweiler, R. Toward a Predictive Understanding of Slow Methyl Group Dynamics in Proteins. Biophys. J. 2011, 101, 910−915. (60) Bowman, G. R. Accurately Modeling Nanosecond Protein Dynamics Requires at Least Microseconds of Simulation. J. Comput. Chem. 2016, 37, 558−566. (61) Yang, D.; Mittermaier, A.; Mok, Y.; Kay, L. A Study of Protein Side-Chain Dynamics from New 2H Auto-Correlation and 13C Crosscorrelation NMR Experiments: Application to the N-terminal SH3 Domain from Drk. J. Mol. Biol. 1998, 276, 939−954. (62) Best, R. B.; Clarke, J.; Karplus, M. The Origin of Protein Sidechain Order Parameter Distributions. J. Am. Chem. Soc. 2004, 126, 7734−7735. (63) Hu, H.; Hermans, J.; Lee, A. L. Relating Side-Chain Mobility in Proteins to Rotameric Transitions: Insights from Molecular Dynamics Simulations and NMR. J. Biomol. NMR 2005, 32, 151−62. 5047
DOI: 10.1021/acs.jpcb.8b02769 J. Phys. Chem. B 2018, 122, 5038−5048
Article
The Journal of Physical Chemistry B (64) Chou, J. J.; Case, D. A.; Bax, A. Insights into the Mobility of Methyl-Bearing Side Chains in Proteins from 3JCC and 3JCN couplings. J. Am. Chem. Soc. 2003, 125, 8959−8966. (65) Lindorff-Larsen, K.; Best, R. B.; Vendruscolo, M. Interpreting Dynamically-Averaged Scalar Couplings in Proteins. J. Biomol. NMR 2005, 32, 273−280. (66) Takemura, K.; Kitao, A. Water Model Tuning for Improved Reproduction of Rotational Diffusion and NMR Spectral Density. J. Phys. Chem. B 2012, 116, 6279−6287. (67) Chen, P.-c.; Hologne, M.; Walker, O. Computing the Rotational Diffusion of Biomolecules via Molecular Dynamics Simulation and Quaternion Orientations. J. Phys. Chem. B 2017, 121, 1812−1823. (68) Chen, P.-c.; Hologne, M.; Walker, O.; Hennig, J. Ab Initio Prediction of NMR Spin Relaxation Parameters from Molecular Dynamics Simulations. J. Chem. Theory Comput. 2018, 14, 1009−1019. (69) Anderson, J. S.; Hernández, G.; LeMaster, D. M. Prediction of Bond Vector Autocorrelation Functions from Larmor FrequencySelective Order Parameter Analysis of NMR Relaxation Data. J. Chem. Theory Comput. 2017, 13, 3276−3289. (70) Anderson, J. S.; LeMaster, D. M. Rotational Velocity Rescaling of Molecular Dynamics Trajectories for Direct Prediction of Protein NMR Relaxation. Biophys. Chem. 2012, 168−169, 28−39.
5048
DOI: 10.1021/acs.jpcb.8b02769 J. Phys. Chem. B 2018, 122, 5038−5048
