Motions and entropies in proteins as seen in NMR relaxation experiments and molecular dynamics simulations.

Subscriber access provided by Iowa State University | Library

Article

Motions and Entropies in Proteins as Seen in NMR Relaxation Experiments and Molecular Dynamics Simulations Olof Allnér, Nicolas Foloppe, and Lennart Nilsson J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/jp506609g • Publication Date (Web): 13 Oct 2014 Downloaded from http://pubs.acs.org on October 18, 2014

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 free 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 accessible to all readers and 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.

The Journal of Physical Chemistry B 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 47

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

Motions and Entropies in Proteins as Seen in NMR Relaxation Experiments and Molecular Dynamics Simulations

Olof Allnér, Nicolas Foloppe1 and Lennart Nilsson*

Karolinska Institutet Department of Biosciences and Nutrition Center for Biosciences SE-141 83 HUDDINGE Sweden

1

Present address: 51 Natal Road, Cambridge CB1 3NY, UK

*

Corresponding author

Telephone: +46-8-524 81099 E-mail: [email protected] Telefax: N/A

1 ACS Paragon Plus Environment


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 Molecular dynamics simulations of E. coli glutaredoxin1 in water have been performed to relate the dynamical parameters and entropy obtained in NMR relaxation experiments with results extracted from simulated trajectory data. NMR relaxation is the most widely used experimental method to obtain data on dynamics of proteins but it is limited to relatively short timescales and to motions of backbone amides, or in some cases 13C-H vectors. By relating the experimental data to the all-atom picture obtained in molecular dynamics simulations, valuable insights on the interpretation of the experiment can be gained. We have estimated the internal dynamics and their time scales by calculating the generalized order parameters (O) for different time windows. We then calculate the quasiharmonic entropy (S) and compare it to the entropy calculated from the NMR-derived generalized order parameter of the amide vectors. Special emphasis is put on characterizing dynamics that are not expressed through the motions of the amide group. The NMR and MD methods suffer from complementary limitations, with NMR being restricted to local vectors and dynamics on a timescale determined by the rotational diffusion of the solute, while in simulations it may be difficult to obtain sufficient sampling to ensure convergence of the results. We also evaluate the amount of sampling obtained with molecular dynamics simulations and how it is affected by the length of individual simulations, by clustering of the sampled conformations. We find that two structural turns act as hinges, allowing the alpha-helix between them to undergo large, long time scale motions that cannot be detected in the time window of the NMR dipolar relaxation experiments. We also show that the entropy obtained from the amide vector does not account for correlated motions of adjacent residues. Finally we show that the sampling in a total of 100 ns molecular dynamics simulation can be increased by around 50% by dividing the trajectory into 10 replicas with different starting velocities.

Keywords: quasiharmonic, diffusion-in-a-cone, correlation


Page 2 of 47

Page 3 of 47

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


Introduction Proteins are the most versatile class of biological macromolecules, performing a wide variety of tasks essential to the cell. Already in the early 1970’s, when only a handful of protein 3D structures were known, it was becoming evident that protein motions, over a range of length- and time-scales, are crucial for proper protein function.1-3 Several aspects of protein function depend on conformational changes, from induced fit and allosteric effects, to performing mechanical work.4 Conformational change is crucial in the many molecular machines5 that have been characterized in great detail in the last decade, such as the ribosome, RNA and DNA polymerases, or ATP synthase, and it may also be important for the evolution of new functions and folds.6 Given that proteins can, and do, change conformation it is of great interest to characterize the relative stabilities of the functionally relevant (macroscopic) states. While the overall thermodynamic parameters characterizing a chemical or conformational equilibrium (the Gibbs free energy difference ∆G=∆H-T∆S, where H, S and T are, respectively, the enthalpy, the entropy and the absolute temperature of the system) are in general possible to measure experimentally, it is not as easy to pinpoint the molecular contributions to the enthalpic (∆H) and entropic (T∆S) terms. The entropy is particularly elusive, and commonly used approaches based on NMR measurements7 of generalized order parameters, although useful, have limitations in that they only provide partial information from the motions of specific bonds (typically backbone amides, or in some cases

13

15

N-H

C-H vectors) ; furthermore correlations between different

vectors are not accounted for. From a thermodynamic perspective, the number of micro states available to a protein within a functional macro state is important even for processes that may not involve large structural changes (e.g. protein stability differences due to point mutations, different ligand affinities),8 since this determines the entropy of the macro state. The entropy can be decomposed in rotational, translational, and configurational (or vibrational) components. Additionally, entropy changes due to solvent and ion release from macromolecular solutes, related to the hydrophobic effect, may have to be considered.



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 free energy can be estimated using a variety of computational methods based on population samples obtained from molecular dynamics (MD) or Monte Carlo simulations.9-12 The accuracy of such estimates depends on the accuracy of the Hamiltonian used to describe the system, and on how well relevant regions of phase-space have been sampled in the simulation, but they do offer the possibility of dissecting the contributions to the free energy in atomic detail. Accurate (and fast) computational free energy estimates would have a range of applications in molecular design, for example in drug discovery.13 In addition to these practical concerns, the entropy and energy terms are of interest in their own right for a deeper understanding of the fundamental properties of biomolecules. Rigorous free energy calculation methods can however be computationally expensive, and more approximate methods have been developed, such as the Linear Interaction Energy approach,14-15 or methods based directly on calculation of the internal energy, complemented with approximations of the solvation free energy using PoissonBoltzmann (MM/PBSA) or generalized Born expressions (MM/GBSA).13, 16-17 The MM/PBSA and MM/GBSA methods need an entropy estimate in order to evaluate the free energy, and there are indications18-19 that these methods require significantly more sampling than was initially hoped for. ∆H could be calculated as the average of the Hamiltonian over the sampled conformations (strictly this is ∆E, but for most proteins the p∆V term is negligible at normal pressures due to small magnitude of ∆V20), and even if this average may have large error bars, it is directly available from a classical or quantum mechanical description of the system. For the entropy the situation is different, since this in principle involves summation over all conformations of the system that are accessible at a given temperature. The translational (Strans) and rotational (Srot) components of the entropy can be calculated using expressions for an ideal poly-atomic gas, which depend on the volume in which the molecule is confined and the mass and moments of inertia of the molecule. ∆Strans and ∆Srot are often neglected if the overall shape and composition of the system is not changing much. Svib, the vibrational entropy, on the other hand depends on internal motions and thus on the atomic interactions in the system, which may change considerably and to a varying degree over the system even if the overall shape remains relatively unperturbed. For a large system such as a protein there are only approximate methods available for the calculation of Svib.


Page 4 of 47

Page 5 of 47

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 most commonly used computational methods to estimate the configurational entropy are based on either a normal mode analysis of a single structure, or a quasi-harmonic analysis of the structural fluctuations in a simulation. In the quasi-harmonic approximation (QHA)21-23 the vibrational frequencies obtained from the eigenvalues of the covariance matrix of positional fluctuations of the N atoms in the system are entered into the expression for the entropy of an N-6 dimensional quantum-mechanical harmonic oscillator.24 In contrast to a standard normal-mode analysis, which is typically applied to the macromolecule in vacuum at 0 K after extensive energy minimization, the quasi-harmonic analysis incorporates effects of the solvent and of an-harmonic motions occurring at T > 0 K. If the conformational distribution sampled by the simulation is multimodal, the QHA overestimates the entropy since it treats the computed coordinate covariance as describing a single wide well in the probability distribution.25-27 The presence of higher order correlations, which are not accounted for in the QHA, also leads to overestimation of the entropy, in particular when the analysis is performed using cartesian coordinates.25, have been suggested recently,27,

29-31

28

Various improved estimators

and show potential to deliver significantly better

configurational entropies at a reasonable increase in computational complexity and cost. While the quasi-harmonic entropy is no “gold-standard” for computing configurational entropies it is commonly used. In the QHA framework it is straightforward to compute the entropy for any subsystem of interest, by simply using the covariance matrix of the atoms in the subsystem, and, in contrast to alternative distribution-oriented methods, it also provides quasi-harmonic normal mode vectors, which are useful in interpreting the underlying motions in the system. Thus the QHA was the method that was chosen for this study. The total free energy, as well as its enthalpic and entropic components, is experimentally measurable using a number of methods, but local contributions in general have to be assessed using indirect approaches, such as analysis of the effect of modifications of the system. NMR provides a view of the vibrational entropy of individual bond-vectors in spin relaxation and residual dipolar coupling experiments. This is achieved through the extraction of generalized order parameters, O,32-33 characterizing the extent of motion of the relevant vector. O2 varies between 0, for completely unrestricted motions, and 1, for a completely rigid system. Hence, with some simple assumptions concerning the nature of the motion described by the vector, for 5 ACS Paragon Plus Environment


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

Page 6 of 47

example that it can be described as freely diffusing within a limiting cone in space (“diffusion-ina-cone” model),34 O can be related to the entropy of the motion.7,

34-38

Several variants and

alternatives to the model free analysis of experimental NMR data have been proposed.39-42 They have also been evaluated with MD simulations using order parameters obtained either directly from the correlation functions calculated from the trajectory, or by applying an NMR-analysis protocol to the spectral density back-calculated from the trajectory.43-45 Generalized order parameters 32-33 for proteins have been extracted from ~100ps MD trajectories using the model free formalism.

