Rotational Dynamics of Proteins from Spin Relaxation Times and Molecular Dynamics Simulations.

Conformational fluctuations and rotational tumbling of proteins can be experimentally accessed with nuclear spin relaxation experiments. However, interpretation of molecular dynamics from the experimental data is often complicated, especially for molecules with anisotropic shape. Here, we apply classical molecular dynamics simulations to interpret the conformational fluctuations and rotational tumbling of proteins with arbitrarily anisotropic shape. The direct calculation of spin relaxation times from simulation data did not reproduce the experimental data. This was successfully corrected by scaling the overall rotational diffusion coefficients around the protein inertia axes with a constant factor. The achieved good agreement with experiments allowed the interpretation of the internal and overall dynamics of proteins with significantly anisotropic shape. The overall rotational diffusion was found to be Brownian, having only a short subdiffusive region below 0.12 ns. The presented methodology can be applied to interpret rotational dynamics and conformation fluctuations of proteins with arbitrary anisotropic shape. However, a water model with more realistic dynamical properties is probably required for intrinsically disordered proteins.

Introduction

Conformational fluctuations and the entropy of proteins play a significant role in their functionality and interactions with other biomolecules. Conformational fluctuations and the overall Brownian tumbling of proteins are experimentally accessible through the spin relaxation times of 15N and 13C nuclei measured with nuclear magnetic resonance (NMR) techniques.[1−7] The spin relaxation rates have been used to, for example, analyze conformational entropies,[1,8−11] binding entropies,[1,12] resolve sampled structures,[3−5,13] and validate molecular dynamics (MD) simulations.[14−21] These analyses are almost exclusively based on the separation of the internal conformational fluctuations and the overall rotational tumbling.[22,23] Also, the isotropic overall diffusion is often assumed, whereas analysis of anisotropic molecules is significantly more complicated.[1,2,24−27] Thus, new approaches are needed to interpret spin relaxation times measured from anisotropic or intrinsically disordered molecules.

Classical MD simulation methods are promising tools to interpret spin relaxation experiments for molecules with significantly anisotropic shape or correlations between internal and overall rotational motions. Practical applications are, however, limited by inaccuracies in the force field descriptions and the available time scales in the simulations.[17−19,28−31] The main issues have been the overestimated overall rotational diffusion of proteins due to inaccuracies in water models[19,29,32] and the insufficient accuracy of correlation functions calculated from single molecules in MD simulations.[30,33]

In this work, we overcome these issues by assuming that the overall rotational dynamics of protein follows anisotropic rigid body diffusion. Diffusion coefficients around inertia axes are directly calculated from angular displacements. The diffusion coefficients are then used to determine the contribution of the overall rotational tumbling to the rotational correlation functions of N–H bonds in the protein backbone. This reduces the required simulation length for the accurate determination of the rotational correlation functions. Furthermore, the overestimated overall Brownian tumbling rates due to the inaccurate water model can be corrected during the correlation function calculation by scaling the diffusion coefficients in all directions with a constant factor. The corrected correlation functions can be used to interpret the spin relaxation experiments for proteins with arbitrarily anisotropic shapes.

The developed approach is demonstrated by interpreting the experimental spin relaxation data of C-terminal domains of TonB proteins from Helicobacter pyroli (HpTonB-92)[34] and from Pseudomonas aeruginosa (PaTonB-96),[35] having 92 and 96 residues, respectively. Both proteins have significantly anisotropic shape, which would complicate the standard spin relaxation data analysis.[1,2,24−27]

Methods

Spin Relaxation Experiments and Rotational Dynamics of Molecules

