Accurate Methyl Group Dynamics in Protein Simulations with AMBER

Subscriber access provided by UNIV OF NEW ENGLAND ARMIDALE

B: Biophysical Chemistry and Biomolecules

Accurate Methyl Group Dynamics in Protein Simulations with AMBER Force Fields Falk Hoffmann, Frans A.A. Mulder, and Lars V. Schäfer J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b02769 • Publication Date (Web): 25 Apr 2018 Downloaded from http://pubs.acs.org on April 26, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Accurate Methyl Group Dynamics in Protein Simulations with AMBER Force Fields Falk Hoffmann,† Frans A. A. Mulder,‡ and Lars V. Schäfer∗,† †Theoretical Chemistry, Ruhr-University Bochum, Germany ‡Interdisciplinary Nanoscience Center (iNANO) and Department of Chemistry, University of Aarhus, Denmark E-mail: [email protected]

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

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 systemspecific 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 2 Saxis 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.

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 S 2

2


Page 2 of 36

Page 3 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 fields 16–22 have enabled nearly quantitative interpretation of protein backbone dynamics measured by NMR 2 relaxation in terms of squared generalized order parameters of amide NH groups, SN H , as

derived from the NMR data using the Lipari-Szabo (LS) formalism. 23,24 NMR relaxation data has not only been used to validate MD force fields a posteriori, but also as target data for optimizing the potential energy functions for 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 methyl groups in

13

13

C H2 2 H

C-, fractionally 2 H-labeled proteins report on reorientation motions of

C–2 H bond vectors on time scales faster than the overall tumbling of the molecule, i.e., from ps up to a few ns 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 S 2 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 best-characterized small

3



Page 4 of 36

proteins, including the availability of NMR data. In addition to generalized order parameters S 2 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 Methods, we calculate the time correlation functions (TCFs) of the methyl C–H bond vectors from the MD simulations, C(t) = hP2 [~µ(0) · µ ~ (t)]i ,

(1)

with the second order Legendre polynomial P2 [x] = (3x2 − 1)/2, the unit vector of the bond µ ~ (t), and h. . .i indicates 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, Z J(ω) =

∞

C(t) e−iωt dt .

(2)

0

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 ns. 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

4


Page 5 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


order parameter is defined as the long-time limit of Cint (t), S 2 = lim Cint (t) . t→∞

(3)

S 2 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 three-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 2 . Realistically capturing the correlation times of methyl group symmetry axis is coined Saxis

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/9 e−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 three-fold axis, assuming ideal tetrahedrality. Showalter and co-workers treated τCH3 as 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 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 field 18–20 (and related AMBER parameter sets) required adjusting the energy barriers of CH3 spinning in all methyl-containing amino acids except THR. This reparametrization was done based on CCSD(T) coupled cluster quantum chemical calculations of isolated dipeptides. Our MD simulation approach enables the prediction of relaxation rates and spectral densi-

5



ties, as well as of LS motional parameters S 2 and τf , without resorting to any information from NMR experiments. Our procedure is inspired by previous work of Bremi an co-workers 27 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 model 38 directly yield rotational correlation times that agree with experiments.

Methods Quantum chemical calculations To calculate the adiabatic energy barriers of methyl group rotation (spinning around its three-fold axis), relaxed potential energy surface (PES) scans were performed for isolated blocked dipeptides, ACE-X-NME, of the methyl-containing amino acids (X = ALA, MET, THR, LEU, ILE, VAL). The M06-2X density functional 39 was used in conjunction with the correlation-consistent double-ζ cc-pVDZ 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 threefold 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, were χ1 was set to 180◦ (trans) and +60◦ (gauche+), respectively. These values for χ1 were chosen, because they 6


Page 6 of 36

Page 7 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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*ILDN 18–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 set 44 and carried over to subsequent parameter sets, such as ff99SB 18 and ff99SB-ILDN, 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 non-bonded interactions. To increment the dihedral angle in 5◦ steps, as was done in the quantum chemical calculations described above, a stiff harmonic dihedral angle restraining potential with a force constant of 20000 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 Results and Discussion, 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 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 16235 TIP4P/2005 38 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 7



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 field 18–20 was used. To keep the temperature constant at 300 K, the velocity rescaling thermostat of Bussi and co-workers 48 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 SETTLE 49 and LINCS 50 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) method 52 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 energy 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