46-50

For short alkanes the effects of coupling between slow

internal motions and fast (in this case) overall tumbling on the corresponding correlation functions obtained from 100ns MD simulations have been discussed.51 More recently 15N-H and 13

C-H order parameters calculated from simulations in the range 10 ns -100 ns have been used to

test force fields52-53 and to analyze the underlying motions in proteins 54-55 and nucleic acids.56-58 A 1.2µs simulation59 of ubiquitin essentially confirms that the shorter MD simulations yield appropriate results for order parameters from NMR relaxation experiments, which are limited in their time-resolution by the 5 ns-10 ns rotational correlation time of the macromolecule. Determination of generalized order parameters using the model-free approach

32-33

is possible

even if internal motions and overall rotational diffusion occur on similar time-scales, provided that the slow internal motions are due to infrequent jumps between states60. Such jump-like motions, and their influence on order parameters calculated from short simulations have been previously described47 and in more recent and longer simulations shown to be largely manageable when the analysis is effectively performed over time-windows that are shorter than the rotational tumbling of the molecule45, 53, 59, 61-63 Adequate sampling, yielding a statistically converged population, is crucial in many problems where simulation approaches are used. Different properties may require very different amounts of sampling,64 and convergence may be difficult to reach in practice,19, transitions.

67

65-66

due in part to rare

The brute force approach of extending the simulation length is not always practical,

and several enhanced-sampling methods have therefore been devised.68 Due to the chaotic character of molecular dynamics simulations69-71 different starting conditions (e.g. velocities, coordinates, system size)72 can be used to initiate a set of statistically independent simulations, 6 ACS Paragon Plus Environment

Page 7 of 47

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


which for a given total computational expense, gives better sampling than one single long trajectory.18,

63, 69, 71, 73-76

Multiple simulations can also help probe different pH regimes by

varying the protonation state of relevant titratable groups, providing an enhanced sampling of states that would be inaccessible at neutral pH, thus leading to an improved characterization of the system.77 In this work we have analyzed molecular dynamics (MD) simulations on reduced E. coli Glutaredoxin1 (GRX1), started from different NMR conformers.78 Glutaredoxins are ubiquitous small globular redox enzymes, typically formed of around 80 amino-acids, which adopt the socalled thioredoxin fold (a central β-sheet between flanking α-helices: βαβ-loop-ββα).79 Glutaredoxins are thiol-disulfide oxidoreductases which contribute to the redox regulation of cell homeostasis and are involved in redox signaling. Thus glutaredoxin active-sites can exist in reduced, oxidized or mixed-disulfide states. The present study used the reduced form of GRX. The dynamics of GRX1 has been studied extensively with NMR80 making GRX1 well suited as a model system for using simulations to evaluate the robustness of estimating NMR dipolar relaxation parameters and conformational sampling. We begin by estimating the global tumbling of the protein to find out the timescales available to study the internal dynamics in an NMR relaxation experiment. We then to investigate the internal dynamics by comparing the generalized order parameters (O) of the

15

N-1H amide spin

interaction vector calculated using different time windows from simulated trajectories, and obtained from NMR relaxation experiments. The order parameters are also back-calculated from the simulated NMR relaxation parameters to assess the robustness of the model free approach. The internal dynamics are further investigated by comparing the entropy estimated from an NMR-experiment using the diffusion in a cone model with the more extensive quasiharmonic method. We put special emphasis on the correlation of motions that the diffusion-in-a-cone model is unable to represent. After this we evaluate how well our simulations have sampled the motions of GRX1 with a special focus on how the division of the total simulation time into one or several trajectories affects the amount of sampled conformations. Finally, we provide a characterization of the



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

motions that are responsible for the large variations in the internal dynamics parameters observed in the different secondary structure elements of GRX1.


Page 8 of 47

Page 9 of 47

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


Methods NMR Relaxation Theory Correlation Functions Since general NMR relaxation theory has been extensively described in existing literature,81 the focus here will be on the equations used in our analysis and particularly on the relaxation of the 15

N-1H amide spin interaction vector. For relaxation to be effective, the spins must oscillate at

frequencies close to the Larmor frequency. The spectral density function, J (ω) , which describes the intensity of oscillations as function of angular frequency, ω, plays a central role in relaxation theory and is given by: ∞

J (ω ) = 2∫ C ( t ) cos (ωt ) dt 0

(1)

where C ( t ) is the correlation function of the spin interaction vector. If the overall and the internal motions of a protein are assumed to be statistically independent as in the model-free approach of Lipari-Szabo,33 the total correlation function, C (t) , can be separated into an overall, CO (t ) , part and an internal, C I (t ) , part: C (t ) = CO (t )CI (t )

(2)

This decorrelation assumption holds if the timescale of global motions is always longer than the timescales of the internal motions of any part of the molecule. The validity of this assumption has been investigated on several occasions41, 56, 60-61, 82 and has been found to generally hold at normal temperatures and pressures. For isotropic overall tumbling, which has been found to be valid for globular proteins,61 Co(t) is given by: CO (t ) = e − t /τ c

(3)

where τ c , the global correlation time, is a measure of the tumbling frequency of the protein in solution.



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

Page 10 of 47

The internal correlation function can be expressed as:

CI (t ) =

P2 ( cos χτ ,τ + t ) r 3 (τ ) r 3 (τ + t )

(4) τ

3x − 1 is the second order Legendre polynomial and χτ ,τ + t is the angle between 2 the vectors at times τ and τ + t . where P2 ( x ) =

2

NMR Relaxation Parameters Three relaxation parameters can be measured in an NMR relaxation experiment, NOE, T1 and T2. NOE, or Nuclear Overhauser Effect, is a consequence of dipole–dipole coupling of nuclear spins and can be used in the determination and refinement of structures. T1, or the spin-lattice relaxation time, can be thought of as a characterization of the return of the spin orientations to their equilibrium positions. T2, or the spin-spin relaxation time, describes the change in magnetic field and is related to the exchange of spin between nuclei. T1, T2 and NOE can be determined from the spectral density J (ω) using the following relations: T1 = a ( 3 J ( ω N ) + J ( ωN − ωH ) + 6 J ( ωN + ωH ) )

T2 =

(5)

a ( 4 J ( 0 ) + 3J ( ωN ) + J ( ωN − ωH ) + 6 J ( ωH ) + J ( ωN − ωH ) + 6 J ( ωN + ωH ) ) (6) 2

NOE = 1 +

γH γN

  6 J ( ω N + ω H ) − J ( ω N − ωH )    6 J ( ωN + ωH ) + 3 J ( ωN ) + J ( ω N − ωH ) 

(7)

Here, γ is the gyromagnetic ratio for N and H, and the prefactor a is: 1µ h  γ Nγ H  a=  0 10  4π 2π 

2

(8)

where µ0 is the permeability of vacuum and h is Planck’s constant. When calculating these parameters all constants were chosen as in the experimental set up used by Kelley et. al.80 with a N-H distance of 0.997 Å, ωH = 500.13 MHz and ωN = 160.0 MHz.


Page 11 of 47

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


Generalized Order Parameters For spatially restricted motions, CI(t) will reach a plateau value for large t. This plateau value is the square of the generalized order parameter O and is estimated by:33, 83

= lim→ =

∑〈 , 〉 = ∑"〈 〉 − 1# (9)

where µ are the Cartesian coordinates of the normalized 15N-1H amide spin interaction vector. The corresponding correlation time for the internal dynamics, τ e , is the effective correlation time and is given by

τe =

τ plat 1 ( CI ( 0 ) − CI ( ∞ ) ) dτ CI ( 0 ) − CI ( ∞ ) ∫0

(10)

where τ plat is the time when the plateau is reached. In the diffusion-in-a-cone model the entropy, S, of the N-H vector motion is related to the generalized order parameter, O, according to34

% / '( = ) + ln , -3 − / 1 + 81

(11)

where A is a model-dependent constant, and kB is the Boltzmann constant.

Molecular Dynamics Protocol The results presented here are based on molecular dynamics simulations of reduced GRX1 (PDB ID: 1EGR)78. For GRX1, two trajectories of length 100 ns and 80 ns together with a set of ten 20 ns trajectories have been used for analysis, and in the calculation of the overall rotational tumbling time of the protein we also used two independent trajectories of 50 ns each for the structurally similar reduced Glutaredoxin 3 (Table 1).

84-85

The trajectories of length 100 ns and

80 ns were started from two different NMR conformers, and snapshots from these trajectories are available from the PDB (IDs: 1UPY, 1UPZ,1UQ0, 1UQ1, 1UQ2, 1UQ3, 1UQN, 1UQ6, 1UQ7, and 1UQH). 11 ACS Paragon Plus Environment


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

Page 12 of 47

All trajectories have been simulated using the same protocol and parameters. The setup has been described previously84-85 and is not discussed in detail here. All simulations have been performed with a 2 fs time step and the TIP3P86 water model was used with periodic boundary conditions in a 50x46x42 Å3 box containing 5 sodium ions to make the system electrically neutral The CHARMM2287 force field was used with a 12 Å cut-off for nonbonded interactions, with the atom-based force shift smoothing function.88 For the 20 ns simulations all simulation data have been used since our analysis of sampling showed no significant difference between the first and second halves of these simulations. The overall rotation and translation of the proteins were removed by superimposing each frame of the trajectories on a reference structure, chosen as the initial experimental structure. The oriented trajectories were used in all analyses, except for global tumbling where only the overall translation was removed from the trajectories. All simulations and analyses have been performed with the software package CHARMM.89-90