Molecular dynamics of the protein backbone residues and spin relaxation experiments can be connected by using the spectral density J(ω)which is the Fourier transformation of the second-order rotational correlation function for N–H bond vectorwhere θ is the N–H bond angle between times t′ and t′ + t and angular brackets refer to the ensemble average. Connection to the experimentally measured spin relaxation times T1, T2 and the nuclear Overhauser effect (NOE) relaxation is given by the Redfield equations[36,37]where ωN and ωH are the Larmor angular frequencies of 15N and 1H, respectively, and the number of bound protons NH = 1 for N–H bonds. The dipolar coupling constant is given bywhere μ0 is the magnetic constant or vacuum permeability, ℏ is the reduced Planck constant, and γN and γH are the gyromagnetic constants of 15N and 1H, respectively. The average cubic length is calculated as ⟨rNH3⟩ = (0.101 nm), and the value of Δσ = −160 ppm is used for the chemical shift anisotropy of N–H bonds in proteins.[37,38]

Spin relaxation experiments are typically interpreted for proteins by assuming that the motions related to the overall Brownian tumbling and conformational fluctuations are independent.[39] The rotational correlation function for each N–H bond can be then written as[1,2,22,23,39]where CI(t) and CO(t) are correlation functions for the internal dynamics and overall rotations, respectively. Conformational fluctuations can be described in this approximation by using the square of the order parameter with respect to molecular axes S2, which are given by the plateau of the internal rotational correlation function. Timescales for the fluctuations can be characterized by using the effective correlation timewhere is the reduced correlation function.[23]

The overall rotational correlation function is often described by approximating the protein as a rigid body. For arbitrarily anisotropic molecules, the correlation functions can be presented as a sum of five exponentials[2,24]where the prefactors A depend on the directions of chemical bonds with respect to the molecular axes[24,26] and the time constants τ are related to the diffusion constants around three principal axes of a molecule (D , D , and D) through equations[2,24]where

The simplest approach to extract molecular dynamics from the experimental data is the original “model-free analysis”,[23] where an isotropic diffusion is assumed for the overall rotation of the protein. This reduces eq to a monoexponential form and the overall rotational dynamics can be described with a single time constant τc. Also, the internal correlation functions for each residue are assumed to decay exponentially with a single time constant τeff toward to the square of the order parameter S2. The three parameters (τc, τeff, and S2) can be then successfully resolved from a fit to the experimental data. However, the number of parameters to be fitted increases if the protein experiences an anisotropic overall diffusion or has several timescales for internal motions. In this case, the fitting becomes often ambiguous, even if the experimental data would be measured with multiple magnetic field strengths.[1,26,40] The anisotropic rotational diffusion is sometimes described with hydrodynamical calculations but they are sensitive to the estimation of the hydration shell around the protein.[41]

A rough estimate for the timescale of overall rotational dynamics is often given by using the T1/T2 ratio.[37] This is based on the assumptions that T1 and T2 are independent of the internal motions and that the overall dynamics is isotropic. The spectral density then reduces toand the correlation time describing the overall rotational motion, τc′, can be estimated by numerically minimizing the equationwith respect to the experimentally measured T1/T2 ratio.

Rotational Dynamics from MD Simulations

A classical MD simulation gives a trajectory for each atom in the system as a function of time. Rotational correlation functions for each bond can be then directly calculated from the trajectories by eq and used to calculate the spin relaxation times through eqs –5. The resulting values can be compared to experimental data to assess simulation model quality[14−21,31,42] and to interpret experiments.[21,42,43]

The direct comparison with experiments is, however, often complicated by the insufficient statistics for the calculated correlation functions and the overestimated rotational diffusion due to inaccuracies in the used water models.[29,30,32] Here, we show that the statistical accuracy of the contribution of the overall tumbling to the correlation functions, CO(t), in eq , can be increased for rigid proteins by directly calculating the diffusion coefficients of the inertia axes. The rotational diffusion coefficients can be related to the timescales τ of the correlation function for anisotropic rigid body rotation in eq by using the relations in eq .[24]