8


Page 8 of 36

Page 9 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of less than 0.2 K over the course of the 100 ns NVE simulations. 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 2 Saxis . In this approach, the internal TCFs of the Cmethyl –Hmethyl vectors (i.e., after removing

overall tumbling) were fitted with six exponentials plus an offset, 27

Cint (t) =

6 X

2 Ai e−t/τi + Slong ,

(4)

i=1

2 where Slong 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 P 2 2 fit to Eq. 4, we used 6i=1 Ai + Slong = 1, 0 ≤ Ai ≤ 1, 0 ≤ Slong ≤ 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 single-exponential function, 28

C(t) = Cint (t) · CO (t) = Cint (t) · e−t/τR ,

9


(5)


Page 10 of 36

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 a LS model, ef f

−t/τc,i 2 −t/τR,i 2 , C(t) = SN + (1 − SN He H) e

(6)

with effective correlation time τcef f = (τR τf )/(τR +τf ). The fit parameters are the correlation 2 times τR,i and τf as well as SN H . 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 website 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 between 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

J(ω) =

6 X i=1

2 τR Slong Ai τ i + , i 1 + (ω τef f )2 1 + (ω τR )2

(7)

with τief f = (τ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) model, 23,24 2 J(ω) = 5

1 2 S τ 9 axis R

1 + (ωτR )2

+