Analysis For all time averaged analyses, the independent trajectories have been combined into one single set of conformations of total length 380ns. For the Root Mean Square Deviation (RMSD) time series, the 100ns trajectory has been used. The global tumbling was determined using three vectors in different directions, Cα37 − Cα81 , Cα27 − Cα76 and Cα19 − Cα67 which were chosen to span the entire protein and be situated in regions

with low mobility. The auto-correlation function of these vectors was calculated, and their average forms the overall correlation function, Co(t). The generalized order parameters were calculated using equation (9) as implemented in the CHARMM program on the N-H amide vector of all non-proline residues. O was determined by dividing the trajectory into shorter time windows of sizes: 1 ns; 10 ns; 50 ns and the combined total set of conformations of 380 ns. The results were then averaged over all available time windows of the same size. In addition to the direct calculation from the trajectories, O was also determined by first calculating the NMR relaxation parameters T1, T2 and NOE, using the same time windows as in the O calculations, with Equations (5-8) as implemented in CHARMM. 12 ACS Paragon Plus Environment

Page 13 of 47

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


These parameters were then analyzed in terms of the Lipari-Szabo model free approach33, 91 using the Modelfree4 software.92-93 The overall tumbling was represented by setting τ c in equation (3) to the approximate experimental value of 6 ns.80 The effective correlation time, τ e , was calculated in a very similar way as O, but instead using equation (10). The presented values were calculated from the consecutive 100 ns simulation of GRX1. The configurational entropy has been determined with two methods, quasiharmonic analysis24 and through the “diffusion-in-a-cone model”.34 The latter is commonly used to analyze experimental NMR data. The quasiharmonic analysis can be carried out on any part of a protein, while the entropy calculated using the diffusion-in-a-cone model was only obtained for the backbone N-H groups of non-proline residues, since this is what is experimentally available for GRX1.80 The rotational entropy is not considered here since it has been found to account for 1 due to the inclusion of the correlation of motions in SComb, is an estimate of how correlated the motions of the residues in the segment are. To assess the convergence of the calculated quasiharmonic entropy over time and the variations between different simulations, the combined backbone entropy of GRX1 was calculated using time windows of different lengths (Figure 4A). The entropies obtained using different simulations sets (including all ten trajectories in the 10 x 20 ns set) exhibit the same dependence on the length of the time window used for the calculation of the entropy, confirming the equivalence of the simulations in that respect. The variations are slightly larger for SΣ(Res) than for SComb but the trends are the same. The curves also begin to level off after an aggregate of 150 ns. The entropy is proportional to the logarithm of the number of accessible states, and if the protein samples the same number of new states per time interval, a logarithmic increase with time is expected for the entropy. A fit of a logarithmic function (Figure 4A, green curve) to the first 100 ns of the computed entropy shows that it is a conservative estimate in the present case. Assuming a 19 ACS Paragon Plus Environment


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

Page 20 of 47

logarithmic increase, and using our calculated entropy increase from 10 ns to 100 ns in Figure 4A (0.34 kcal/mol/K), the entropy would be expected to increase by 1.7 kcal/mol/K if this trend continues for a 105-fold increase of simulation time (i.e. to ~10ms). When this time scale is approached, the issue becomes more complicated since we are likely to begin sampling also unfolded conformations, whereas our focus here is on the folded state of the protein. The upper and lower curves in Figure 4A have the same shape, and are separated by ~1.5 kcal/mol/K at 380 ns, indicating that the difference between them is converged for backbone conformations within the native state ensemble. In the following entropy analyses, the total 380 ns of GRX1 simulations has been divided into three parts and the average is presented together with the standard deviation. In Table 2, SΣ(Res) and SComb are presented together with their ratio for different atom selections of GRX1. One sees a considerably higher degree of correlated motions for backbone atoms and N-H groups than for the more independently moving side chains. Neglecting the correlations between N-H groups would result in an overestimation of the entropy by about 50%. The sum of backbone and side chains entropy is larger than when all atoms are analyzed together because of correlated motions within the residues. The QHA entropy of all backbone atoms is between 2.23 (SΣ(Res)) and 2.18 (SComb) times higher than for the N-H atoms only, indicating that the two atoms of the N-H vector have a slightly higher average entropy per atom than the other backbone atoms, possibly due to a greater mobility of the amide hydrogen. In Figure 4B the SΣ(Res)/SComb ratio (calculated for the backbone amide groups for different parts of the protein) is presented above the corresponding secondary structure segment. It can be seen that parts with a high degree of correlated motions (high ratio) are parts having low order parameters (indicating long time-scale motions), like the turns at Arg8-Tyr13 and the second αhelix. Parts undergoing fast motions, like the turn at Arg28-Phe31, have a relatively low correlation of motions. Figure 4B also displays a comparison of the entropy calculated with the diffusion-in-cone method and the quasiharmonic entropy of

15

N-1H vectors. To allow comparison of the entropies

calculated with the two methods, the constant A in equation 11 was chosen such that the entropy from the NMR data set was the same as that from the QHA for residue Ala21, which has the 20 ACS Paragon Plus Environment

Page 21 of 47

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


smallest entropy, with a very low error interval (± 3.5x10-4 for the QHA set and ± 3.3x10-4 for the NMR set), in both datasets. Comparing the amide group entropy calculated with the quasiharmonic and diffusion in cone methods, it is evident that they do not give the same relative entropies across the different regions of GRX1. While the quasiharmonic method indicates more entropy in almost all regions compared to the common minimum at Ala21, the deviation between the two methods is greater for some regions. Regions with large scale motions on long time scales (low O2 values with slow convergence) have comparatively higher entropies with the QHA method than with the simpler diffusion in cone model. However, segments with these types of motions also exhibit a high degree of correlations within the segment, which reduces the entropy. This is illustrated when comparing the active site at Arg8-Tyr14 and the turn at Arg28-Phe31. The two regions have similar entropies as predicted by the diffusion in a cone model. Arg8-Tyr14 has larger motions on a slower time scale but the motions are more correlated (SΣ(Res)/SComb = 1.24) than for Arg28Phe31 (SΣ(Res)/SComb=1.14), and the net result is that the entropy per residue calculated with the QHA method for Arg8-Tyr14 is almost twice that for Arg28-Phe31.

Sampling Efficiency – Single vs. Multiple Trajectories Efficient sampling of the conformational space of biological macromolecules is an essential aspect of successful molecular dynamics simulations. The most straightforward way of obtaining better sampling is of course to run longer simulations; other methods include umbrella sampling, replica exchange and running simulations at elevated temperature. Another factor that may play a role is how the total amount of simulation is obtained. One single, long trajectory might not cover the same conformational space as several shorter simulations do together. We have compared the sampling of the conformational space of GRX1 obtained by 100 ns of simulation originating from either a single 100 ns simulation versus parts of ten separate 20 ns simulations. These trajectories have been grouped into five 100 ns segments of data: 1) 1 x 100 ns of consecutive simulation; 2) the first 5 x 20 ns trajectories; 3) the last 5 x 20 ns trajectories; 4) 10 x10 ns from the first half of each 20 ns trajectory; 5) 10 x 10 ns from the last half of each 20 ns trajectory.



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

Page 22 of 47

The conformational sampling has been quantified by assigning each frame of the trajectories to a cluster according to the backbone torsion angles Φ and Ψ of the protein. The structural differences between the obtained clusters mainly lie in mobile regions of turns and the second αhelix, the residues of which exhibit segment RMSDs (calculated using all heavy atoms for superposition) between clusters of typically 3-5 Å but up to 12 Å for some clusters. The clusters have been sorted and numbered according to their population sizes (Figure 5) and the order or numbering does not represent the same structural conformation for the different simulation sets. The distribution width of the clusters is an indication of the conformational space sampled by each 100 ns aggregate of simulation. The 100 ns of data stemming from a single consecutive simulation are distributed into only six clusters while 100 ns coming from five or ten different simulations are distributed into nine to twelve clusters. This indicates a more efficient sampling of possible conformations when running several shorter simulations compared to one long. The simulation segment consisting of the second halves of the ten 20 ns simulation is distributed into 10 clusters while the first halves cover 8 clusters, indicating only a slight improvement of the sampling by allowing an equilibration time of equal length as the production time. A qualitative comparison of conformational sampling of GRX1 obtained with the different simulation lengths and parts can be performed by a pair-wise RMSD analysis of the backbone atoms. The data used for the pair-wise RMSD analyses have, for the sake of clarity, been arranged into the clusters described above and the clusters are arranged in the order of appearance in the original trajectories; note that this removes the time ordering in the presented data (Figure 6). A large number of clusters is not necessarily equivalent to sampling of conformations far apart in conformational space (high RMSD). This can be illustrated by looking at the RMSD within the sets of clusters obtained with the different simulation schemes. By comparing the 6 clusters of the consecutive 100 ns simulation (bottom left corner of Figure 6A, B), the 12 clusters of the first 5 x 20 ns set (Figure 5A upper right corner) and the 9 clusters of the second 5 x 20 ns set (Figure 6B upper right corner), we can see that the distribution of data into many clusters does not automatically mean large RMSD of the sampled conformations. The 12 clusters of the first 5 x 20 ns set display smaller internal RMSD than the 9 clusters of the second 5 x 20 ns set and the 6