The rotational diffusion coefficients are calculated by fitting a linear slope to the square angle deviation of the inertia axes (see below). Lag times up to one hundredth of the total simulation length were used. This is expected to be the maximum lag time for the good statistics of rotational dynamics analyzed from a single molecule in MD simulations.[33] Error bars for the diffusion coefficients were defined to include results when the lag time was varied with ±1 ns. This requires less simulation data for the good statistics than a direct fit of the multiexponential sum in eq to the rotational correlation function calculated from MD simulation. In addition, the overestimated rotational diffusion due to the water model[19,29,32] can be corrected by scaling the diffusion coefficients around all inertia axes by a constant factor. This approach takes into account the anisotropic shape of the molecule. This is a significant advancement to the previous studies, which assume isotropic rotational diffusion with a single exponential rotational correlation function[10,15−17,44] or use order parameters to compare simulations with experimental data.[14,17,18,44]

The practical analysis can be divided into seven steps as follows:where ⟨Δα2⟩, ⟨Δβ2⟩, and ⟨Δγ2⟩ are the mean square angle deviations of the rotation around inertia axes from the longest protein inertia axis to the shortest, respectively.

The total rotational correlation functions C(t) for N–H bond vectors in a protein are directly calculated from the MD simulation trajectory by applying eq .

The rotational correlation functions for internal dynamics CI(t) are calculated from the MD simulation trajectory by removing the overall rotation of the protein.

The overall and internal motions are assumed to be independent and the overall rotational correlation function is calculated from eq as CO(t) = C(t)/CI(t).

The mean square angle deviations of rotation around protein inertia axes are calculated from the MD simulation trajectory.

rotational diffusion constants D, D, and D around inertia axes are calculated by fitting a straight line to the mean square angle deviations

contribution of the overall rotational tumbling to all correlation functions is assumed to follow eq with the timescales τ calculated from the rotational diffusion constants by using the relations in eq . Weighting factors A are determined by fitting the equation to the overall rotational correlation functions calculated from MD simulations in step 3.

The new correlation functions are calculated by substituting internal correlation functions, CI(t), from step 2 and anisotropic rigid body rotational correlation functions, CO(t), from step 6 to eq giving

These correlation functions are then used to calculate spin relaxation times from eqs –5. The incorrect overall rotational diffusion due to a water model can be corrected at this point by scaling the rotational diffusion coefficients, that is, timescales τ, with a constant factor before calculating new correlation functions from eq . Here, we determine the optimal scaling factors separately for each system. Scaling factors between 1 and 4 are explored with the spacing of 0.1 and the value giving the best agreement with the experimental spin relaxation data is selected to be the optimal scaling factor.

Simulation and Analysis Details

All simulations were performed using Gromacs 5[54] software and Amber ff99SB-ILDN[55] force field for proteins. The protein was solvated to tip3p,[56] tip4p,[56] or OPC4[57] water models. Initial structures were taken from the lowest-energy NMR structures of HpTonB-92 (PDB code: 5LW8)[34] and PaTonB-96 (PDB code: 6FIP).[35] The results from different initial conformations of both proteins with the tip3p water model are shown in Section S2 in the Supporting Information. The temperature was coupled to the desired value with the v-rescale thermostat,[58] and the pressure was isotropically set to 1 bar using a Parrinello–Rahman barostat.[59] Time step was 2 fs, Lennard-Jones interactions were cut off at 1.0 nm, particle mesh Ewald[60,61] was used for electrostatics, and LINCS was used to constrain all bond lengths.[62] The simulated systems are listed in Table with the references giving access to the trajectories and the related simulation files. Equilibration of the trajectories was followed by monitoring the protein root-mean-square-deviation, inertia tensor eigenvalues, and rotation angles. Sufficient amount of data was omitted from the beginning of simulation trajectories to remove the significant fluctuations in these parameters. If such fluctuations were not observed, the first 10 ns of the trajectory was omitted as an equilibration period.

Table 1

Simulated Systems and Rotational Diffusion Coefficients (rad2·107/s) Calculated from Simulations