2 ) τcef f (1 − 91 Saxis

1 + (ωτcef f )2

,

(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: h i 2 2 1. Grid search within Saxis ∈ 0, 0.01, ..., 2Slong and τf ∈ 0 ps, 1 ps, ..., 2τfefitf . Here, we P 2 used Slong and τfefitf = 6i=1 Ai τi as starting points for the fitting parameters.

10


Page 11 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2 2. Direct (gradient-based) fit of Saxis and τf with the results from the previous grid search 2 ; τf ) of the fit. as starting points. This yielded the final parameter pair (Saxis

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 3 R(Dz ) = 16

qQe2 h ¯

2 [J(ωD ) + 4J(2ωD )] ,

2 1 qQe2 [9J(0) + 15J(ωD ) + 6J(2ωD )] , R(Dy ) = 32 h ¯ X Rmodel (i) − Ranalytical (i) 2 2 χ = . Rerror (i)

(9) (10) (11)

i∈{Dz ,Dy }

2 The fitting parameter pair (Saxis ; τ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 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 either 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 Supporting Information shows that the TCFs are correctly captured up to lag-times of 10 ns, demonstrating the validity of the applied assumptions. 11



Page 12 of 36

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.2 kJ/mol in the original force field, whereas for ILE-Cγ , the barrier is overestimated by 5.1 kJ/mol (Table 1 and Figure 1; the potential energy curves for all side-chains are 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 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). Table 1: Potential energy barriers (in kJ/mol) of methyl group spinning in blocked dipeptides. 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. ALA MET THR VAL Original FF 15.5 9.0 11.0 18.4/17.3 Reparametrized FF 14.2 7.2 11.0 13.1/12.1 CCSD(T) 14.2 7.1 11.4 14.0/11.5

12


LEU ILE 16.8/16.2 17.4/13.5 13.9/13.3 12.4/10.7 14.1/12.9 12.2/10.7

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


V (kJ/mol)

Page 13 of 36

16 ALA 14 12 10 8 6 CCSD(T)/cc-pVTZ original FF 4 reparametrized FF 2 M06-2X/cc-pVDZ 0 -180 -160 -140 -120 -100 -80 dihedral N-Cα-Cβ-H (degrees)

-60

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 M062X/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.

Table 2: Difference of the force constants kdih of the dihedral energy terms for methyl group spinning between the reparametrized and the original force field. 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-CHmethyl dihedral angle in MET, kdih = 1.39467 kJ/mol in the original force field. THR was not reparametrized. 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

13



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 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-cc-pVQZ quadruple-ζ basis sets. For ethane in vacuo, M06-2X/cc-pVDZ 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 methylcontaining 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 three-fold axis has the general form Vdih = kdih [1 + cos(3φ)]. Thus, each term contributes ca. 2 kdih to the methyl rotation barrier. 36 For example, for ALA, reducing kdih by 0.06964 kJ/mol (Table 2) should reduce the barrier by ca. 2 · 9 · 0.06964 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 non-bonded (electrostatic and van der Waals) interactions. However, because changes to the non-bonded interactions

14


Page 14 of 36

Page 15 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 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 2 , and correlation times of fast internal motions, τf , and then turn to order parameters, Saxis 2 deuterium relaxation rates and spectral densities. We calculated the LS parameters Saxis

and τf from our MD simulations with the spectral density mapping approach described in Methods, 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 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 (RMSD) are summarized in Table 3. The results obtained from our simulations in the NVE ensemble are very similar and consistent to 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). 2 Our results show that concerning Saxis , both force fields agree equally well with the NMR 2 experiment, with an RMSD in Saxis of about 0.1. This agreement is considered encouraging,

given that obtaining accurate methyl order parameters from MD simulations is considered 2 much harder than for backbone amides, for which RMSDs in SN H of about 0.05 between

MD and NMR are typically achieved. For ubiquitin, previous MD simulation studies re2 ported RMSDs from the NMR-derived 54 Saxis of 0.125 (Long and co-workers, obtained from

1 µs simulations with the AMBER ff99SBnmr1-ILDN force field 59 ), ca. 0.15 (Bowman, 60 15



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 co-workers 54 and Liao and co2 workers. 58 Interestingly, the RMSD in Saxis 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 2 values derived from Lee and co-workers being on sets seem to be systematic, with the Saxis

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 2 values vary substantially between different types warranted. For example, as typical Saxis

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 2 of this work is not on Saxis , 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 rotamer transitions around all χ 2 2 angles of a side-chain is expected to have a strong influence on Saxis , with low Saxis values

indicating population of more than a single rotamer well. 61–63 Insufficient sampling in our 2 MD simulations would thus be expected to lead to a systematic overestimation of Saxis , 2 which indeed seems to be the case for the low order parameters (Saxis < 0.4), which are all

16


Page 16 of 36

Page 17 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 2 (here: ten 100 ns simulations) and with a single molecule. The noticeable deviations in Saxis

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 2 for Cγ1 /Cγ2 of 0.59/0.68, respectively, which are considerably MD simulations predict Saxis

higher than the NMR deuterium relaxation values of 0.38/0.34. 58 Likewise, for ILE13-Cγ2 , 2 = 0.75, whereas NMR relaxation yields 0.52. VAL70 and ILE13-Cγ2 have MD predicts Saxis

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 side-chain 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 rotamer transitions in VAL70 might be underestimated in the MD simulations, and rotamer 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). 2 Averaging Cint (t) between 25–50 ns yields Saxis = 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 ILE13-Cγ2 and both VAL70 methyl groups were also overestimated in the MD simulations of Liao and co-workers. 58 2 Figure 2B shows that compared to Saxis , the improvements due to force field reparametriza-

tion 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

17



Page 18 of 36

A)

B)

2 (A) and corresponding correlation time τf (B) Figure 2: Methyl axis order parameter Saxis from MD simulations with the AMBER ff99SB*-ILDN force field without (red dots) and with adjusted methyl rotation barriers (blue triangles) 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) 2 and root mean square deviation (RMSD) of methyl order parameter Saxis between NMR and MD. Simulation RP Original FF (NPT) 0.91 Reparametrized FF (NPT) 0.85 Reparametrized FF (NVE) 0.85 18

RS 0.91 0.88 0.86


RMSD 0.10 0.13 0.13

Page 19 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 ps time scale while the typical tumbling time of a small protein is a few ns. In 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 side chains 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 ms, were systematically lower than the corresponding S 2 values derived from deuterium relaxation. 64,65 This indicates rotamer 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 S 2 in future work. In any case, it is thus beneficial not to restrict the comparison of MD and NMR only to S 2 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

19



A)

B)

Figure 3: Relaxation rates RM D (Dy ) (A) and RM D (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) 21.3 21.0 Reparametrized FF (NPT) 15.8 9.3 Reparametrized FF (NVE) 16.6 9.4

20


Page 20 of 36