Page 23 of 47

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 the consecutive 100 ns display an intra simulation RMSD somewhere in between the two sets of shorter simulations. Figure 6A-B shows that the clusters from the consecutive 100 ns simulations (bottom left corner of Figure 6A-B) and the 5 x 20 ns simulations (upper right corner of Figure 6A-B) show less deviations within themselves than compared to each other (upper left and bottom right corners). This indicates that the differences in simulation lengths result in the sampling of different parts of the conformational space. This view is confirmed in Figure 6C, when the two 5 x 20 sets of simulations are compared. Apart from one apparent outlying cluster (visible at trajectory position 171-183 ns in Figure 6B-C and divided between 160-170 ns and 77 ns in Figure 6D) the deviations are of the same magnitude within the simulation sets as between them. In Figure 6D, the first 10 ns of the ten equivalent simulations (8 clusters, lower left corner of Figure 6D) are compared to the last 10 ns (10 clusters, upper right corner of Figure 6D) of the same simulations. Apart from the outlying cluster, which has been divided between the two sets, no major differences in the sampled parts of the conformations space can be seen for either ends of the simulations, suggesting an equilibrated system fluctuating around an average. The pair wise RMSD data in Figure 6 can also be presented by a 2-D Sammon projection (Figure 7).98, 101 The points are arranged in distinct clusters, confirming the previous results showing that the protein structure in the simulated trajectories is distributed into a few discrete backbone conformations. The regions in Figure 7 not visited in any simulation (white spots) may be energetically unfavorable but could also result from incomplete sampling. Some of the clustered conformations are only visited in one of the three (sets) of trajectories, indicating a different and therefore incomplete sampling by each of the three trajectories. Many of the clusters lie isolated with only few connections to other conformations, suggesting that the conformations are relatively long lived and only visited once, or very few times. This limitation in the sampling of otherwise stable conformations may reflect high energetic barriers.

RMSD and structural dynamics The overall RMSD of the trajectories has been presented previously84 and here we focus on a characterization of the motions seen for structural elements in the 100 ns GRX1 trajectory (Figure 8). The 100 ns trajectory was chosen for its length and more pronounced dynamics 23 ACS Paragon Plus Environment


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

Page 24 of 47

compared to the other trajectories. To study and characterize the dynamics, the GRX1 protein was divided into its secondary structure elements, which are analyzed separately (Figure 8A and B). To visualize the dynamics, the trajectory was divided into 10 ns time segments, from each of which an average structure is calculated. These average structures are then superimposed to illustrate the mobility of each structural segment. The protein was divided into three parts (residues 1-37, containing two parallel β-sheets with an α-helix; 38-60, containing an α-helix surrounded by two turns; 61-85, containing two anti-parallel β-sheets followed by α-helix) and the resulting structures are presented in Figure 8C. A comparison between the experimental structure (inset Figure 1) and the average (over the entire 100 ns trajectory) simulated structure (Figure 8B) shows that the simulated structure does not deviate significantly from the experimental NMR structure over the 100 ns and that the secondary structures are well preserved throughout the simulations. The exception is the last βsheet in GRX1, which loses its inter-strand hydrogen bonds in some instances. The RMSD of individual secondary structures (Figure 8A) shows that most of the structural changes occur in the middle segment and especially in the central α-helix (magenta) and solvent exposed turns (light green) while the other structural elements display limited changes compared to the experimental structure.78 The differences in RMSD are also reflected in the superimposed time-averages displayed in Figure 8C. The increased mobility of turns and unstructured parts compared to α-helices and βsheets, quantified by the order parameter and entropy analysis, are illustrated in these images. What is most striking is the large motions of the α-helix in the middle segment, including residues 38-60. The time scales of the dynamics for different parts of GRX1 revealed in the order parameter and entropy analysis are confirmed here. Residues 1-37 and 61-85, with mainly well converged S2values, show stable RMSDs consistent with a conformation fluctuating around a well defined energy minimum with only very short-lived deviating conformations, and towards the end of the simulation a larger conformational change of the C-terminal α-helix lasting around 15 ns (Figure 8A, black). On the other hand, for the structural elements of the middle part of the protein (yellow, magenta, light green), no well defined average exists and during our 100 ns trajectory, 24 ACS Paragon Plus Environment

Page 25 of 47

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 2nd α-helix visits conformations that lie several Å apart in space. While the conformational changes occur rapidly (< 1ns) the resulting conformations have life times of up to tens of nanoseconds, making them difficult to detect in NMR experiments. These motions also explain the patterns seen in Figure 7 where a few clusters of conformations are visited sequentially with very limited “sampling exchanges” between them. Root mean square fluctuation analysis of backbone atoms in GRX1 (data not shown) confirms the trends described above. Peaks in fluctuations exist for segments of residues in turns and undefined secondary structure elements and for residues in the second α-helix.



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

Page 26 of 47

Discussion The complex behavior of biological macromolecules, with functionally important motions occurring on widely separated timescales, is next to impossible to fully characterize by any single technique, and to approach the goal of a more complete description we have to rely on complementary methods, such as NMR dipolar relaxation and molecular dynamics simulation. These methods are both well-suited to characterize structural dynamics in the nanosecond timescale for molecules of molecular weight up to around 50 kDa. Information obtainable from the two methods partially overlaps, but they do have different limitations. Molecular dynamics simulations can now routinely be performed on the ~µs timescale for a single protein containing a few hundred amino acids. While this yields information about the motions of each individual atom throughout the trajectory of the simulation, it is still just one instance of the trajectory and sampling of some properties may not be fully converged even on this timescale. In the NMR experiments convergence is not an issue (there is a very large number of molecules in the NMRtube), but rather the observable windows in time (up to ≈10 ns) and in structure (individual bond vectors, e.g. backbone N-H groups in a protein) are limiting. We have analyzed these aspects of the dynamics of a relatively small globular protein using molecular dynamics simulations of an aggregate 0.38 µs of simulation time, making two kinds of comparisons: 1/ to published experimental NMR data on the same protein; 2/ between direct analysis of the motions using either all information available in the simulated trajectory or the information that would have been available in an NMR dipolar relaxation experiment. Assessing convergence of molecular dynamics simulations is an important, but non-trivial, issue, which depends on both the system and the property of interest. Configurational entropy is a notoriously ill-behaved property since it depends on the sum of configurations visited during the simulation, a number that will increase with the length of the simulation. In a first approximation the entropy increases with the logarithm of the simulation length. In reality a protein in its native state can only access a finite number of micro states, so the configurational entropy of the native state as probed by a simulation is bounded, and will not increase indefinitely with simulation length; in particular the logarithmic limit implies that simulation times several orders of magnitude shorter than the time scale of protein unfolding will be sufficient to yield meaningful entropy difference estimates (e.g. between different regions in a protein). Below we discuss some 26 ACS Paragon Plus Environment

Page 27 of 47

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


aspects of these issues that we have analyzed from our simulated trajectories of glutaredoxin. The main result is that for this small globular protein we find that secondary structure elements exhibit motions with an amplitude of several Å on a ~10 ns timescale, which require ~100 ns long simulations to yield converged entropy difference estimates in the sense discussed above. These motions do not change the overall fold or the local structure of the active site of the protein. The quantitative details certainly depend on the protein as well as on the force field used, but qualitatively we believe this to be a general result. Furthermore, entropies obtained from NMR backbone N-H generalized order parameters are strongly affected by not including correlation effects, and by being blind to motions on the time scale of these secondary structure element motions. The QHA overestimates entropies when the analysis includes rapid conformational transitions, of the kind exhibited by the secondary structure elements in the simulations. Internal Dynamics. Generalized order parameters calculated directly from the correlation function agree very well with the back calculated O2 from the modelfree analysis of T1, T2, and NOE values from the spectral density function. Residues with converged O2 values above 0.85 are generally found in the middle of α-helices or β-strands, such as Cys14-Glu27 and Tyr72Lys80, engaged in multiple hydrogen-bond interactions with other parts of the protein. These residues have very limited steric freedom and the highly confined environment allows only small scale, high frequency vibrations. The short time scales involved and the absence of any large scale, correlated motions make these regions well described by the order parameters of amide vectors, when a time window smaller than the global tumbling is used (Figure 3). A few regions display low O2 values combined with relatively fast convergence, for instance for the turn at Asp29-Phe31. This combination indicates large motions that occur on rapid time scales. With the Asp29-Phe31 turn (cyan in Figure 8), the large amplitude of the motions is made possible by exposure to the solvent and lack of tight interactions with other parts of the protein. At the same time, the turn is anchored to the more rigid α-helix and β-sheet connected by the turn. This imposes high forces to the vibrational motions of the turn, resulting in short time scales. As illustrated with the entropy Asp29-Phe31 turn, dynamics on fast time scales combined with a low degree of correlated motions can be accurately described by the order parameters of individual amide vectors. 27 ACS Paragon Plus Environment


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