protein	watera	T (K)b	ts (ns)c	ta (ns)d	Dx	Dy	Dz	D∥/D⊥e	Davf	filesg
PaTonB-96	tip3p	298	400	300	4.2 ± 0.1	4.4 ± 0.1	10.4 ± 0.1	2.42 ± 0.1	6.4 ± 0.1	(45)
PaTonB-96	tip4p	298	400	390	1.81 ± 0.01	2.06 ± 0.03	4.55 ± 0.03	2.35 ± 0.04	2.80 ± 0.02	(46)
PaTonB-96	tip4p	310	400	390	2.60 ± 0.02	2.22 ± 0.05	5.0 ± 0.1	2.07 ± 0.09	3.26 ± 0.07	(47)
PaTonB-96	OPC4	310	1200	1190	2.01 ± 0.01	2.19 ± 0.01	5.01 ± 0.03	2.39 ± 0.02	3.07 ± 0.01	(48)
HpTonB-92	tip3p	310	570	370	8.25 ± 0.05	7.67 ± 0.06	15.9 ± 0.3	1.99 ± 0.06	10.6 ± 0.2	(49)
HpTonB-92	tip3p	303	800	790	6.24 ± 0.02	7.04 ± 0.03	11.9 ± 0.2	1.80 ± 0.03	8.40 ± 0.07	(50)
HpTonB-92	tip4p	310	470	370	3.6 ± 0.1	3.24 ± 0.01	6.3 ± 0.3	1.8 ± 0.1	4.4 ± 0.2	(51)
HpTonB-92	tip4p	303	400	200	2.7 ± 0.1	2.71 ± 0.02	5.6 ± 0.5	2.1 ± 0.2	3.7 ± 0.2	(52)
HpTonB-92	OPC4	310	800	790	2.85 ± 0.01	2.70 ± 0.01	5.56 ± 0.01	2.00 ± 0.01	3.70 ± 0.01	(53)

Water model used in simulation.

Simulation temperature.

Total simulation time.

Analyzed simulation time.

Citation to a repository containing the simulation data.

The rotational correlation functions are calculated with gmx rotacf from Gromacs package.[63] The overall rotation was removed for CI(t) calculation by using a fit option of the gmx trjconv tool in Gromacs package.[63] The order parameters S2 were determined by averaging the rotational correlation functions from the oriented trajectory, CI(t), over the lag times above 50 ns. The effective correlation times were then calculated by eq . Inertia axes of proteins were calculated with the compute_inertia_tensor function from MDTraj python library.[64]

Spectral density was calculated by fitting a sum of 471 exponentials with timescales from 1 ps to 50 ns with logarithmic spacingto the new correlation function from eq by using the lsqnonneg routine in MATLAB.[65] The Fourier transform was then calculated by using the analytical function for the sum of exponentials

A similar approach has been previously used for the lamellar lipid and surfactant systems in combination with solid-state NMR experiments.[66,67] All computer programs used for the analysis are available from ref (68).

Spin Relaxation Experiments

NMR experiments were recorded on a Bruker Avance III HD NMR spectrometer operated at 1H frequency of 850.4 MHz equipped with a cryogenic probe head. The longitudinal (T1), transverse (T2), and 1H–15N-heteronuclear NOE spin relaxation times for the backbone 15N atoms of HpTonB-92[34] were collected at 303 K using the well-established NMR pulse sequences described previously.[37,69,70] The similarly detected spin relaxation data for PaTonB-96 at 298 K are also reported in another publication.[35] The T1 and T2 relaxation times were measured using the following series of the delays: 10, 50, 100, 200, 300, 500, 800, 1000, 1200, and 2000 ms for T1 and 16, 64, 96, 128, 156, 196, 224, and 256 ms for T2. Recycle delays of 3.0 and 2.0 s were used for T1 and T2 experiments, respectively. The relaxation rates (R1 = 1/T1, R2 = 1/T2) were calculated as an exponential fit of a single exponential decay to peak intensity values: I(t) = I0 exp(−t/T1) or I(t) = I0 exp(−t/T2), where I(t) is the peak volume at a time t. The 15N{1H}-NOE measurements were carried out with a recycling delay of 5.1 s with and without saturation of the amide protons. The 15N{1H}-NOE values were derived from the volumes of the heteronuclear single-quantum coherence (HSQC) peaks using the equation of ν = I/I0. The relaxation data were processed and analyzed using Bruker Dynamic Center software (version 2.1.8).