Page 21 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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-ns) 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 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) 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), which was also used in the experimental NMR study. 58 Figure 4 compares these MD-based spectral densities 2 and to the experimental curves, which were constructed from the reported values for Saxis

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 compares the spectral densities from NMR experiments in 99.9% D2 O 58 with those from MD simulations in 100% H2 O. 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, D2 O has a higher viscosity than H2 O and hence τR is somewhat longer. In their

21



ALA28-CB

ILE61-CG2

LEU50-CD2

MET1-CE

VAL17-CG2

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 2 line, and the fitting parameters Saxis and τf are reported. The spectral densities from NMR, obtained with a global tumbling time of 5.0 ns in D2 O, 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 ns 57 from 15 N relaxation in 90% H2 O/10% D2 O. 22


Page 22 of 36

Page 23 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


NMR deuterium relaxation study, Liao and co-workers 58 used the ratio of the viscosities of D2 O and H2 O to obtain the global tumbling correlation time of 5.0 ns from the value of 4.16 ns measured in previous

15

N relaxation experiments in 90% H2 O/10% D2 O. 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 2 time of 4.16 ns in H2 O (assuming that Saxis 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 D2 O/H2 O viscosity ratio to mimic the slower tumbling in D2 O 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-ns) 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 NMRderived LS parameters 58 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

23



ALA28-CB

ILE61-CG2

LEU50-CD2

MET1-CE

VAL17-CG2

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, respec2 tively. The internal TCFs from NMR, calculated from Saxis and τf reported by Liao and 58 co-workers, are shown in green. 24


Page 24 of 36

Page 25 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


to reproduce site-specific spectral densities of backbone amide bonds in ubiquitin and protein G. 66 An alternative route was followed by Chen and co-workers, 67,68 who i) also focused exclusively on backbone amide bond motions (and not on side-chains), and ii) describe 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 can not be expected to yield realistic rotational diffusion times at different conditions, such as higher temperature or pressure. Recently, Anderson and co-workers 69 described an MD simulation approach to predict time correlation functions of backbone amide and non-methyl 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.

Conclusions In the present work, we describe an MD simulation approach to directly calculate deuterium relaxation rates, spectral densities, as well as the associated Lipari-Szabo dynamic param2 eters Saxis 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

25



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 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 sidechain 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.

Supporting Information Available Supporting Information. Tables with relaxation rates. Validation of factorization of total TCF into internal and overall TCFs and of their functional forms. Potential energy curves for 2 methyl rotation in all side-chains. Correlation plot of Saxis from NPT and NVE simulations.

Plot of the rotameric transitions of the VAL70 side-chain. Plots of spectral densities from MD and NMR. Plots of TCFs from MD and NMR. This information is available free of charge via the Internet at http://pubs.acs.org

Acknowledgement 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 26


Page 26 of 36

Page 27 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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. (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.

27



(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. Prot. 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 RelaxationBased 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 All-Atom 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.

28


Page 28 of 36

Page 29 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(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 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 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 Self-Consistent 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. Analysis of Experimental Results. J. Am. Chem. Soc. 1982, 104, 4559–4570. (25) Li, D.-W.; Brüschweiler, R. NMR-Based Protein Potentials. Angew. Chemie Int. Ed. 2010, 49, 6778–6780.

29



(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 TimeScale 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

13

C-Labeled and Fractionally 2 H-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.

30


Page 30 of 36

Page 31 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(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

15

N 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 2 H 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 Jr., T. H. 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 Jr., T. H.; 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 2000, 40, 389–408.

31



(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

32


Page 32 of 36

Page 33 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 2 H and

13

C 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

15

N 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: 2 H 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.

33



(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 2 H Auto-Correlation and 13 C Cross-correlation 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. (64) Chou, J. J.; Case, D. A.; Bax, A. Insights into the Mobility of Methyl-Bearing Side Chains in Proteins from 3 JCC and 3 JCN 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.

34


Page 34 of 36

Page 35 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(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 Frequency-Selective 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.

35



Graphical TOC Entry

Original force ﬁeld NMR

Reparametrized force ﬁeld

36


Page 36 of 36

Accurate Methyl Group Dynamics in Protein Simulations with AMBER

Recommend Documents