Page 28 of 47

In reduced GRX1, mainly three regions display low, slow converging order parameters, Arg8Tyr13, Arg39-Val59 and the C-terminal segment. The region Arg8-Tyr13 is, like the turn at Asp29-Phe31, largely exposed to the solvent and lacks contacts to other parts of the protein, allowing greater dynamics. But unlike Arg8-Tyr13, Asp29-Phe31 is not rigidly anchored on the C-terminal side, which allows motions on longer time scales. The most dynamic region of GRX1 is Arg39-Val59, which includes an α-helix between two turns. Here, the turns at Arg40-Gly42 and Val55-Glu56 surrounding the 2nd α-helix work as hinges,102 allowing the α-helix to move without losing its own internal structure. As seen in Figure 8A, these conformational changes are rare and only occur a few times during a 100 ns trajectory. The C-terminal end undergoes considerably more motions than the N-terminal end and is almost completely disordered (O2 = 0) on long time scales. This can be explained by the increased flexibility given by the turn at Glu81Asn82 and the undefined secondary structure, with significantly less tertiary interactions, of the C-terminal residues. It is clear that the study of individual amide vectors on timescales less than the global tumbling can not fully describe rare conformational changes or motions on slow time scales. These motions will remain undetected as is evident from the experimental O2-values in Figure 3. Entropy. The quasi harmonic entropies of the amide vectors are overall larger, and also show a larger variation along the peptide sequence, than the values obtained with the diffusion-in-a-cone model. This underestimation of entropy relative to QHA is, at least in part, due to slow motions with relatively large amplitude for local regions in the structure (Figure 8). This is for example seen for the turn at Arg8-Tyr13 and the hinge at Arg39-Val59. These regions have similar O2 values and hence similar entropies as predicted by the diffusion-in-a-cone model, but the slow large scale motions within the Arg8-Tyr13 and Arg39-Val59 regions result in an entropy increase from the minimum at Ala21 which is more than twice as large with the QHA model as with the diffusion-in-a-cone model. However, these regions are the ones with the most correlated motions (the ratio SΣ(Res)/SComb is the largest here), and hence methods that do not account for correlations will overestimate the entropy. It seems that with the diffusion-in-a-cone model, the underestimation of the entropy for residues involved in long time scale motions is, to some degree, compensated by the overestimation of correlated motions to the benefit of. One should also be aware that QHA may overestimate the entropy of regions with several energy minima.25 28 ACS Paragon Plus Environment

Page 29 of 47

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


Sampling. We have compared the sampling efficiency of a total of 100 ns of molecular dynamics trajectories obtained from simulations of different individual lengths, including variants keeping only the first or second half of the same simulations. To quantitatively asses the amount of sampling produced by the different approaches, we clustered the conformations of the different sampling aggregates according to the Φ and Ψ backbone torsions (Figure 5). While clustering does not quantify the absolute amount of sampling, it is a good method to compare the sampling across protocols.68 In addition, to get a qualitative picture of the sampled conformations, we have performed a pair-wise RMSD analysis of the backbone atoms between trajectory frames (Figure 6). This allowed us to estimate the similarity, in terms of position in space, between distinct timewindows of one simulation of between different simulations. We also investigated how the sampled conformations are visited as a function of time and how much exchange there is between those during the trajectories (Figure 7). This is done by projecting the RMSD of time dependent trajectories into 2 dimensions and connecting the data point by order in time. Figure 5 shows that allocating the total simulation time to multiple trajectories yields a larger amount of sampled conformations than if allocated to a single simulation. This is true for the division into 5 x 20 ns as well as 10 x 10 ns of simulations. Keeping in mind that the clustering of backbone torsions is not directly correlated to a wide distribution in conformational space (RMSD) and that only a small part of the total conformational space is visited in simulations of this time scale, a clear increase in relative sampling is observed with multiple simulations. Judging from the number of clusters, up to 100% more sampling can be obtained for the same computational cost by running multiple simulations with different starting velocities. Taking into account the cost advantages of parallelized computing compared to linear single trajectories, the benefits become even clearer. One may envisage that even more efficient sampling would be achieved by using slightly different starting conformers,18 for instance drawn from an NMR structure bundle. The limit for how finely a protein simulation may be divided into several separate trajectories is governed by the time scale of the processes studied. Each simulation must be at least longer than the time scale of the slowest process studied. Even the motions of single side chains may be quite complicated. For example, Arg8 in GRX1 may modulate the pKa of the reactive Cys11 by coming close to Cys11 and stabilize its thiolate form. In the 100ns simulation it takes ~50 ns 29 ACS Paragon Plus Environment


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

Page 30 of 47

before this happens, and while the final transition is very fast, it is preceded for ~10 ns by a state were Arg8 has approached Cys11, but not yet made contact. The present characterization of the motions of GRX1 shows that, although many conformational changes occur infrequently, the actual process of change is generally fast. This is illustrated by examining the large scale changes of the Ile38-Pro60 segment depicted in Figure 8A and 8C-2. Most conformational changes of this region occur on timescale of a few nanoseconds but last for order of magnitudes longer. This feature of the studied protein is illustrated in a more general way in Figure 7. In the 2-D projection of the RMSD between states it is clear that the trajectories mainly consist of well defined clusters of many points in series (long times) connected by a few points (short times) between them.

Supporting Information A table of the best fitting model for each residue from the Modelfree4 analysis of the GRX1 simulations is available as Supporting Information. This material is available free of charge via the Internet at http://pubs.acs.org.


Page 31 of 47

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


Tables Protein GRX1 GRX3

Simulation time (ns) a 100 + 82 + 10x20 b 2x50 2

Table 1. Length of simulations analyzed in this paper. a The total 380 ns of simulations was used when calculating order parameters and entropy. b when calculating the global tumbling an additional 100 ns data of the similarly structured protein GRX3 was used.

SΣ(Res) SComb [kcal/mol/K] All Atoms 7.92 ± 0.13 5.90 ± 0.07 Backbone 3.44 ± 0.09 1.97 ± 0.04 Side chains 5.70 ± 0.07 4.64 ± 0.04 1.54 ± 0.04 0.90 ± 0.02 N-H N/A N-H (NMR) 1.36 ± 0.02

SΣ(Res)/SComb 1.34 1.75 1.23 1.78 N/A

Table 2. Total entropy of different parts of GRX1 calculated from all available simulation data as the sum of individual residues and the combined total. The quasiharmonic entropy has been calculated for different parts of the protein while the NMR entropy (diffusion in a cone model) is presented only for the N-H vector. The NMR entropy is obtained setting the constant A in equation 11 to a value such that the entropy for the residue with the smallest entropy using equation 11 is the same as the quasiharmonic N-H entropy. The total 380 ns simulation data have been divided into three equally long sets and the presented values are the average and standard deviations of the parts.



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

Page 32 of 47

Figures

Figure 1. Correlation function of global orientation of GXR1 and GRX3. Black curve is combined average of all long simulations of GRX1 and GRX3. Red curve is the fitted first order exponential 2 = 0.94 7 8/".9 + 0.28 . Inset picture shows GRX1 with the vectors used for determining the global tumbling in magenta and residues, for which correlation functions are presented in Figure 2, in green stick representation.


Page 33 of 47

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


Figure 2. A) Examples of internal correlation functions of residues with short (green line, Asp29, τe=0.15ns), long (blue line, Leu83, τe=0.1.9 ns) and unconverged (red line, Thr44) effective correlation times; see inset in Figure 1 for location of these residues within the GRX1 structure. B) Overall rotational tumbling (black line, from exponential fit in figure 1); internal correlation functions for Asp29 (solid green) and Thr44 (solid red); total correlation functions calculated directly from the trajectory before removal of overall rotation for Asp29 (green dots) and Thr44 (red dots); total correlation function obtained as the product of the overall tumbling and the internal correlation function for Asp29 (green dashes) and Thr44 (red dashes).



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

Figure 3. Comparison of order parameters (O2) per residue of GRX1 for different time windows calculated directly from simulation (A) or back calculated with Modelfree4 (B). Magenta is the average experimental values with associated standard deviation. Black curves are obtained over 1 ns time windows, red over 10 ns, green over 50 ns and blue over all 380 ns.


Page 34 of 47

Page 35 of 47

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


A

B

B



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

Figure 4. A) Quasi-harmonic entropy of the backbone atoms as function the length of the timewindow used for computing the entropy for three simulations. Black: single 82 ns. Red: single 100ns; Blue: combined 10 x 20 ns. The lower curves (squares) are SComb and the upper (circles) are SΣ(Res). The last points at 380 ns in the blue curves are the combined data of all simulations. Error bars are standard deviations between individual time windows. The green line is the function 0.60+0.36ln(t) obtained after fitting to the points in the interval 1ns – 100ns of the blue curve. B). Entropy per residue of GRX1. The two top curves are quasiharmonic entropy for all backbone atoms (black) and 15N-1H group (red) in each residue. The blue curve is the entropy of the 15N-1H group estimate by the diffusion in a cone model. The entropy is presented per residue and to allow comparison, the minimum value of each dataset (Ala21 for both) has been translated to zero. Each curve is obtained as the average of three equally sized parts of the available data and error bars were estimated by the standard deviation between them. The bottom values are the ratios of entropy calculated individually per residue in relation to the summed entropy of the segment as a whole, SΣ(Res)/SComb.