Results and Discussion

Global Rotational Dynamics of the Protein

The mean square angle deviations for the rotation of the PaTonB-96 protein around inertia axes in the simulation with the OPC4 water model are shown in Figure . This is the longest simulation data set in this work (1.2 μs), and the linear behavior of the mean square angle deviations is observed for the lag times up to one hundredth of the total simulation length (12 ns), which is expected to be the maximum lag time for the good statistics of rotational dynamics analyzed from a single molecule in MD simulations.[33] Deviations from the linear behavior are only seen with the lag times longer than this limit, as also demonstrated for the shorter simulations with tip4p water at two different temperatures in Figures S1 and S2 in the Supporting Information. The plots with log–log scale in Figures , S1, and S2 reveal a weakly subdiffusive region only below very short timescales of approximately 0.12 ns. Thus, we conclude that the protein experiences the Brownian rotational tumbling with a good approximation. The diffusion coefficients can be then calculated from the slope of the mean square angle deviations according to eq by using the lag times less than one hundredth of the total MD simulation length. The error bars were calculated by varying the lag time with 1 ns to both directions. The data from HpTonB-92 protein (not shown) led to similar conclusions.

Figure 1

Mean square angle deviations of the rotation around inertia tensor axes calculated from PaTonB-96 simulation with the OPC water model. The data are shown with linear (top) and logarithmic scale (bottom).

The resulting rotational diffusion constants from different simulations are summarized in Table . As expected, the rotational diffusion coefficients increase with the temperature and the decreasing size of a protein. The values are, however, larger than expected from the experimental T1/T2 ratio analyzed with eq and from the previously reported values for proteins with similar sizes,[71] especially when tip3p water model is used. Similar results were previously explained by the overestimated water self-diffusion of the tip3p water model.[19,29,32]

The analysis leading to the new correlation functions in eq (see Methods section) is exemplified in Figure for three residues located in different domains of PaTonB-96 with different characteristic rotational dynamics. The flexible C-terminus is represented by the residue 341, more rigid β-sheet by the residue 331, and a flexible loop between two β-strands by residue 322 (see the labeling in Figure ). The total correlation functions C(t) of all residues in Figure (top, solid lines) decay toward zero within ∼10–50 ns. The internal correlation functions CI(t) in Figure (middle) decay to a plateau value, which defines the square of the order parameter S2. As expected, the internal correlation function for residue 331 in the rigid β-sheet rapidly decays to the largest order parameter value, whereas the correlation functions of the residues in the loop and C-terminus decay slower to the smaller order parameter values because of the larger conformational ensemble sampled by these regions.

Figure 2

Rotational correlation functions calculated from MD simulations of PaTonB-96 with the tip4p water model at 298 K for residues at different regions. (Top) Total correlation functions C(t) calculated from MD simulation (solid lines) and new correlation functions determined from eqs and 8 by using rotational diffusion constants and fitted prefactors (dashed lines); (middle) correlation functions for internal motions calculated from simulation with removed overall protein rotation; and (bottom) correlation function for overall motions determined as CO(t) = C(t)/CI(t) (solid lines) and by fitting to eq with timescales from rotational diffusion coefficients in Table (dashed lines).

Figure 6

(A) Structures sampled by PaTonB-96 from MD simulations with tip4p at 298 K (100 structures from 400 ns long trajectory). Secondary structures are color-labeled with Visual Molecular dynamics;[73,74] α-helixes and β-strands are red and blue, respectively. Residues 246–251, 320–326, and 338–342 with increased internal dynamics are yellow and α-helix fluctuations between two orientations (residues 266–270) are violet in the left column. Terminal ends are labeled with N and C. The structure from left is rotated with approximately 100° to the figure on right. (B) Spin relaxation times from experiments (circles) and tip4p simulations (squares) with rotational diffusion coefficients divided by a constant factor of 1.2 at 298 K. Order parameters and effective internal correlation times calculated from simulations.

The overall rotational correlation functions, CO(t) = C(t)/CI(t), are shown in Figure (bottom, solid lines). Also, the correlation functions of anisotropic rigid body rotation from eq are shown in Figure (bottom, dashed lines). The timescales for the latter, τ, are given by the rotational diffusion coefficients from the simulation and the relations in eq . The prefactors, A, are determined by fitting eq to the overall rotation correlation functions, CO(t), calculated from the MD simulation. The new correlation functions, determined from eq and shown in Figure (top, dashed lines), are indistinguishable from the correlation functions calculated from the original MD simulations with the lag times shorter than one hundredth of the total simulation time (approximately 4–12 ns), which is the maximum lag time for the good statistics in single-molecule MD simulations.[33] This suggests that the anisotropic rigid body diffusion model (eq ) and the separation of internal and global motions (eq ) are good approximations for the proteins studied in this work. The analytical description of the overall rotation with eq in the new correlation functions clearly reduces the statistical fluctuations with the long lag times in Figure . The effect is most visible for the flexible C-terminus (residue 341) having the smallest, thus the least detectable, contribution from the overall rotation of the protein due to the small order parameters.

Global Rotational Dynamics in Simulations and Experiments

Spin relaxation times of HpTonB-92 are compared between the experiments and simulations using two different water models in Figure . The simulation with tip3p water model underestimates the T1/T2 ratios, suggesting too fast overall rotational diffusion dynamics.[72] This is in agreement with the previous study, where the overestimated rotational diffusion was attributed to the self-diffusion of tip3p.[19,29,32] On the other hand, simulation results with tip4p water model show better agreement with the experimental data in Figure .

Figure 3

15N spin relaxation times for HpTonB-92 from experimental data (circles) and MD simulations with different water models (squares).

To see if the discrepancy in spin relaxation times for simulations with tip3p water model could be explained by the overestimated overall diffusion of the protein, the diffusion coefficients were divided by the optimal scaling factor before applying eq to calculate the new correlation functions. The scaling factor value of 2.9 gave the best agreement with the experimental spin relaxation data. The spin relaxation times calculated from the new correlation functions after scaling the rotational diffusion coefficients with the optimal scaling factor value are shown in Figure .

Figure 4

(A) Structures of HpTonB-92 from the MD simulations with tip3p at 303 K (100 structures taken from 400 ns long trajectory). Secondary structures are color-labeled with Visual Molecular dynamics;[73,74] α-helices are highlighted in red and β-strands in blue. Terminal ends are labeled with N and C. The structure from left is rotated with approximately 150° to the figure on right. (B) Spin relaxation times from experiments (circles) and tip3p simulations (squares) with rotational diffusion coefficients divided by a constant factor of 2.9 at 303 K. Order parameters and effective internal correlation times calculated from simulations.

Similar comparison for the spin relaxation times of PaTonB-96 between experiments and simulations with tip3p, tip4p, and OPC4 water models is shown in Figure . The experimentation of the OPC4 water model was inspired by the recent study reporting significant improvements in lipid monolayer simulations when this water model was used.[75] The underestimation of T1/T2 ratio was also observed in the simulations of PaTonB-96 with tip4p and OPC4 water models when compared with the experiments. The discrepancy is, however, less severe than with tip3p, suggesting that the required scaling factor for the overall rotational diffusion should be smaller for tip4p and OPC4 water models. Indeed, the spin relaxation times calculated from PaTonB-96 simulation with the tip4p water model were found to be in good agreement with the experiments in Figure when the diffusion coefficients were divided with a constant factor of 1.2, which is smaller than 2.9 used for the tip3p simulation of HpTonB-92 above. The scaling factors used to correct the overall rotational diffusion of different proteins with different water models are shown in Table S1 together with the corresponding coefficients for self-diffusion of water.[29,57] Notably, the effect of 12 °C temperature difference on the spin relaxation times from tip4p simulations in Figure is significantly smaller than the observed differences between simulations and experiments or the changes due to the scaling of the diffusion coefficient.