Page 36 of 47

Page 37 of 47

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


Figure 5. Number of conformational clusters produced by 100 ns sampling obtained by different simulation lengths and parts. Black bars refers to one continuous 100 ns simulation, red and blue to two equivalent 5×20 ns sets and green and magenta refers to the first and second half of the 10×20 ns sets respectively. The clusters have been sorted and numbered according to their respective population size and the numbering does not represent any particular structural conformation.



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

Page 38 of 47

Figure 6. Pairwise RMSD between conformations produced by 100 ns of sampling obtained by different simulation lengths and parts. The compared simulations in each matrix are: A) 1×100 ns– 5×20 ns (1); B) 1×100 ns – 5×20 ns (2); C) 5×20 ns (1) - 5×20 ns (2); D) 10×10 ns (1st-half) 10×10 ns (2nd-half). The matrices are organized with the first 100 ns trajectory at 0-100 and the second, comparison trajectory following at 100-200. The diagonal from (0,0) to (200,200) is the symmetry line with complementary pairs on each side.


Page 39 of 47

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


Figure 7. 2-D Sammon projection of backbone RMSD for three continuous simulations. The points are colored according the trajectory, red: 100 ns, green: 82 ns and blue: 10×20 ns. Points consecutive in time are connected by lines.



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

Page 40 of 47

Figure 8. Root Mean Square Deviation of secondary structure elements (segment RMSD) in relation to the experimental structure37 and dynamics of structural segments of GRX1. A) RMSD time series from the 100 ns simulation of backbone atoms of GRX1. The protein has been divided into its secondary structures, which are colored according to the sequence in the legend. B) Average simulated structure colored according to the legend in A C) Dynamics of structural segments of GRX1 visualized by ten average structures (from 10 ns each) segments overlaid on top of each other. The three pictures contain residues: 1-37 (C-1); 38-60 (C-2); 38-60 (C-3).


Page 41 of 47

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


References 1. Campbell, I. D.; Dobson, C. M.; Moore, G. R.; Perkins, S. J.; Williams, R. J. P., Temperature Dependent Molecular Motion of a Tyrosine Residue of Ferrocytochrome C. Febs Letters 1976, 70, 96-100. 2. Huber, R.; Deisenhofer, J.; Colman, P. M.; Matsushima, M.; Palm, W., Crystallographic Structure Studies of an Igg Molecule and an Fc Fragment. Nature 1976, 264, 415-420. 3. Perutz, M. F., Stereochemistry of Cooperative Effects in Haemoglobin: Haem-Haem Interaction and the Problem of Allostery. Nature 1970, 228, 726-734. 4. Karplus, M.; McCammon, J. A., Molecular Dynamics Simulations of Biomolecules. Nature Struct. Biol. 2002, 9, 646-652. 5. Tama, F.; Brooks, C. L., Symmetry, Form, and Shape: Guiding Principles for Robustness in Macromolecular Machines. Annu. Rev. Biophys. Biomol. Struct. 2006, 35, 115-133. 6. Tokuriki, N.; Tawfik, D. S., Protein Dynamism and Evolvability. Science 2009, 324, 203207. 7. Stone, M. J., Nmr Relaxation Studies of the Role of Conformational Entropy in Protein Stability and Ligand Binding. Acc. Chem. Res. 2001, 34, 379-388. 8. Gibbs, A. C., Elements and Modulation of Functional Dynamics. J. Med. Chem. 2014. 9. Zhou, H.-X.; Gilson, M. K., Theory of Free Energy and Entropy in Noncovalent Binding. Chemical Reviews 2009, 109, 4092-4107. 10. Straatsma, T. P.; Mccammon, J. A., Computational Alchemy. Annu. Rev. Phys. Chem. 1992, 43, 407-435. 11. Roux, B., The Calculation of the Potential of Mean Force Using Computer Simulations. Comp. Phys. Comm. 1995, 91, 275-282. 12. Deng, Y.; Roux, B., Computations of Standard Binding Free Energies with Molecular Dynamics Simulations. J. Phys. Chem. B 2009, 113, 2234-2246. 13. Foloppe, N.; Hubbard, R., Towards Predictive Ligand Design with Free-Energy Based Computational Methods? Curr. Med. Chem. 2006, 13, 3583-3608. 14. Åqvist, J.; Medina, C.; Samuelsson, J. E., New Method for Predicting Binding-Affinity in Computer-Aided Drug Design. Prot. Eng. 1994, 7, 385-391. 15. Åqvist, J.; Luzhkov, V. B.; Brandsdal, B. O., Ligand Binding Affinities from Md Simulations. Acc. Chem. Res. 2002, 35, 358-365. 16. Gohlke, H.; Kiel, C.; Case, D. A., Insights into Protein-Protein Binding by Binding Free Energy Calculation and Free Energy Decomposition for the Ras-Raf and Ras-Ralgds Complexes. J. Mol. Biol. 2003, 330, 891-913. 17. Kollman, P. A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S.; Chong, L.; Lee, M.; Lee, T.; Duan, Y.; Wang, W.; et al., Calculating Structures and Free Energies of Complex Molecules: Combining Molecular Mechanics and Continuum Models. Acc. Chem. Res. 2000, 33, 889-897. 18. Genheden, S.; Diehl, C.; Akke, M.; Ryde, U., Starting-Condition Dependence of Order Parameters Derived from Molecular Dynamics Simulations. J. Chem. Theory Comp. 2010, 6, 2176-2190. 19. Genheden, S.; Ryde, U., How to Obtain Statistically Converged Mm/Gbsa Results. J. Comp. Chem. 2010, 31, 837-846. 41 ACS Paragon Plus Environment


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

20. Roche, J.; Caro, J. A.; Norberto, D. R.; Barthe, P.; Roumestand, C.; Schlessman, J. L.; Garcia, A. E.; García-Moreno E., B.; Royer, C. A., Cavities Determine the Pressure Unfolding of Proteins. Proc. Natl. Acad. Sci. USA 2012, 109, 6945-6950. 21. Karplus, M.; Kushick, J. N., Method for Estimating the Configurational Entropy of Macromolecules. Macromolecules 1981, 14, 325-332. 22. Teeter, M. M.; Case, D. A., Harmonic and Quasiharmonic Descriptions of Crambin. J. Phys. Chem 1990, 94, 8091-8097. 23. Brooks, B. R.; Janezic, D.; Karplus, M., Harmonic-Analysis of Large Systems .1. Methodology. J. Comp. Chem. 1995, 16, 1522-1542. 24. Andricioaei, I.; Karplus, M., On the Calculation of Entropy from Covariance Matrices of the Atomic Fluctuations. J. Chem. Phys. 2001, 115, 6289-6292. 25. Chang, C. E.; Chen, W.; Gilson, M. K., Evaluating the Accuracy of the Quasiharmonic Approximation. J. Chem. Theory Comp. 2005, 1, 1017-1028. 26. Hnizdo, V.; Darian, E.; Fedorowicz, A.; Demchuk, E.; Li, S.; Singh, H., NearestNeighbor Nonparametric Method for Estimating the Configurational Entropy of Complex Molecules. J. Comp. Chem. 2007, 28, 655-668. 27. Li, D.-W.; Brüschweiler, R., In Silico Relationship between Configurational Entropy and Soft Degrees of Freedom in Proteins and Peptides. Physical Review Letters 2009, 102, 118108. 28. Baron, R.; van Gunsteren, W. F.; Hünenberger, P. H., Estimating the Configurational Entropy from Molecular Dynamics Simulations: Anharmonicity and Correlation Corrections to the Quasi-Harmonic Approximation. Trends Phys. Chem. 2006, 11, 87-122. 29. Hnizdo, V.; Tan, J.; Killian, B. J.; Gilson, M. K., Efficient Calculation of Configurational Entropy from Molecular Simulations by Combining the Mutual-Information Expansion and Nearest-Neighbor Methods. J. Comp. Chem. 2008, 29, 1605-1614. 30. Hensen, U.; Lange, O. F.; Grubmüller, H., Estimating Absolute Configurational Entropies of Macromolecules: The Minimally Coupled Subspace Approach. PLoS ONE 2010, 5, e9179. 31. Sharp, K., Calculation of Molecular Entropies Using Temperature Integration. J. Chem. Theory Comp. 2012, 9, 1164-1172. 32. Halle, B.; Wennerström, H., Interpretation of Magnetic Resonance Data from Water Nuclei in Heterogeneous Systems. J. Chem. Phys. 1981, 75, 1928-1943. 33. Lipari, G.; Szabo, A., Model-Free Approach to the Interpretation of Nuclear MagneticResonance Relaxation in Macromolecules .1. Theory and Range of Validity. J. Am. Chem. Soc. 1982, 104, 4546-4559. 34. Yang, D. W.; 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. 35. Akke, M.; Brüschweiler, R.; Palmer Iii, A. G., Nmr Order Parameters and Free Energy: An Analytical Approach and Its Application to Cooperative Ca2+ Binding by Calbindin D9k. J. Am. Chem. Soc. 1993, 115, 9832-9833. 36. Yang, D. W.; Mok, Y. K.; FormanKay, J. D.; Farrow, N. A.; Kay, L. E., Contributions to Protein Entropy and Heat Capacity from Bond Vector Motions Measured by Nmr Spin Relaxation. J. Mol. Biol. 1997, 272, 790-804. 37. Li, Z.; Raychaudhuri, S.; Wand, A. J., Insights into the Local Residual Entropy of Proteins Provided by Nmr Relaxation. Protein Science 1996, 5, 2647-2650. 38. Spyracopoulos, L.; Sykes, B. D., Thermodynamic Insights into Proteins from Nmr Spin Relaxation Studies. Curr. Op. Struct. Biol. 2001, 11, 555-559. 42 ACS Paragon Plus Environment