Figure 5

Plots of experimental (circles) and simulated (squares) spin relaxation times for PaTonB-96.

The scaling of the overall rotational diffusion coefficients with a constant factor led to a good agreement with the experimental spin relaxation data for both systems simulated with different water models, as seen in Figures and 6. The good agreement with experiments suggests that the scaled rotational diffusion coefficients from MD simulations can be considered as an interpretation of the anisotropic rotational motion in NMR experiments. The scaled rotational diffusion coefficients from the simulations giving the best agreement with the experimental data are summarized in Table . In contrast to the unscaled diffusion constants in Table , these results are in line with the previously reported values for proteins with similar sizes.[71] Also, the timescales, τc′, estimated from eq are close to the average diffusion coefficient, τc = (6Dav)−1, in Table .

Table 2

Rotational Diffusion Coefficients (rad2·107/s) Giving the Best Agreement with Experimental Spin Relaxation dataa

	HpTonB-92	PaTonB-96
Dx	2.15 ± 0.01	1.51 ± 0.01
Dy	2.43 ± 0.01	1.72 ± 0.03
Dz	4.10 ± 0.01	3.79 ± 0.03
Dav	2.90 ± 0.03	2.30 ± 0.02
τc (ns)b	5.7 ± 0.1	7.2 ± 0.1
τc′ (ns)c	5.8 ± 0.1	6.9 ± 0.1

For HpTonB-92 construct, the values calculated from simulation with tip3p were scaled with 2.9 (spin relaxation data in Figure ), and for PaTonB-96, the values from tip4p simulation at 298 K were scaled by 1.2 (spin relaxation data in Figure ).

τc = (6Dav)−1.

Average overall residues given by eq .

Interpretation of Protein Internal Relaxation from MD Simulations

The good agreement of the spin relaxation times between the simulations with the scaled overall rotational diffusion coefficients and the experiments (Figures and 6) suggests that the simulations can be used to interpret the internal mobility of proteins from the experimental data.

Only small variations between different residues are observed for spin relaxation times of HpTonB-92 in Figure . This indicates a rather rigid protein structure, which is also seen in the MD simulation snapshots overlayed in Figure A. Only few residues in the terminal ends show slightly enhanced conformational fluctuations in the MD simulation and in spin relaxation data. In addition, some deviations from the average spin relaxation times are observed in the experimental data close to residues 210–222. Simulations of HpTonB-92 do not offer any explanation for this observation; however, the similar region in PaTonB-96 simulation shows fluctuations between two orientations of α-helix.[35] Exceptionally low order parameters and long effective correlation times are observed in simulations for residues 245–250 of HpTonB-92. Moreover, short T1 times are experimentally observed close to this region, but the interpretation is not straightforward as the low T1 times are not reproduced by MD simulations.

PaTonB-96 exhibits more internal mobility and the segments with enhanced conformational fluctuations are labeled with yellow color in Figure . The larger number of sampled conformations in both terminal ends is characterized by the low order parameters and long effective internal correlation times observed in the simulations. Enhanced conformational fluctuations are also observed for residues 320–326, which correspond to the loop between two β-strands. MD simulations predict low order parameters and long internal effective correlation times also close to residues 266–271, which can be explained by two different orientations sampled by the α-helix in this region (color-labeled with violet in Figure A). The orientational fluctuations of the similar short helix could also explain the above mentioned deviations of spin relaxation times for residues 210–222 of HpTonB-92.[34]

MD simulations can be used to analyze different components contributing to the rotational dynamics of individual N–H bonds. In this work, we have fitted a sum of 471 different timescales to the correlation functions according to eq . Most of the prefactors (α in eq ) are zero in all correlation functions after the fitting; thus, the timescales τ corresponding to nonzero prefactors are considered as the components contributing to the total relaxation process of each N–H bond. The prefactors are shown in Figure for the same residues of PaTonB-96, which were used to exemplify the correlation functions in Figure . As expected for residue 322 in the rigid β-sheet with large order parameter, the rotational relaxation is dominated by timescales of ∼5.5 and ∼8 ns, matching with the protein overall rotation. Also, the dynamics of residue 322 in the flexible loop of PaTonB is dominated by the timescales around ∼8 ns corresponding to the protein overall rotation; however, the fast motions from internal mobility are more evident than for the rigid β-sheet residue. This is in agreement with smaller order parameter observed in the flexible loop residues. On the other hand, the rotational dynamics of residue 341 in the flexible N-terminus of PaTonB is dominated by timescales below 3 ns, most likely related to the internal motion of the protein. Contributions from timescales around ∼13 ns to the dynamics of residue 341 probably arise from the slow conformational fluctuations of the N-terminus, rather than the overall rotational dynamics. This supports the conclusion that the large amount of sampled conformations lead to the small order parameters and large effective correlation times observed in Figure . Although the separation of rotational dynamics of individual N–H bonds to different components gives intuitively understandable results, it should be kept in mind that it is based on the fitting of a multiexponential sum to the simulation data and the solution of such fit is not unique.

Figure 7

Prefactors α corresponding to different timescales τ resulting from a fit of eq to correlation functions from MD simulation of PaTonB-96 at 298 K. The used correlation functions give a good agreement with experimental spin relaxation times as shown in Figure .

Conclusions

The experimental spin relaxation data for protein backbone N–H bonds were successfully reproduced by using the classical MD simulations for two different small domains. Thus, the simulation trajectories give an atomic resolution interpretation for protein dynamics measured with NMR experiments. Interpretation of the overall and internal dynamics was demonstrated for two proteins with anisotropic molecular shape and some flexible regions. Interpretation of the 15N spin relaxation data measured from such proteins has been very challenging with the previously available methods.[26,69]

The overall rotation of the studied proteins was found to be Brownian, having only a small subdiffusive behavior with short timescales below ∼0.12 ns, which could be contrasted with crowded environments, where anomalous diffusion is expected to be more significant.[76] The direct analysis of classical MD trajectories did not, however, reproduce the experimental 15N spin relaxation data. Comparison between the rotational diffusion coefficients and spin relaxation times between simulations and experiments suggested that the overall Brownian tumbling of proteins is too rapid in the simulations, in agreement with the previous report suggesting that the discrepancy arises from the inaccuracies in water models.[19,29,32] Scaling down the anisotropic diffusion coefficients in the simulation data led to a good agreement with the experimental data. Overall rotational diffusion coefficients were overestimated by a factor of ∼3 in the HpTonB-92 simulations with the tip3p water model, in line with previous studies.[28−30,77] Simulations with tip4p and OPC4 water models gave the spin relaxation times in reasonable agreement with experiments with scaling factors of ∼1–1.2, which are significantly less than that for tip3p. The scaling factors for different proteins with different water models are summarized in Table S1.

The similarity between the correlation functions from the original MD trajectory and the new correlation functions from eq suggests that the usage of the inertia axes and the separation of internal and the overall rotational motions (eq ) are good approximations for the above investigated proteins. This is in line with the previous studies of other proteins with well-defined structure.[10,29] However, it remains to be seen how well this and other related approaches[28,30,78] will succeed for intrinsically disordered proteins without the well-defined shape. Because the correction of the incorrect overall rotational diffusion due to the water model may become highly complicated for such proteins, it may be necessary to employ a water model giving correct overall rotational diffusion coefficients for biomolecules.[19,32]

As further demonstrated in ref (35), the approach presented in this work can be used to interpret the rotational dynamics of proteins with anisotropic shape from 15N spin relaxation data measured only with one magnetic field strength. This is a significant advancement over currently available methods, which may not be applicable in such cases, even though experimental data would be measured with multiple magnetic field strengths.