Page 42 of 47

Page 43 of 47

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


39. Idiyatullin, D.; Daragan, V. A.; Mayo, K. H., 15nh Backbone Dynamics of Protein Gb1: Comparison of Order Parameters and Correlation Times Derived Using Various "Model-Free" Approaches. J. Phys. Chem. B 2003, 107, 2602-2609. 40. LeMaster, D. M., Nmr Relaxation Order Parameter Analysis of the Dynamics of Protein Side Chains. J. Am. Chem. Soc. 1999, 121, 1726-1742. 41. Tugarinov, V.; Liang, Z.; Shapiro, Y. E.; Freed, J. H.; Meirovitch, E., A Structural ModeCoupling Approach to 15n Nmr Relaxation in Proteins. J. Am. Chem. Soc. 2001, 123, 3055-63. 42. Brüschweiler, R.; Wright, P. E., Nmr Order Parameters of Biomolecules: A New Analytical Representation and Application to the Gaussian Axial Fluctuation Model. J. Am. Chem. Soc. 1994, 116, 8426-8427. 43. Best, R. B.; Vendruscolo, M., Determination of Protein Structures Consistent with Nmr Order Parameters. J. Am. Chem. Soc. 2004, 126, 8090-8091. 44. Chen, J. H.; Brooks, C. L.; Wright, P. E., Model-Free Analysis of Protein Dynamics: Assessment of Accuracy and Model Selection Protocols Based on Molecular Dynamics Simulation. J. Biomol. NMR 2004, 29, 243-257. 45. Nederveen, A. J.; Bonvin, A. M. J. J., Nmr Relaxation and Internal Dynamics of Ubiquitin from a 0.2 Mu S Md Simulation. J. Chem. Theory Comp. 2005, 1, 363-374. 46. Chandrasekhar, I.; Clore, G. M.; Szabo, A.; Gronenborn, A. M.; Brooks, B. R., A 500-Ps Molecular-Dynamics Simulation Study of Interleukin-1-Beta in Water - Correlation with Nuclear-Magnetic-Resonance Spectroscopy and Crystallography. Journal of Molecular Biology 1992, 226, 239-250. 47. Eriksson, M. A. L.; Berglund, H.; Härd, T.; Nilsson, L., A Comparison of Nitrogen-15 Nmr Relaxation Measurements with a Molecular Dynamics Simulation: Backbone Dynamics of the Glucocorticoid Receptor DNA-Binding Domain. Proteins 1993, 17, 375-90. 48. Lipari, G.; Szabo, A.; Levy, R. M., Protein Dynamics and Nmr Relaxation: Comparison of Simulations with Experiment. Nature 1982, 300, 197-198. 49. Olejniczak, E. T.; Dobson, C. M.; Karplus, M.; Levy, R. M., Motional Averaging of Proton Nuclear Overhauser Effects in Proteins - Predictions from a Molecular-Dynamics Simulation of Lysozyme. J. Am. Chem. Soc. 1984, 106, 1923-1930. 50. Post, C. B., Internal Motional Averaging and 3-Dimensional Structure Determination by Nuclear-Magnetic-Resonance. Journal of Molecular Biology 1992, 224, 1087-1101. 51. Levy, R. M.; Karplus, M.; Wolynes, P. G., Nmr Relaxation Parameters in Molecules with Internal Motion - Exact Langevin Trajectory Results Compared with Simplified Relaxation Models. J. Am. Chem. Soc. 1981, 103, 5998-6011. 52. Buck, M.; Bouguet-Bonnet, S.; Pastor, R. W.; MacKerell, A. D., Jr., Importance of the Cmap Correction to the Charmm22 Protein Force Field: Dynamics of Hen Lysozyme. Biophys. J. 2006, 90, L36-38. 53. Showalter, S. A.; Brüschweiler, R., Validation of Molecular Dynamics Simulations of Biomolecules Using Nmr Spin Relaxation as Benchmarks: Application to the Amber99sb Force Field. J. Chem. Theory Comp. 2007, 3, 961-975. 54. Markwick, P. R. L.; Bouvignies, G.; Blackledge, M., Exploring Multiple Timescale Motions in Protein Gb3 Using Accelerated Molecular Dynamics and Nmr Spectroscopy. J Am Chem Soc 2007, 129, 4724-4730. 55. 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. 43 ACS Paragon Plus Environment


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

56. Villa, A.; Stock, G., What Nmr Relaxation Can Tell Us About the Internal Motion of an Rna Hairpin: A Molecular Dynamics Simulation Study. J. Chem. Theory Comp. 2006, 2, 12281236. 57. Duchardt, E.; Nilsson, L.; Schleucher, J., Cytosine Ribose Flexibility in DNA: A Combined Nmr 13c Spin Relaxation and Molecular Dynamics Simulation Study. Nucl. Acids Res. 2008, 36, 4211-4219. 58. Musselman, C.; Zhang, Q.; Al-Hashimi, H.; Andricioaei, I., Referencing Strategy for the Direct Comparison of Nuclear Magnetic Resonance and Molecular Dynamics Motional Parameters in Rna. J. Phys. Chem. B 2009, 114, 929-939. 59. Maragakis, P.; Lindorff-Larsen, K.; Eastwood, M. P.; Dror, R. O.; Klepeis, J. L.; Arkin, I. T.; Jensen, M. O.; Xu, H.; Trbovic, N.; Friesner, R. A.; Palmer, A. G.; Shaw, D. E., 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. 60. Halle, B., The Physical Basis of Model-Free Analysis of Nmr Relaxation Data from Proteins and Complex Fluids. J Chem Phys 2009, 131, 224507. 61. Wong, V.; Case, D. A., Evaluating Rotational Diffusion from Protein Md Simulations. J. Phys. Chem. B 2008, 112, 6013-6024. 62. Macek, P.; Novak, P.; Zidek, L.; Sklenar, V., Backbone Motions of Free and PheromoneBound Major Urinary Protein I Studied by Molecular Dynamics Simulation. J. Phys. Chem. B 2007, 111, 5731-5739. 63. Koller, A. N.; Schwalbe, H.; Gohlke, H., Starting Structure Dependence of Nmr Order Parameters Derived from Md Simulations: Implications for Judging Force-Field Quality. Biophysical Journal 2008, 95, L4-L6. 64. Clarage, J. B.; Romo, T.; Andrews, B. K.; Pettitt, B. M.; Phillips, G. N., A Sampling Problem in Molecular Dynamics Simulations of Macromolecules. Proc. Natl. Acad. Sci. USA 1995, 92, 3288-3292. 65. Genheden, S.; Akke, M.; Ryde, U., Conformational Entropies and Order Parameters: Convergence, Reproducibility, and Transferability. J. Chem. Theory Comp. 2014, 10, 432-438. 66. Genheden, S.; Ryde, U., Will Molecular Dynamics Simulations of Proteins Ever Reach Equilibrium? Phys. Chem. Chem. Phys. 2012, 14, 8662-8677. 67. Yang, W.; Nymeyer, H.; Zhou, H.-X.; Berg, B.; Brüschweiler, R., Quantitative Computer Simulations of Biomolecules: A Snapshot. J. Comp. Chem. 2008, 29, 668-672. 68. Zuckerman, D. M., Equilibrium Sampling in Biomolecular Simulations. Annu. Rev. Biophys. 2011, 40, 41-62. 69. Caves, L. S. D.; Evanseck, J. D.; Karplus, M., Locally Accessible Conformations of Proteins: Multiple Molecular Dynamics Simulations of Crambin. Protein Sci. 1998, 7, 649-666. 70. Lahiri, A.; Nilsson, L., Examining the Characteristics of Chaos in Biomolecular Dynamics: A Random Matrix Approximation. Chem. Phys. Lett. 1999, 311, 459-466. 71. Andrews, B. K.; Romo, T.; Clarage, J. B.; Pettitt, B. M.; Phillips, G. N., Characterizing Global Substates of Myoglobin. Struct. Fold. Des. 1998, 6, 587-594. 72. Genheden, S.; Ryde, U., A Comparison of Different Initialization Protocols to Obtain Statistically Independent Molecular Dynamics Simulations. J. Comp. Chem. 2011, 32, 187-195. 73. Elofsson, A.; Nilsson, L., How Consistent Are Molecular Dynamics Simulations Comparing Structure and Dynamics in Reduced and Oxidized Escherichia-Coli Thioredoxin. J. Mol. Biol. 1993, 233, 766-780. 44 ACS Paragon Plus Environment

Page 44 of 47

Page 45 of 47

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


74. Straub, J. E.; Rashkin, A. B.; Thirumalai, D., Dynamics in Rugged Energy Landscapes with Applications to the S-Peptide and Ribonuclease A. J Am Chem Soc 1994, 116, 2049-2063. 75. Liu, L.; Gronenborn, A. M.; Bahar, I., Longer Simulations Sample Larger Subspaces of Conformations While Maintaining Robust Mechanisms of Motion. Proteins: Structure, Function, and Bioinformatics 2012, 80, 616-625. 76. Roy, J.; Laughton, C. A., Long-Timescale Molecular-Dynamics Simulations of the Major Urinary Protein Provide Atomistic Interpretations of the Unusual Thermodynamics of Ligand Binding. Biophys. J. 2010, 99, 218-226. 77. Nilsson, L.; Karshikoff, A., Multiple Ph Regime Molecular Dynamics Simulation for Pk Calculations. PLoS ONE 2011, 6, e20116. 78. Sodano, P.; Xia, T. H.; Bushweller, J. H.; Bjornberg, O.; Holmgren, A.; Billeter, M.; Wuthrich, K., Sequence-Specific H-1-Nmr Assignments and Determination of the 3-Dimensional Structure of Reduced Escherichia-Coli Glutaredoxin. J. Mol. Biol. 1991, 221, 1311-1324. 79. Martin, J. L., Thioredoxin - a Fold for All Reasons. Structure 1995, 3, 245-250. 80. Kelley, J. J.; Caputo, T. M.; Eaton, S. F.; Laue, T. M.; Bushweller, J. H., Comparison of Backbone Dynamics of Reduced and Oxidized Escherichia Coli Glutaredoxin-1 Using N-15 Nmr Relaxation Measurements. Biochemistry 1997, 36, 5029-5044. 81. Keeler, J., Understanding Nmr Spectroscopy. Wilye: Chichester, England, 2005. 82. Henzler-Wildman, K. A.; Lei, M.; Thai, V.; Kerns, S. J.; Karplus, M.; Kern, D., A Hierarchy of Timescales in Protein Dynamics Is Linked to Enzyme Catalysis. Nature 2007, 450, 913-U27. 83. Henry, E. R.; Szabo, A., Influence of Vibrational Motion on Solid-State Line-Shapes and Nmr Relaxation. J Chem Phys 1985, 82, 4753-4761. 84. Foloppe, N.; Nilsson, L., The Glutaredoxin -C-P-Y-C- Motif: Influence of Peripheral Residues. Structure 2004, 12, 289-300. 85. Foloppe, N.; Sagemark, J.; Nordstrand, K.; Berndt, K. D.; Nilsson, L., Structure, Dynamics and Electrostatics of the Active Site of Glutaredoxin 3 from Escherichia Coli: Comparison with Functionally Related Proteins. J. Mol. Biol. 2001, 310, 449-470. 86. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J.; Impey, R. W.; Klein, M. L., Comparison of Simple Potential Functions for Simulating Liquid Water. J Chem Phys 1983, 79, 926-935. 87. MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; et al., All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586-3616. 88. Steinbach, P. J.; Brooks, B. R., New Spherical-Cutoff Methods for Long-Range Forces in Macromolecular Simulation. J. Comp. Chem. 1994, 15, 667-683. 89. Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M., Charmm: A Program for Macromolecular Energy, Minimization, and Dynamics Calculations. J. Comp. Chem. 1983, 4, 187-217. 90. Brooks, B. R.; Brooks, C. L.; Mackerell, A. D.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; et al., Charmm: The Biomolecular Simulation Program. J. Comp. Chem. 2009, 30, 1545-1614. 91. Lipari, G.; Szabo, A., Model-Free Approach to the Interpretation of Nuclear MagneticResonance Relaxation in Macromolecules .2. Analysis of Experimental Results. J. Am. Chem. Soc. 1982, 104, 4559-4570.



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

Page 46 of 47

92. Palmer, A. G.; Rance, M.; Wright, P. E., Intramolecular Motions of a Zinc Finger DNABinding Domain from Xfin Characterized by Proton-Detected Natural Abundance C-12 Heteronuclear Nmr-Spectroscopy. J. Am. Chem. Soc. 1991, 113, 4371-4380. 93. Mandel, A. M.; Akke, M.; Palmer, A. G., Backbone Dynamics of Escherichia-Coli Ribonuclease Hi - Correlations with Structure and Function in an Active Enzyme. J. Mol. Biol. 1995, 246, 144-163. 94. Wrabl, J. O.; Shortle, D.; Woolf, T. B., Correlation between Changes in Nuclear Magnetic Resonance Order Parameters and Conformational Entropy: Molecular Dynamics Simulations of Native and Denatured Staphylococcal Nuclease. Prot. Struct. Func. Gen. 2000, 38, 123-133. 95. Karpen, M. E.; Tobias, D. J.; Brooks, C. L., Statistical Clustering-Techniques for the Analysis of Long Molecular-Dynamics Trajectories - Analysis of 2.2-Ns Trajectories of Ypgdv. Biochem. 1993, 32, 412-420. 96. Carpenter, G. A.; Grossberg, S., Art-2 - Self-Organization of Stable Category Recognition Codes for Analog Input Patterns. App. Opt. 1987, 26, 4919-4930. 97. Levitt, M., Molecular-Dynamics of Native Protein .2. Analysis and Nature of Motion. J. Mol. Biol. 1983, 168, 621-657. 98. Sammon, J. W., A Nonlinear Mapping for Data Structure Analysis. IEEE. Trans. Comput. 1969, C 18, 401-409. 99. Stone, M. J.; Chandrasekhar, K.; Holmgren, A.; Wright, P. E.; Dyson, H. J., Comparison of Backbone and Tryptophan Side-Chain Dynamics of Reduced and Oxidized Escherichia-Coli Thioredoxin Using N-15 Nmr Relaxation Measurements. Biochem. 1993, 32, 426-435. 100. Mark, P.; Nilsson, L., Structure and Dynamics of the Tip3p, Spc and Spc/E Water Models at 298k. J. Phys. Chem. A 2001, 105, 9954-9960. 101. Levitt, M., Computer Simulation of DNA Double Helix Dynamics. Cold Spring Harbor Symp. Quant. Biol. 1983, 47, 251-261. 102. Johnson, E.; Showalter, S. A.; Bruschweiler, R., A Multifaceted Approach to the Interpretation of Nmr Order Parameters: A Case Study of a Dynamic α−Helix. J. Phys. Chem. B 2008, 112, 6203-6210.


Page 47 of 47

Table of content graphic

3.4 3.2 3.0

RMSD Residue:60--80

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.8 2.6 2.4 2.2 2.0 1.8 1.6 1.4 1.2 1.0 0.8 0.6 0.4 0

20

40

60

80

100

Time (ns)


General order parameter based correlation analysis of protein backbone motions between experimental NMR relaxation measurements and molecular dynamics simulations.

Charge templates in aromatic plus ionic liquid systems revisited: NMR experiments and molecular dynamics simulations.

Slow motions in microcrystalline proteins as observed by MAS-dependent 15N rotating-frame NMR relaxation.

Gorge Motions of Acetylcholinesterase Revealed by Microsecond Molecular Dynamics Simulations.

Collective motions in proteins: a covariance analysis of atomic fluctuations in molecular dynamics and normal mode simulations.

Nonadiabatic molecular dynamics simulations: synergies between theory and experiments.

Local isotropic diffusion approximation for coupled internal and overall molecular motions in NMR spin relaxation.

Polymer Brushes under Shear: Molecular Dynamics Simulations Compared to Experiments.

Molecular Dynamics and NMR Shed Light on Motions Underpinning Dynamical Transitions in Biomolecules.

Concerted dynamic motions of an FABP4 model and its ligands revealed by microsecond molecular dynamics simulations.

Consecutive non-natural PZ nucleobase pairs in DNA impact helical structure as seen in 50 μs molecular dynamics simulations.

Reply to Jensen and Blackledge: Dual quantifications of intrinsically disordered proteins by NMR ensembles and molecular dynamics simulations.

Conformation of endothelin in aqueous ethylene glycol determined by 1H-NMR and molecular dynamics simulations.

Molecular dynamics simulations in biology.

Lipid Models for United-Atom Molecular Dynamics Simulations of Proteins.

Trypsinogen activation as observed in accelerated molecular dynamics simulations.

Distance-Based Configurational Entropy of Proteins from Molecular Dynamics Simulations.

Molecular dynamics study of secondary structure motions in proteins: application to myohemerythrin.

Molecular dynamics simulations of the STAS domains of rat prestin and human pendrin reveal conformational motions in conserved flexible regions.

Molecular dynamics simulations of proton transverse relaxation times in suspensions of magnetic nanoparticles.

The liquid surface of chiral ionic liquids as seen from molecular dynamics simulations combined with intrinsic analysis.

Searching for protein binding sites from Molecular Dynamics simulations and paramagnetic fragment-based NMR studies.

Interactions of peripheral proteins with model membranes as viewed by molecular dynamics simulations.

Solution structure of the Lewis x oligosaccharide determined by NMR spectroscopy and molecular dynamics simulations.