Ab-initio Prediction of NMR Spin-Relaxation Parameters from

Subscriber access provided by UNIVERSITY OF ADELAIDE LIBRARIES

Article

Ab-initio Prediction of NMR Spin-Relaxation Parameters from Molecular Dynamics Simulations Po-chia Chen, Maggy Hologne, Olivier Walker, and Janosch Hennig J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00750 • Publication Date (Web): 02 Jan 2018 Downloaded from http://pubs.acs.org on January 3, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a 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.

Journal of Chemical Theory and Computation 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 19 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

Journal of Chemical Theory and Computation

Ab-initio Prediction of NMR Spin-Relaxation Parameters from Molecular Dynamics Simulations Po-chia Chen,∗,† Maggy Hologne,‡ Olivier Walker,‡ and Janosch Hennig† † Structural and Computational Biology Unit, EMBL Heidelberg, Meyerhofstrasse 1 69117 Heidelberg, Germany ‡ Université de Lyon, CNRS, Université Claude Bernard Lyon1, Ens de Lyon, Institut des Sciences Analytiques, UMR 5280, 5 rue de la Doua, F-69100 Villeurbanne, France E-mail: [email protected]

Abstract 1 H-15 N

NMR spin relaxation and relaxation dispersion experiments can reveal the time scale and extent of protein motions across the ps–ms range, where the ps–ns dynamics revealed by fundamental quantities R1 , R2 , and heteronuclear N OE can be well-sampled by molecular dynamics simulations (MD). Although the principles of relaxation prediction from simulation are well-established, numerous NMR–MD comparisons have hitherto focused upon the aspect of order parameters S 2 due to common artefacts in the prediction of transient dynamics. We therefore summarize here all necessary components and highlight existing and proposed solutions, such as the inclusion of quantum mechanical zero-point vibrational corrections, and separate MD convergence of global and local motions in coarse-grained and all-atom forcefields, respectively. To test the accuracy of MD prediction, two model proteins GB3 and Ubiquitin are used to validate five atomistic forcefields against published NMR data, supplemented by the coarsegrained forcefield MARTINI+EN. In Amber and CHARMM-type forcefields, quantitative agreement was achieved for structured elements with minimum adjustment of global parameters. Deviations from experiment occur in flexible loops and termini, indicating differences in both the extent and time-scale of backbone motions. Lack of systematic patterns and watermodel dependence suggest that modelling of the local environment limits prediction accuracy. Nevertheless, qualitative accuracy in a 2 µs-CHARMM36m Stam2 VHS-domain simulation demonstrates the potential of MD-based interpretation in combination with NMR-measured dynamics, increasing the utility of spin-relaxation in integrative structural biology.

Introduction

heteronuclear N OE, 3,4 and other parameters. 5,6 Whereas µs-ms dynamics are obtained by relaxation dispersion experiments. 7 These fundamental quantities characterizing heavy atom-proton relaxation rates reveal information about the timescales and extent of motions affecting the local backbone, which can be further integrated with molecular structure modelling 8–10 and dynamics simulations (MD) 11 to

Solution nuclear magnetic resonance spectroscopy (NMR) 1,2 is a powerful technique to observe biomolecular dynamics across psms timescales: Faster ps-ns dynamics are obtained by measuring 15 N spin-relaxation parameters such as longitudinal-relaxation R1 , transverse-relaxation R2 , 1 H-15 N steady-state

ACS Paragon Plus Environment

1

Journal of Chemical Theory and Computation 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

provide a more complete understanding of protein function, as influenced by temperature dependence, 12,13 intermolecular interactions, 14,15 protein disorder, 13,16–18 and enzymatic mechanisms. 7,19–22 The natural complementarity and access to similar timescales provides a direct method for cross-validation: of MD forcefields by spinrelaxation prediction, 11,23,24 and also of experimental model-free assumptions 25,26 by simulation. 27 Validation requires the translation of disparate experimental and simulation data onto a common basis for comparison. This ranges from the prediction of raw experimental measurements (i.e. R1 , R2 and N OE), 11,17,28,29 derivation of common analytical quantities such as order parameters S 2 , 30–35 to evaluating the heavy atom-proton vector autocorrelation functions C(t) and their Fourier-transformed spectral density J(ω), 36–38 the last of which constitutes the structural mechanism underlying spin relaxation. S 2 is a popular choice and directly obtainable from both NMR and MD data, and describes the degree of local order/disorder by measuring the rigidity of the N–H dipole via the relative contribution of local motions CI (t) versus global tumbling Drot . However, the popular Lipari-Szabo form of S 2 standard in NMR utilises a single time-constant and does not accurately represent potentially rich transient dynamics observed in simulation, while more extensive evaluations of S 2 require multi-field measurements to increase the descriptive power of fitted dynamics. The mapping of J(ω) from experiment can provide a more comprehensive picture when measurements at multiple magnetic fields are available, although is in practice often replaced by reduced-J(ω) variants 33,39,40 using a single field. Lastly, direct prediction of experimental measurements shifts all interpretive burden onto the simulation side, which can avoid mismatches between assumptions made in NMR versus MD analysis. This last option requires an estimation of quantum-mechanical phenomena such as chemical-shift anisotropy (CSA) and zero-point vibrations. MD prediction of spin-relaxation requires simulations of sufficient length to achieve sampling convergence, as well as corrections to re-

Page 2 of 19

move artefacts arising from, e.g., water modeldependence of Drot 29,39,41 and contribution of quantum-mechanical (QM) vibrational components. 42 Convergence by brute force is a simple matter of conducting ns–µs simulations to reliably compute C(t) and Drot . 41,43,44 On the other hand, water models are intimately tied to respective MD forcefields – alterations to water model behavior generally require extensive validation of its downstream effect upon dynamics and energetics, which hinders the popularisation of currently proposed solutions. 29,39 Until a concerted effort can be made, common practice is to substitute computed Drot with external knowledge either derived from experiment, 9,10 computed from static structures, 8 or simulated using coarse-grained MD. 44 This work aims to advance the state of affairs in the ab-initio prediction of NMR spinrelaxation and thus encourage broader utilisation of all the information it can offer, with similar goals as the recently-published Anderson et al. 45 whom utilize rotational velocity rescaling to correct Drot . We highlight three issues pertaining to quantitative prediction: First, the inclusion of QM-vibrational corrections (ζ) 42 sometimes ignored by MD studies, which introduce systematic baseline artefacts relative to experiment. Second, the replacement of Drot in all-atom MD to circumvent water-model artefacts and reduce convergence requirements. The coarse-grained MARTINI forcefield has been used here to supply independent and accurate Drot -predictions for folded proteins. 44 Third, the utilisation of an MD-based formalism that improves the modeling of Drot anisotropy over the standard modelfree formalism. These steps permit the computation of R1 , R2 , and N OE from two MD simulations alone, and exposes its methods to direct experimental scrutiny. The above approach is demonstrated exhaustively in two model proteins Ubiquitin (Ubq) and GB3, which are used to analyse the descriptions of local ps-ns dynamics according to five popular all-atom forcefields: Amber14sb, Amber99sb*-ILDN, CHARMM22*, CHARMM36m, and Gromos54A7. The accuracy obtained serves as a basis to discuss the


2

Page 3 of 19 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


tions 48 using the protein Stam2-VHS domain, whose dynamics contains a complete range between ordered and disordered dynamics. We hope that this contribution encourages the utilisation of NMR spin-relaxation in the repertoire of integrative structural biology collaborations alongside other MD–experiment combinations such as FRET 49,50 and small-angle scattering. 51–53

role of ζ and its relation to other global parameters such as 15 N chemical-shift anisotropy. 46,47 Spin-relaxation comparisons also expose inaccuracies with respect to the conformational fluctuations in semi-ordered motifs, which can be targeted in future forcefield development. To highlight the potential for future combined NMR and MD studies, we show a promising qualitative fidelity of CHARMM36m simula-

R1 =fDD [J(ωH − ωN ) + 3J(ωN ) + 6J(ωH + ωN )] + fCSA J(ωN ), R2 =0.5 ∗ fDD [4J(0) + J(ωH − ωN ) + 3J(ωN ) + 6J(ωH ) + 6J(ωH + ωN )] 1 + fCSA [4J(0) + 3J(ωN ] , 6 fDD γH N OE =1 + [6J(ωH + ωN ) − J(ωH − ωN )] , R1 γN 1 µ0 ~γN γH 2 fDD = ( ), 3 10 4πrN H 2 2 2 ∆CSA , fCSA = ωN 15

Theory and Methods

(1)

(2) (3) (4) (5)

measured in experiment can be classified into three distinct sources: global tumbling of the host domain CO (t), its local motions CI (t), as well as QM zero-point vibrations described by a constant correction factor ζ, 42 initially set as (1.02/1.041)6 ∼ 0.89. Although forcefield parametrisation do include modelling of the vN H out-of-plane and bending motions that affect C(t), classical mechanics cannot represent QM phenomena where, e.g., experimental C(t = 0) is less than unity. With the assumption that all three motions are uncorrelated, the relationship between J(ω) and C(t) becomes:

We summarise established theories of spinrelaxation 54 and model-free analysis 25,26 below with modifications relevant to MD predictions (c.f. d’Auvergne 55 §1.3 and §6.2 for a modern derivation). The rate of spin-relaxation in the backbone amide 15 N nucleus is primarily dependent on motions of its dipolar coupling to the bonded 1 H and the 15 N chemical shift anisotropy ∆CSA . The rotational motion of this NH-bond vN H determines an autocorrelation decay C(t), whose Fourier transform J(ω) evaluated at combinations of Larmor frequencies ωH and ωN underly the measured relaxation rates R1 , R2 and N OE as given in Eq. 1 through 5. 56,57 Here, ∆CSA is defined to be -170 ppm oriented along vN H , 46,58 and the vN H magnitude is defined to be 1.02 ˚ A according to neutron scattering 59 and NMR dipolar coupling. 60 Since vN H reorients over time, rather than selecting a single average vector we consider J(ω) to be an ensemble average over all vN H sampled during MD-trajectories. The vN H -motions

Z

J(ω) =

∞

dt cos(ωt)C(t) 0

C(t) = ζCO (t)CI (t).

,

(6)

vˆN H

(7)

Following the notation of Lee et al. 57 and for purposes of computational ease, we choose to absorb the 2/5 pre-factor in common J(ω) definitions 61 into fDD and fCSA in Eqs 4 and 5. This pre-factor arises from halving of the integral in Eq. 6 and the full solution of CO (t) in


3


Eq. 7 for spheroid tumbling. It is worth noting here that many proteins exhibit ns-timescale dynamics, and thus the nominal overlap of CO (t) and CI (t) technically violates the Lipari-Szabo separability criterion. Nevertheless, the appropriateness of this division has been repeatedly demonstrated for stable globular proteins, 28,32,41,62 and is thus maintained here. In this work, we approximate CO (t) with the global tumbling of an axially–symmetric rigid body where the rotational diffusion tensor Drot is comprised of axial components D⊥ and Dk , or equivalently isotropic tumbling Diso and axial anisotropy D∆ : 54,57 CO (t) =

3 X

Ai e−Di t ,

Page 4 of 19

vN H with the Drot frame: CLS (t) = S 2 + (1 − S 2 )e−kI t

(11)

Analytical transformation of single exponentials in C(t) thus yields the corresponding ensemble spectral density hJLS (ω)i, without QMcorrections: ∞

Z hJLS (ω)i =

dt cos(ωt)CO (t)CLS (t)

,

0

=

* 3 X

Ai

i=1

(1 − S 2 )(Di + kI ) S 2 Di + 2 2 (Di + kI )2 + ω 2 Di + ω

+ . (12)

In practice, the full set of internal motions influencing vN H in flexible residues cannot be wellapproximated by one or two components (Fig. S1). Thus, Eq. 11 must be expanded to capture sources across multiple timescales: 4,26,34,39

(8)

i=1

D1 D2 D3 A1

= 6D⊥ = 5D⊥ + Dk = 2D⊥ + 4Dk , = (3e2k − 1)2 /4

A2 = 3e2k (1 − e2k )

CI (t) = S 2 +

(9)

N X

αj e−kj t ,

(13)

j

where (10)

S2 +

A3 = 0.75(1 − e2k )2 .

N X

αj + αf = 1.

(14)

j

We include αf to capture timescales that occur faster than the trajectory sampling density (10 ps in this work). Eq. 13 is fitted to the raw computed CI (t) with increasing N from 1 ∼ 4 to identify the number of extractable components. Due to the fit merging contributions at similar timescales into composite timeconstants, we find that N is rarely greater than three. This fit is cut-off at τmem to exclude longer timescales that are hidden by Drot and thus unobservable in experiment, which is chosen as approximately twice the global tumbling time τiso = (6Diso )−1 . 33,34 The ensemble J(ω) after the above modifications is thus:

The axisymmetric form is taken instead of the fully-anisotropic form due to accuracy limitations in the brute-force computation of the rhombicity Dρ , 44 and its small value in many proteins. Values of Ai are determined based on the relative orientation of vN H within the coordinate frame of Drot , where ek is the axial component of vN H directed along the unique axis Dk . For GB3 and Ubq, the fullyanisotropic Drot is computed in a separate coarse-grained simulation using the MARTINI forcefield plus elastic networks, 63,64 by fitting simple-exponentials 41 to the auto-correlation of quaternion molecular orientations. 44 Values are then transferred to the all-atom representation with the same reference frame and an assumption that its conformational dependence is small. The Lipari-Szabo model–free formalism describes CI (t) by a single fast exponential with kI Drot and the complementary order parameter S 2 governing the residual alignment of

* hJ(ω)i =

ζ

3 X i=i

" Ai

X αj (Di + kj ) S 2 Di + (Di + kj )2 + ω 2 Di2 + ω 2

#+

j

(15)

The computation of Eq. 15 from MD is thus used to predict the spin-relaxation parameters in Eqs. 1 to 3.


4

Page 5 of 19

overall tumbling

local disorder

2.5

R1 R2 NOE

2

R1

R1

2

-1

rigid sphere isotropic model axisymm. model

1.5

rigid sphere isotropic model axisymm. model

1 8

10

1 0

7

5

R2

R2

6

6

4

5

3 0.8

4 0.8

NOE

-1

0.6

NOE

R1 / R2 [s ]

100

hetNOE [a.u.]

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 0.4

-3

0.1

1

10 τ [ns]

100

0.6 0.4

0

10

20

30 # Resid

40

50

0

10

20

30

40 50 # Resid

60

70

Figure 1: Primer on NMR spin-relaxation at proton frequency fH =600.133 MHz using theoretical and computed models. (left) Expected spin-relaxation values when vN H -reorientation is dominated by a single timescale τ , e.g. when embedded in a rigid-sphere undergoing tumbling, or associated with a disordered coil. (middle) Computed spin-relaxation for GB3 based on construction of J(ω) from coarse-grained global tumbling (Table S1) and all-atom CHARMM36m local dynamics. R1 , R2 and N OE have been plotted for three conditions: isotropic tumbling without or with local motion (black dashedline or orange circles), and axisymmetric tumbling with local motion (green diamonds). Secondary structures observed in simulations have been included for visual reference: coils (grey), α-helices (red), β-sheets (yellow arrows), β-bridges (tan-spikes), turns (cyan bulges), and 3-10 helices (pink). (right) Equivalent values for Ubq.

Results and Discussions

upon the timescale of slower motions – this can reveal ns-motions such as conformational rearrangements distinct from Diso . Further analysis and comparison enable the evaluation of additional influences such Rex denoting conformational exchange from slower µs-ms dynamics, and ligand-binding influences.

Basic characteristics of NMR spinrelaxation In order to contextualise the general readership, we illustrate an interpretation of spin-relaxation in Fig. 1 using CHARMM36m-simulated GB3 and Ubq. The dynamics affecting vN H –reorientation includes the baseline contribution from Diso that is equivalent to rigid-sphere tumbling, and additional contributions to local dynamics further modulated by the alignment of each dipole against Drot anisotropy (Fig. 1). We note that the inclusion of ζ exerts a systematic decrease of both R1 and R2 (compare Figs. S2 and S3). The qualitative lineshape reports upon differences in local dynamics: Less-flexible residues exhibit little deviation from the rigid sphere baseline other than via tumbling anisotropy. In contrast, more-flexible residues with significant ps-motions exhibit decreases in all values, in particular N OE. Drot anisotropy manifests where vN H is aligned with its long and short axes, noting that the dipoles in GB3 helices are aligned with Dk and those in GB3 sheets are aligned with D⊥ . After accounting for this, the relative R1 and R2 ratios inform

Structural-dependence of local dynamics predictions Having introduced the main characteristics of spin-relaxation, we now use experimental data to validate simulated backbone dynamics. In Figs. 2 and 3, the overall results are summarised using the performance of Amber14sb and CHARMM36m simulations of GB3 and Ubq dynamics respectively, and after χ-fitting of Diso and ζ to remove residual baseline errors. For the rest of this work we choose to report raw relaxations in the main text, while including fitted local dynamics of Eq. 13 in the Supporting Information. The complete dataset for five forcefields, with and without ζ-corrections and fitting are included in Figs. S2 through S7. The local dynamics parameters are included in Figs. S8 and S9. The two forcefields are representative of the


5


3

R1

R1

3

2

6

5

5

R2

6

4

4

3

3

1

1

0.5

0

10

20

30 # Resid

Expt. Charmm36m

1

NOE

R2 NOE

2

Expt. Amber14sb

1

40

50

0.5

60

0

10

20

30 # Resid

40

50

60

Figure 2: Computed relaxation parameters for GB3 at 600.133 MHz, based on the construction of J(ω) from coarsegrained global tumbling (Table S1) and all-atom local motions. (Left) The global rotational diffusion Drot extracted from 10,µs MARTINI simulations, superimposed upon its reference frame 1P7F.pdb. 65 (Center) NMR–relaxation computed via combining Drot with local motions from 500 ns Amber14sb simulations. Diso and ζ have been adjusted by χ Powellminimisation to remove slight baseline deviations. (Right) As above, but using 500 ns CHARMM36m simulations instead. The most frequently observed secondary-structures in respective simulations have been included for visual reference.

3

R1

2 Expt. Amber14sb

1 9 8 7 6 5 4 3

R2

R2

R1

3

2 Expt. Charmm36

1 9 8 7 6 5 4 3 1

NOE

1

NOE

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 19

0.5

0

10

20

30

40 50 # Resid

60

70

0.5

0

10

20

30

40 50 # Resid

60

70

Figure 3: Computed relaxation parameters for Ubq at 600.133 MHz, based on the construction of J(ω) from coarsegrained global tumbling (Table S1) and all-atom local motions. (Left) The global rotational diffusion Drot extracted from 10,µs MARTINI simulations, superimposed upon its reference frame 1UBQ.pdb. 66 (Center) NMR–relaxation computed via combining Drot with local motions from 500 ns Amber14sb simulations. Diso and ζ have been adjusted by χ Powellminimisation to remove slight baseline deviations. (Right) As above, but using 500 ns CHARMM36m simulations instead. The most frequently observed secondary-structures in respective simulations have been included for visual reference.


6

Page 7 of 19 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


performance of tested Amber and CHARMM families, able to quantitatively reproduce the spin-relaxation profiles in the majority of ordered regions. The matching of secondary structure features – e.g. the β-sheet saw-tooth pattern of GB3 E15-K19 and K50-T55, and the 3-turn α-helix pattern of Ubq T22-E34, – indicate that the combination of MARTINI global tumbling with all-atom local dynamics is sufficient to describe all dynamics observed in experiment, including perturbations via Drot anisotropy. Notably, this permits independent quantification of Rex contributions, e.g. analysed via R1ρ experiments for Ubq I25, N25, T55, and V70. 67 Forcefields can also partially model residues with significant backbone flexibility. Here, we include corresponding CI (t) for GB3 residues mentioned below as an interpretative aid in Fig. S10. For example, GB3-G41 undergoes fast and significant vibrations, while its structural analogue GB3-G14 is stabilised by internal hydrogen-bonds. However, tested forcefields underestimate the latter stabilisation and report significant drops in relaxation whether by itself or as part of an over-flexible β-turn. Our survey across all GB3 and Ubq β-turns reveal no systematic patterns in prediction accuracy: GB3-L12 turns are well described except for Amber14sb’s superfluous flexibility, but GB3-D47 turns are universally too rigid. On the other hand, experimental relaxations of Ubq-G10 turns cannot be simultaneously satisfied by forcefield predictions. This suggests that atomistic MD does not yet fully describe all factors influencing glycine motion and loop dynamics. Although forcefields do correctly identify structural effects such as stabilisation of coil-elements by adjacent secondary structures (c.f. GB3-A20 & V39, and Ubq-G35), existing parametrisation targets do not appear to cover transient dynamics that are reported by spinrelaxation. It is worth noting here that the extended Ubq L45-T55 stretch represents current frontiers in the modelling of semi-ordered behaviour, spanning multiple β-turns and bridges partially stabilising an otherwise coiled conformation. The significant qualitative differences between prediction and experiment indicates

the presence of conformational rearrangements on the ns-scale that are not observed in-silico. Therefore, we encourage caution in the interpretation of similar elements in future work where it impacts biological function. To summarise, MD simulations achieve semiquantitative success in replicating protein dynamics. Given that GB3 and Ubq are model proteins, this performance may deteriorate in more flexible systems. We have confirmed replicability of the Amber14sb results under small changes in temperature and water models (Fig. S11) – noting that the modified SPC/Eb water proposed by Takemura and Kitao 39 indeed mitigates the Drot artefacts otherwise produced by TIP3p water and additionally produces slightly improved β-turn dynamics. This highlights the influence of hydration models not just in global tumbling, but also in local dynamics. However, the slow convergence of D∆ ultimately limits brute-force Drot computation within all-atom MD. Overall quantification of force field accuracy Although forcefield performance is strongly influenced by local structure, for the proposes of general guidance in forcefield selection an overall comparison is desired. For this purpose, we obtained published GB3 relaxation data measured at five different proton frequencies ranging over 400–800 MHz, 68 used to exhaustively test experimental agreement with χ-agreement as the standard metric (Fig. 4). The performance of the Amber and CHARMMtype force fields are again nearly equivalent, noting the marginally better performance of Amber99sb*-ILDN and CHARMM22* variants that have been additionally tuned using long time-scale simulations. The poor performance of GROMOS–54A7 is due to increased psflexibility throughout the protein, leading to strong fluctuations in relaxation (Fig. S4). Although this forcefield was also tuned to reproduce certain NMR-observables of the protein backbone, 69 it appears that the 3 J-couplings and N OE-violations involved do not inform on N–H reorientation tested here. Nevertheless, this disparity of performance should instead be viewed as a reaping of benefits from the exten-


7


Amber14sb Amber99sb*-ILDN CHARMM22* CHARMM36m GROMOS 54A7

15

χ

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 8 of 19

10

5

0

R1

R2 NOE 400 MHz

R1

R2 NOE 500 MHz

R1

R2 NOE 600 MHz

R1

R2 NOE 700 MHz

R1

R2 NOE 800 MHz

Figure 4: σ–scaled χ agreement with experimental GB3 relaxation at 5 proton frequencies. Black bars indicate agreement where MARTINI-Diso and standard ζ were used, while colored bars indicate agreement where both values have been adjusted to minimise the aggregated R1/R2/NOE χ at a particular frequency. Note the relative height of χ (R1 ,500 MHz) is due to lower experimental uncertainties.

mura and Kitao 39 ), where systematic deviations of J(ω) remain even after authors’ corrections for Diso . In contrast, insertion of ζ is sufficient to remove almost all baseline artefacts in this work, such that further adjustments to global parameters are minuscule in comparison. It is worth noting that ζ is often included in studies reporting S 2 , where MD is known to otherwise overestimate experimental values. We note that an alternative to fitting ζ is to increase rN H from 1.02 to ∼1.04 ˚ A reflecting the average effect of vibrations. 42 While this is adopted by e.g. Anderson et al. 45 , we believe it conflicts with the philosophical separation of structure and dynamics in MD forcefields. A second alternative to ζ is to modify average ∆CSA , and a similar fitting was conducted with χ-values for all conditions reported in Table 1. A slightly-worse χ-agreement could be achieved by lowering ∆CSA to ∼-135 ppm for both GB3 and Ubq. This significantly contradicts consensus values previously reported, 46,68,73 and should be rejected. Thus ζ-corrections appear to be the principal factor underlying previous overestimation of spin-relaxation parameters, with residual variations stemming from forcefield inaccuracies and the influence of the local environment on ζ and ∆CSA .

sive and ongoing parametrisation work of both Amber and CHARMM communities, 11,23 particularly on the focused refinement of backbone dihedrals. 70,71 It is worth noting here explicit improvements by Amber99sb for Ubq 11 and GB3. 31 We also note that the residual differences of the best forcefields appear to approach limits imposed by the modelling of ∆CSA as a constant rather than incorporating site-specific variations. According to the 5-field GB3 measurements 68 in which a range of ∆CSA models was employed, the average deviation of experimental S 2 is ∼0.025 and comparable to observed deviations in our NOE measurements (Fig. S4 and S7). Influence of global parameters on baseline deviations We return to the question of the necessity of ζ in quantitative spin-relaxation prediction, and compare χ-agreement under a number of assumptions in Table 1. As mentioned above in passing, the absence of QMcorrections leads to consistent over-prediction of R1 and R2 in GB3 and Ubq (Figs. S2 and S5). Similar artefacts were discovered in published literature, including the yeast transcription factor GCN4 (Fig. S5 of Robustelli et al. 33 ) and an independent Ubq/GB3 study (Fig. 5 of Take-


8

Page 9 of 19 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


Table 1: Influence of parameter optimisation on R1 , R2 , and N OE-agreement. GB3 values including uncertainties are averages from separate minimisation against fH =400,600,700, and 800 MHz. GB3 Diso standardised to 297 K and 9%D2 O. Experimental data for GB3 72 and Ubq have been generously provided by respective authors. Default ζ is adopted from (1.02/1.041)6 . 42 Dashes (–) indicate use of fixed default values. † : Values from 10 µs-MARTINI simulations.

Amber14sb

CHARMM36m GB3

Parameters fitted

Reference

None

Diso

Diso ζ

Diso ∆CSA

None

Diso

Diso ζ

Diso ∆CSA

Diso [107 s−1 ]

4.99±0.16

4.93±0.23†

4.87±0.04

4.97±0.05

5.02±0.10

4.93±0.23†

4.94±0.05

5.00±0.05

5.03±0.07

∆CSA [ppm]

-170

–

–

–

-183±5

–

–

–

-177±3

ζ

0.890

–

–

0.918±0.006

–

–

–

0.903±0.004

–

χ

n.a.

5.88±0.65

5.84±0.64

5.34±0.51

5.21±0.81

5.19±0.83

5.00±0.79

Parameters fitted

Reference

None

Diso

5.31±0.49 Diso ∆CSA

None

Diso

4.99±0.79 Diso ∆CSA

Diso [107 s−1 ]

4.99±0.16

4.93±0.23†

5.17±0.20

4.86±0.09

4.93±0.23†

5.25±0.20

4.89±0.11

∆CSA [ppm]

-170

–

–

-135±25

–

–

-127±31

ζ

1.00

–

–

–

–

–

–

χ

n.a.

7.99±0.35

7.57±0.56

5.55±0.37

8.98±0.57

8.37±1.00

5.33±0.67

None

Diso

Diso ζ

Parameters fitted

Reference

None

Diso

Ubiquitin Diso Diso ζ ∆CSA

Diso [107 s−1 ]

3.84±0.01

3.74†

3.97

3.93

3.96

3.74†

4.01

3.97

4.01

∆CSA [ppm]

-170

–

–

–

-171

–

–

–

-170

ζ

0.890

–

–

0.897

–

–

–

0.898

–

χ

n.a.

11.8

10.4

10.4

14.2

12.9

Reference

None

Diso

None

Diso

12.7 Diso ∆CSA

12.9

Parameters fitted

10.1 Diso ∆CSA

Diso [107 s−1 ]

3.84±0.01

3.74†

3.64

3.85

3.74†

3.69

3.88

∆CSA [ppm]

-170

–

–

-135

–

–

-135

ζ

1.00

–

–

–

–

–

–

χ

n.a.

19.0

18.8

8.04

19.4

19.4

10.7

Diso ∆CSA

Ubq-binding during lysosomal degradation, 74 and exhibits here in the apo-state slow conformational fluctuations on the order ∼ 10 ns that suggests conformational-selection as its binding mechanism. The dynamics shown by MD is qualitatively if not quantitatively in agreement with experiment, partially confirming the observations. A surprisingly gradual transition from order to disorder is also observed at the C-terminus from Drot -dominated dynamics at E146 to partial-disorder at Q156. The MD interpretation of this transition attributes this to a stabilising effect of F150-binding into a hydrophobic cleft on the helix-bundle. The lower flexibility of double prolines P151/P152 also reduces the reorientation rate of downstream residues. On the other hand, simulations do not predict a similar partial-ordering effect at the

Dynamics of mixed ordered-disordered systems Although spin-relaxation experiments can also be applied to study protein disorder, MD approaches currently lack sufficiently well-characterised forcefields in either coarse-grained or all-atom contexts. The recently published CHARMM36m is a significant improvement on this front, and to illustrate its potential application we demonstrate a 2 µssimulation of Stam2-VHS in CHARMM36m alone, which is a helix-bundle domain containing disordered termini. After removing baseline artefacts, experimental and computed spin-relaxation reveal nsscale dynamics characterised by increased R1 and decreased R2 (Fig. 5 and Fig. S12) relative to the slower global tumbling. The loopregion spanning E26-D32 is associated with


9


Page 10 of 19

Figure 5: Depiction of VHS and relaxation predictions. (left) All-atom based Drot from 2 µs CHARMM36m simulations, superimposed in 1X5B.pdb with the C-terminal is oriented upwards. Protein illustrated in yellow cartoon and surface representations with a number of coil residues depicted as sticks. (right) Experimental and simulation relaxation profiles, sorted from top to bottom as R1 , R2 , N OE, and component-wise breakdown of fitted CI (t) into S 2 and αj of eq. 13 sorted according to timescales kj . Secondary structures have been included above for visual context. In the raw relaxation comparison, simulated data was colored green to indicate that its experimental counterpart (black) could be observed, and cyan otherwise. Drot and ζ was adjusted to minimise baseline discrepancies.

N-terminus, where an analogous gradual transition is found by experiment. This suggests that the exposed hydrophobic residues M8 and F11 may exhibit similar stabilising interactions that were not discovered by MD sampling, and not present in initial NMR-derived structures. To summarise the demonstration, the accuracy of CHARMM36m appears to be sufficient now to produce meaningful interpretations of spin-relaxation dynamics, although still limited by feasibility of required convergence. The uncertainty of D∆ in particular limits quantitative accuracy for structured regions, although this hindrance can be mitigated in comparative studies.

of this common basis shifts the burden of translation between the two extremes of structurallytransparent parameters such as S 2 and C(t) versus parameters directly accessible to experiment. S 2 in particular is intuitive and easily accessible from both NMR and MD. However, the 1 − S 2 components in the NMR-standard Lipari-Szabo formalism assumes the dominance of one or two components: while this is suitable for rigid structural elements, it often fails to describe the rich transient dynamics characteristic of flexible residues under the influence of the complex local environment. Rather than applying a simplifying model to the complex dynamics described by MD simulations, it would be more informative to conduct sufficient independent NMR measurements so as to construct improved dynamics models. Here, more refined spin-relaxation quantities can be utilised, such as the exchange-free transverse relaxation Rdd 6 and the 15 N CSA/15 N-1 H dipolar longitudinal ηz and transverse ηxy cross-correlation, 5 which are also able to disambiguate influences of Rex and ∆CSA -variations that are not cur-

Summary The direct comparison between MD- and NMRmeasured dynamics involves conversion from primary data – trajectories in MD and relaxation in NMR, – into some common basis on which interpretations can be made. The choice


10

Page 11 of 19 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


rently modelled within MD. In lieu of extensive efforts on the experimental front, the ubiquitous R1 , R2 , and heteronuclear N OE form a valid alternative basis of comparison that does not enforce a comparison between model-free and model-based (MD) formalisms. Here, simulation is capable of sufficiently sampling ps-ns motions, Drot anisotropy, and vN H orientations that are together necessary to achieve quantitatative comparisons with experiment and expose forcefield limitations, Rex contributions, and ∆CSA variations. In the rigid, structural elements of biomolecules, good agreements can be readily obtained by utilising Amber- and CHARMM-type forcefields and minimal adjustment of global parameters. This demonstrates the maturity of their development with respect to folded protein modelling. In regions of significant flexibility, discrepancies reveal the influence of both forcefield artefacts and local environment upon the dynamics of various β-turns and coils connecting structured elements. In particular, the slower ns-scale motions characteristic of conformational fluctuations are not consistently reproduced. Given that the GB3 and Ubq tested here are considered to be model proteins with supreme stability, a more in-depth investigation into other systems will be required to fully investigate forcefield quality with respect to semiordered behaviour. Here, spin-relaxation is well-suited to furnish information on timescale and extent of motions governing semi-ordered, and disordered behaviour for ongoing forcefield development. To further illustrate this, recent improvements in CHARMM36m with respect to disordered peptides was applied to simulate VHS, demonstrating that MD can qualitatively capture all dynamics between order and disorder. One should also note that concommitant modifications to water-models 75 may be necessary to further improve accuracy of disordered protein models. We also devoted significant space for ζ corrections to highlight its varied and inconsistent adoption across literature. Strictly speaking, the zero-point vibrations experienced by vN H implies CI (0) < 1 in experiment, and thus both the Lipari-Szabo model-free formalism and clas-

sical dynamics in MD are incomplete. While the QM potential energy surface is a core part of forcefield parametrisation, the quantum versus classical dynamics occurring upon this surface are fundamentally different. This subtely has likely contributed to persistence of baseline artefacts in literature, notwithstanding cancellation of errors via e.g. insufficient sampling. We also note that ζ-corrections are inevitably entangled with the ps-dynamics described by a given MD setup. Under-stabilization of the Gromos54A7 protein backbone leads to artificially large ps-scale motions, which the fitting process compensates for via a 5∼10% increase of ζ relative to other forcefields (Table S2). Given this, variations in ζ may occur due to protein structure, forcefield, and even simulation parameters – just as ∆CSA is dependent upon the local shielding environment. We note here that ±20 ppm variations of ∆CSA exert a ∼ 5% effect on relaxation values for rigid segments (Fig. 4 of Fushman and Cowburn 76 ). For flexible residues, both ζ and ∆CSA may exhibit larger deviations from existing averages modelled using stable proteins. Therefore, a focussed study using well-known systems containing flexibility will be required to investigate how to further improve relaxation prediction and disentangle transient dynamics from the above parameter variations. In summary, we have refined a combined coarse-grained and all-atom MD approach to predict ab initio experimental NMR spinrelaxation at the ps-ns level. By explicitly deriving such predictions from dynamics alone, this tool is able to independently interpret experimentally-observed dynamics with sampled atomistic mechanisms, generating a valuable integrated model useful to guide further research. The use of coarse-grained axisymmetric global tumbling effectively mitigates resource-limitations upon simulation convergence as well as all-atom water-model artefacts, permitting full use of atomistic simulations at feasible timescales. With the addition of ζ and fitting, these advances are sufficient to quantify experimental tumbling measurements. As a glimpse of the future potential of in silico predictions using all-atom simulation data alone,


11


we demonstrated the case of V HS where MD furnishes reasonable predictions and explanations for significant ns-level motions in semiordered residues. Thus, with further increases in affordable computational power, we hope to see full utilisation of spin-relaxation in combined experimental and computational studies of biomolecular structural dynamics.

script Drawer).

Page 12 of 19

(https://github.com/zharmad/SSS-

Simulation Setup Source coordinates for proteins have been obtained from the Protein Databank for the GB3 domain (1P7F), 65 Ubq (1XQQ), 84 and Stam2 VHS-domain (1X5B). Protein-specific values below are reported respectively in this order. The size of dodecahedron simulation boxes have been determined by applying 1.5 nm buffer to the protein coordinates. This produces system sizes of ∼20,000, ∼22,000, and ∼115,000 atoms, which are reduced to 563, 1774, and 9443 atoms after application of MARTINI coarse-graining. The MARTINI-EN VHS system is not used in this work due to its inability to model protein disorder. MARTINI equilibration involves 2500-step energy minimisation of the coarse-grained coordinates in vacuum, a second energy minimisation after solvation in 150 mM ions, 100 ps stepwise thermalization to 300 K in 200 ps NVT simulations, and finally 200 ps-NPT simulations. All-atom equilibration of GB3 and Ubq involves solvation in 150 mM NaCl, removal of isolated intruding water molecules, 2500-step energy minimisation, followed by thermalization and equilibration in 200 ps NPT simulations with 1000 kJ mol−1 nm−2 restraints on all protein atoms. In the case of VHS where more careful equilibration was justified, a gradual relaxation of side-chain then backbone restraints over 2 ns was conducted from 5000 to 0 kJ mol−1 nm−2 . Other relevant parameters such as temperature and pressure controls are provided in Table 3, noting that separate thermostats were applied to protein and solvent atoms. Simulation trajectories have been recorded every 10 ps for analysis. Since convergence of autocorrelation dynamics by bruteforce sampling requires simulation lengths of the order 100×τ , 43 we conducted 500 ns allatom simulations for GB3 and Ubq, but 2 µs simulations for VHS. Whereas 10 µs-MARTINI simulations as conducted for GB3 and UBq.

Materials and Methods NMR Data Collection Primary measurements have been taken from published work for GB3 68,72 and VHS. 14 Ubq measurements have been generously provided (see Acknowledgements). Computational Clusters We acknowledge the following computational clusters for providing the necessary resources to compute µslevel simulations: the CiNES cluster in Montpellier, the EMBL cluster in Heidelberg, and the GWDG cluster in Göttingen. Software, algorithm and forcefield distributions All MD simulations in this study have been conducted using GROMACS 5.1.x 77 software. The following forcefields are either available directly or downloadable from the GROMACS website: Amber14sb, 78 Amber99sb*-IDLN, 79 CHARMM22*, 80 and Gromos54A7. 69 MARTINI (v2.2) 63,64 and CHARMM36m (November 2016) 48 are available from respective developer websites in GROMACS-compatible form. MARTINI forcefield is used with elastic networks to restrain conformations to near-native states. Forcefieldrelated parameters, including cut-off schemes and timestepping, are shown in Table 2. All bond lengths have been constrained via the LINCS algorithm. Quaternion-based computation of Drot were conducted using a forked distribution 44 of PLUMED-2. 81 Protein visualisations were produced in VMD. 82 Secondary structures are computed according to STRIDE 83 definitions, averaged across conformations observed in respective all-atom simulations, and visualized using a simple python

Computation of CO (t) and CI (t) For the purposes of defining atoms associated


12

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


with the reference frame, we include all atoms of GB3 and Ubq, and residues 15 to 146 of VHS. The use of all-atoms of ordered domains instead of defining rigid structured elements minimizes human interpretation and forcefield-dependence to decisions based on on residue-level disorder. This produces negligible difference of The computation of Drot using our PLUMED-fork 44 (https://github.com/zharmad/plumed2) is based on the autocorrelation of the molecular orientation expressed as a quaternion Q, whose vector components contain information on the axial elements of Drot after appropriate frame transformations. The numerical resolution of CI (t) and CO (t) was set to every 10 ps, which implicitly assumes that fast motions < 10 ps in CI (t) was effectively collapsed into αf not visible in J(ω). Based on the predicted isotropic tumbling of GB3 and Ubq at 3∼4 ns, τmem is chosen to be 10 ns at which both CO (t) and CI (t) are cut-off. We note here that a range of τrot 33 ∼ τmem ∼ 5τrot 34 exists in literature, although increasing τmem to 20 ns results only in minor changes to resultant relaxation. Our current implementation of spin-relaxation computation workflows will be made publicly available (https://github.com/zharmad/SpinRelax), based on Python2.7 and utilising MDTraj 85 to read trajectory data. Code-refactoring for MDTraj inclusion is being considered.

access to the GWDG Cluster.

Supporting Information The supporting information document contains additional results from coarse-grained simulations, full five-forcefield comparisons, and fitted parameters. Table S1 records the computed Drot used in spin-relaxation prediction. Table S2 records adjusted Diso and ζ values to match the experimental baseline relaxation, for each all-atom forcefield prediction. Figures S1 to S3 records the five-forcefield spinrelaxation predictions of GB3, under three different assumptions of global parameter values. Figures S4 to S6 records the five-forcefield spinrelaxation predictions of Ubq, under three different assumptions of global parameter values. Figure S7 shows the simulated and fitted CI (t) for a small subset of GB3 residues for the five forcefields as an example of ordered and semiordered dynamics. Figure S8 collates additional results when Amber14sb is applied with the SPC/Eb water model. This information is available free of charge via the Internet at http://pubs.acs.org

References (1) Ishima, R.; Torchia, D. A. Protein Dynamics from NMR. Nat. Struct. Mol. Biol. 2000, 7, 740–743. (2) Kovermann, M.; Rogne, P.; Wolf-Watz, M. Protein Dynamics and Function from Solution State NMR Spectroscopy. Q. Rev. Biophys. 2016, 49 .

Acknowledgements We are grateful to Prof D. Fushman for providing the relaxation data for Ubq and GB3 at different spectrometer fields. The Centre Informatique National de l’Enseignement Supérieur (CINES) and EMBL HPC cluster are also acknowledged for providing computational ressources. OW and MH are supported by the French National Research Agency as part of the ”OH risque” programme (Metadyn, ANR-14OHRI-0006-01). PC is supported by EMBL and a EU Marie Curie Actions Cofund grant. We also thank Dr. Bernd Simon for valuable discussions on relaxation theory and its computation, and Dr. Jochen Hub for providing

(3) Br¨ uschweiler, R. New Approaches to the Dynamic Interpretation and Prediction of NMR Relaxation Data from Proteins. Curr. Opin. Struct. Biol. 2003, 13, 175–183. (4) Morin, S. A Practical Guide to Protein Dynamics from 15N Spin Relaxation in Solution. Prog. Nucl. Magn. Reson. Spectrosc. 2011, 59, 245–262. (5) Kroenke, C. D.; Rance, M.; Palmer, A. G. Variability of the 15N Chemical Shift Anisotropy in Escherichia Coli Ribonuclease H in Solution. J. Am. Chem. Soc. 1999, 121, 10119–10125. (6) Hansen, D. F.; Yang, D.; Feng, H.; Zhou, Z.; Wiesner, S.; Bai, Y.; Kay, L. E. An ExchangeFree Measure of 15N Transverse Relaxation: An


13


Page 14 of 19

Table 2: Forcefield-related parameters. Twin-range switching functions are indicated by double values. Forcefield Water-model Timestep [fs] Cut-off Scheme Pairlist [nm] Electrostatics Short-range cut-off [nm] Van-der-Waals Cut-offs [nm] Dispersion Corr.

MARTINI CG 20 Group 1.4 Shift 1.2 Shift 0.9/1.2 None

Amber14sb TIP3p

1.0 Potential-shift 1.0

Amber99sb*-ILDN CHARMM22* CHARMM36m TIP3p or SPC/Eb 39 TIP3p TIP3p 4 with virtual sites 86 Verlet n.a. Particle-Mesh Ewald 87 1.0 1.2 1.2 Potential-shift Force-Switch Force-Switch 1.0 1.0/1.2 1.0/1.2 Energy+Pressure

Gromos54A7 TIP3p

1.2 Potential-shift 1.2

Table 3: Simulation-related parameters. Formalism Step Thermostat τt [ps] Tref [K] Barostat τp [ps] Pref [atm]

coarse-grained Equilibration Production V-rescale 88 1.0 1.0 300 300 Parrinello-Rahman 90 12.0 5.0 1.0 1.0

Equilibration V-rescale 1.0 300 Berendsen 91 2.5 1.0

all-atom Production Nose-Hoover 89 2.5 300 Parrinello-Rahman 5.0 1.0

NMR Spectroscopy Application to the Study of a Folding Intermediate with Pervasive Chemical Exchange. J. Am. Chem. Soc. 2007, 129, 11468– 11479.

Griesinger, C.; Zweckstetter, M.; Luchinat, C. Long-Range Correlated Dynamics in Intrinsically Disordered Proteins. J. Am. Chem. Soc. 2014, 136, 16201–16209.

(7) Palmer III, A. G. Chemical Exchange in Biomacromolecules: Past, Present, and Future. J. Magn. Reson. 2014, 241, 3–17.

(14) Lange, A.; Hoeller, D.; Wienk, H.; Marcillat, O.; Lancelin, J.-M.; Walker, O. NMR Reveals a Different Mode of Binding of the Stam2 VHS Domain to Ubiquitin and Diubiquitin,. Biochemistry 2011, 50, 48–62.

(8) Garcıá de la Torre, J.; Huertas, M. L.; Carrasco, B. HYDRONMR: Prediction of NMR Relaxation of Globular Proteins from Atomic-Level Structures and Hydrodynamic Calculations. J. Magn. Reson. 2000, 147, 138–146.

(15) Rezaei-Ghaleh, N.; Klama, F.; Munari, F.; Zweckstetter, M. Predicting the Rotational Tumbling of Dynamic Multidomain Proteins and Supramolecular Complexes. Angew. Chem. Int. Ed. 2013, 52, 11410–11414.

(9) Dosset, P.; Hus, J.-C.; Blackledge, M.; Marion, D. Efficient Analysis of Macromolecular Rotational Diffusion from Heteronuclear Relaxation Data. J. Biomol. NMR 2000, 16, 23–28.

(16) Cino, E. A.; Karttunen, M.; Choy, W.-Y. Effects of Molecular Crowding on the Dynamics of Intrinsically Disordered Proteins. PLoS ONE 2012, 7, e49876.

(10) Walker, O.; Varadan, R.; Fushman, D. Efficient and Accurate Determination of the Overall Rotational Diffusion Tensor of a Molecule from 15N Relaxation Data Using Computer Program ROTDIF. J. Magn. Reson. 2004, 168, 336–345.

(17) Salvi, N.; Abyzov, A.; Blackledge, M. MultiTimescale Dynamics in Intrinsically Disordered Proteins from NMR Relaxation and Molecular Simulation. J. Phys. Chem. Lett. 2016, 7, 2483– 2489.

(11) Showalter, S. A.; Br¨ uschweiler, R. Validation of Molecular Dynamics Simulations of Biomolecules Using NMR Spin Relaxation as Benchmarks: Application to the AMBER99SB Force Field. J. Chem. Theory Comput. 2007, 3, 961–975.

(18) Mushtaq, A. U.; Park, J. S.; Bae, S.-H.; Kim, H.Y.; Yeo, K. J.; Hwang, E.; Lee, K. Y.; Jee, J.G.; Cheong, H.-K.; Jeon, Y. H. Ligand-Mediated Folding of the OmpA Periplasmic Domain from Acinetobacter Baumannii. Biophys. J. 2017, 112, 2089–2098.

(12) Chang, S.-L.; Tjandra, N. Temperature Dependence of Protein Backbone Motion from Carbonyl 13C and Amide 15N NMR Relaxation. J. Magn. Reson. 2005, 174, 43–53.

(19) Eisenmesser, E. Z.; Bosco, D. A.; Akke, M.; Kern, D. Enzyme Dynamics During Catalysis. Science 2002, 295, 1520–1523.

(13) Parigi, G.; Rezaei-Ghaleh, N.; Giachetti, A.; Becker, S.; Fernandez, C.; Blackledge, M.;


14

Page 15 of 19 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) Grey, M. J.; Wang, C.; Palmer, A. G. Disulfide Bond Isomerization in Basic Pancreatic Trypsin Inhibitor: Multisite Chemical Exchange Quantified by CPMG Relaxation Dispersion and Chemical Shift Modeling. J. Am. Chem. Soc. 2003, 125, 14324–14335.

(30) Palmer, A. G.; Case, D. A. Molecular Dynamics Analysis of NMR Relaxation in a Zinc-Finger Peptide. J. Am. Chem. Soc. 1992, 114, 9059–9067. (31) 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.

(21) Hass, M. A. S.; Thuesen, M. H.; Christensen, H. E. M.; Led, J. J. Characterization of M s-ms Dynamics of Proteins Using a Combined Analysis of 15N NMR Relaxation and Chemical Shift: Conformational Exchange in Plastocyanin Induced by Histidine Protonations. J. Am. Chem. Soc. 2004, 126, 753–765.

(32) Maragakis, P.; Lindorff-Larsen, K.; Eastwood, M. P.; Dror, R. O.; Klepeis, J. L.; Arkin, I. T.; Jensen, M. Ø.; Xu, H.; Trbovic, N.; Friesner, R. A. et al. Microsecond Molecular Dynamics Simulation Shows Effect of Slow Loop Dynamics on Backbone Amide Order Parameters of Proteins. J. Phys. Chem. B 2008, 112, 6155–6158.

(22) Vallurupalli, P.; Chakrabarti, N.; Pomès, R.; Kay, L. E. Atomistic Picture of Conformational Exchange in a T4 Lysozyme Cavity Mutant: An Experiment-Guided Molecular Dynamics Study. Chem. Sci. 2016, 7, 3602–3613.

(33) Robustelli, P.; Trbovic, N.; Friesner, R. A.; Palmer, A. G. Conformational Dynamics of the Partially Disordered Yeast Transcription Factor GCN4. J. Chem. Theory Comput. 2013, 9, 5190– 5200.

(23) Huang, J.; MacKerell, A. D. CHARMM36 AllAtom Additive Protein Force Field: Validation Based on Comparison to NMR Data. J. Comput. Chem. 2013, 34, 2135–2145.

(34) Gu, Y.; Li, D.-W.; Br¨ uschweiler, R. NMR Order Parameter Determination from Long Molecular Dynamics Trajectories for Objective Comparison with Experiment. J. Chem. Theory Comput. 2014, 10, 2599–2607.

(24) Chalmers, G.; Glushka, J. N.; Foley, B. L.; Woods, R. J.; Prestegard, J. H. Direct NOE Simulation from Long MD Trajectories. J. Magn. Reson. 2016, 265, 1–9.

(35) Genheden, S. Effect of Solvent Model When Probing Protein Dynamics with Molecular Dynamics. J. Mol. Graph. Model. 2017, 71, 80–87.

(25) Lipari, G.; Szabo, A. Model-Free Approach to the Interpretation of Nuclear Magnetic Resonance Relaxation in Macromolecules. 1. Theory and Range of Validity. J. Am. Chem. Soc. 1982, 104, 4546– 4559.

(36) Farrow, N. A.; Zhang, O.; Szabo, A.; Torchia, D. A.; Kay, L. E. Spectral Density Function Mapping Using 15N Relaxation Data Exclusively. J Biomol NMR 1995, 6, 153–162.

(26) Clore, G. M.; Szabo, A.; Bax, A.; Kay, L. E.; Driscoll, P. C.; Gronenborn, A. M. Deviations from the Simple Two-Parameter Model-Free Approach to the Interpretation of Nitrogen-15 Nuclear Magnetic Relaxation of Proteins. J. Am. Chem. Soc. 1990, 112, 4989–4991.

(37) Atkinson, R. A.; Lefèvre, J.-F. Reduced Spectral Density Mapping for Proteins: Validity for Studies of 13C Relaxation. J Biomol NMR 1999, 13, 83– 88.

(27) Peter, C.; Daura, X.; van Gunsteren, W. F. Calculation of NMR-Relaxation Parameters for Flexible Molecules from Molecular Dynamics Simulations. J. Biomol. NMR 2001, 20, 297–310.

(38) Kadeˇrávek, P.; Zapletal, V.; Rabatinová, A.; ˇ ıdek, L. Spectral DenKrásn´ y, L.; Sklenáˇr, V.; Z´ sity Mapping Protocols for Analysis of Molecular Motions in Disordered Proteins. J Biomol NMR 2014, 58, 193–207.

(28) Prompers, J. J.; Br¨ uschweiler, R. General Framework for Studying the Dynamics of Folded and Nonfolded Proteins by NMR Relaxation Spectroscopy and MD Simulation. J. Am. Chem. Soc. 2002, 124, 4522–4534.

(39) Takemura, K.; Kitao, A. Water Model Tuning for Improved Reproduction of Rotational Diffusion and NMR Spectral Density. J. Phys. Chem. B 2012, 116, 6279–6287.

(29) Anderson, J. S.; LeMaster, D. M. Rotational Velocity Rescaling of Molecular Dynamics Trajectories for Direct Prediction of Protein NMR Relaxation. Biophys. Chem. 2012, 168–169, 28–39.

(40) Lindorff-Larsen, K.; Trbovic, N.; Maragakis, P.; Piana, S.; Shaw, D. E. Structure and Dynamics of an Unfolded Protein Examined by Molecular Dynamics Simulation. J. Am. Chem. Soc. 2012, 134, 3787–3791.


15


(41) Wong, V.; Case, D. A. Evaluating Rotational Diffusion from Protein MD Simulations. J. Phys. Chem. B 2008, 112, 6013–6024.

Page 16 of 19

(53) Björling, A.; Niebling, S.; Marcellini, M.; van der Spoel, D.; Westenhoff, S. Deciphering Solution Scattering Data with Experimentally Guided Molecular Dynamics Simulations. J. Chem. Theory Comput. 2015, 11, 780–787.

(42) Case, D. A. Calculations of NMR Dipolar Coupling Strengths in Model Peptides. J. Biomol. NMR 1999, 15, 95–102.

(54) Woessner, D. E. Nuclear Spin Relaxation in Ellipsoids Undergoing Rotational Brownian Motion. J. Chem. Phys. 1962, 37, 647–654.

(43) Lu, C.-Y.; Bout, D. A. V. Effect of Finite Trajectory Length on the Correlation Function Analysis of Single Molecule Data. J. Chem. Phys. 2006, 125, 124701.

(55) d’Auvergne, E. J. Protein Dynamics: A Study of the Model-Free Analysis of NMR Relaxation Data. Ph.D. thesis, University of Melbourne, 2006.

(44) Chen, P.-c.; Hologne, M.; Walker, O. Computing the Rotational Diffusion of Biomolecules via Molecular Dynamics Simulation and Quaternion Orientations. J. Phys. Chem. B 2017,

(56) Abragam, A. The Principles of Nuclear Magnetism; Oxford University Press, 1961. (57) Lee, L. K.; Rance, M.; Chazin, W. J.; Palmer, A. G. Rotational Diffusion Anisotropy of Proteins from Simultaneous Analysis of 15N and 13Cα Nuclear Spin Relaxation. J Biomol NMR 1997, 9, 287–298.

(45) Anderson, J. S.; Hernández, G.; LeMaster, D. M. Prediction of Bond Vector Autocorrelation Functions from Larmor Frequency-Selective Order Parameter Analysis of NMR Relaxation Data. J. Chem. Theory Comput. 2017, 13, 3276–3289.

(58) Lee, A. L.; Wand, A. J. Assessing Potential Bias in the Determination of Rotational Correlation Times of Proteins by NMR Relaxation. J Biomol NMR 1999, 13, 101–112.

(46) Tjandra, N.; Szabo, A.; Bax, A. Protein Backbone Dynamics and 15N Chemical Shift Anisotropy from Quantitative Measurement of Relaxation Interference Effects. J. Am. Chem. Soc. 1996, 118, 6986–6991.

(59) Kvick, ˚ A.; Al-Karaghouli, A. R.; Koetzle, T. F. Deformation Electron Density of α-Glycylglycine at 82 K. I. The Neutron Diffraction Study. Acta Crystallogr. B 1977, 33, 3796–3801.

(47) Fushman, D.; Tjandra, N.; Cowburn, D. Direct Measurement of 15N Chemical Shift Anisotropy in Solution. J. Am. Chem. Soc. 1998, 120, 10947– 10952.

(60) Yao, L.; Vögeli, B.; Ying, J.; Bax, A. NMR Determination of Amide N-H Equilibrium Bond Length from Concerted Dipolar Coupling Measurements. J. Am. Chem. Soc. 2008, 130, 16518–16520.

(48) Huang, J.; Rauscher, S.; Nawrocki, G.; Ran, T.; Feig, M.; de Groot, B. L.; Grubm¨ uller, H.; MacKerell Jr, A. D. CHARMM36m: An Improved Force Field for Folded and Intrinsically Disordered Proteins. Nat. Methods 2017, 14, 71–73.

(61) Ghose, R.; Fushman, D.; Cowburn, D. Determination of the Rotational Diffusion Tensor of Macromolecules in Solution from NMR Relaxation Data with a Combination of Exact and Approximate Methods—Application to the Determination of Interdomain Orientation in Multidomain Proteins. J. Magn. Reson. 2001, 149, 204–217.

(49) Kalinin, S.; Peulen, T.; Sindbert, S.; Rothwell, P. J.; Berger, S.; Restle, T.; Goody, R. S.; Gohlke, H.; Seidel, C. A. M. A Toolkit and Benchmark Study for FRET-Restrained High-Precision Structural Modeling. Nat. Methods 2012, 9, 1218– 1225.

(62) Johnson, E. Separability between Overall and Internal Motion: A Protein Folding Problem. Proteins 2012, 80, 2645–2651.

(50) Hoefling, M.; Grubm¨ uller, H. In Silico FRET from Simulated Dye Dynamics. Comp. Phys. Commun. 2013, 184, 841–852.

(63) Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; de Vries, A. H. The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations. J. Phys. Chem. B 2007, 111, 7812– 7824.

(51) Park, S.; Bardhan, J. P.; Roux, B.; Makowski, L. Simulated X-Ray Scattering of Protein Solutions Using Explicit-Solvent Models. J. Chem. Phys. 2009, 130, 134114.

(64) Monticelli, L.; Kandasamy, S. K.; Periole, X.; Larson, R. G.; Tieleman, D. P.; Marrink, S.-J. The MARTINI Coarse-Grained Force Field: Extension to Proteins. J. Chem. Theory Comput. 2008, 4, 819–834.

(52) Chen, P.-c.; Hub, J. S. Interpretation of Solution X-Ray Scattering by Explicit-Solvent Molecular Dynamics. Biophys. J. 2015, 108, 2573–2584.


16

Page 17 of 19 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


(65) Ulmer, T. S.; Ramirez, B. E.; Delaglio, F.; Bax, A. Evaluation of Backbone Proton Positions and Dynamics in a Small Protein by Liquid Crystal NMR Spectroscopy. J. Am. Chem. Soc. 2003, 125, 9179–9191.

(75) Bomblies, R.; Luitz, M. P.; Scanu, S.; Madl, T.; Zacharias, M. Transient Helicity in Intrinsically Disordered Axin-1 Studied by NMR Spectroscopy and Molecular Dynamics Simulations. PLoS ONE 2017, 12, e0174337.

(66) Vijay-Kumar, S.; Bugg, C. E.; Cook, W. J. Structure of Ubiquitin Refined at 1.8˚ Aresolution. J. Mol. Biol. 1987, 194, 531–544.

(76) Fushman, D.; Cowburn, D. In Methods in Enzymology; James, T. L., Dötsch, V., Schmitz, U., Eds.; Nuclear Magnetic Resonance of Biological Macromolecules - Part B; Academic Press, 2001; Vol. 339; pp 109–126.

(67) Massi, F.; Grey, M. J.; Palmer, A. G. Microsecond Timescale Backbone Conformational Dynamics in Ubiquitin Studied with NMR R1ρ Relaxation Experiments. Protein Sci. 2005, 14, 735–742.

(77) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High Performance Molecular Simulations through Multi-Level Parallelism from Laptops to Supercomputers. SoftwareX 2015, 1–2, 19–25.

(68) Hall, J. B.; Fushman, D. Variability of the 15N Chemical Shielding Tensors in the B3 Domain of Protein G from 15N Relaxation Measurements at Several Fields. Implications for Backbone Order Parameters. J. Am. Chem. Soc. 2006, 128, 7855– 7870.

(78) Gil-Ley, A.; Bottaro, S.; Bussi, G. Empirical Corrections to the Amber RNA Force Field with Target Metadynamics. J. Chem. Theory Comput. 2016, 12, 2790–2798.

(69) Schmid, N.; Eichenberger, A. P.; Choutko, A.; Riniker, S.; Winger, M.; Mark, A. E.; van Gunsteren, W. F. Definition and Testing of the GROMOS Force-Field Versions 54A7 and 54B7. Eur Biophys J 2011, 40, 843.

(79) Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved Side-Chain Torsion Potentials for the Amber ff99SB Protein Force Field. Proteins 2010, 78, 1950–1958.

(70) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of Multiple Amber Force Fields and Development of Improved Protein Backbone Parameters. Proteins 2006, 65, 712–725.

(80) Piana, S.; Lindorff-Larsen, K.; Shaw, D. E. How Robust Are Protein Folding Simulations with Respect to Force Field Parameterization? Biophys. J. 2011, 100, L47–L49.

(71) Best, R. B.; Mittal, J.; Feig, M.; MacKerell, A. D., Jr. Inclusion of Many-Body Effects in the Additive CHARMM Protein CMAP Potential Results in Enhanced Cooperativity of Alpha-Helix and Beta-Hairpin Formation. Biophys. J. 2012, 103, 1045–1051.

(81) Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New Feathers for an Old Bird. Comput. Phys. Commun. 2014, 185, 604–613. (82) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graph. 1996, 14, 33–38.

(72) Hall, J. B.; Fushman, D. Characterization of the Overall and Local Dynamics of a Protein with Intermediate Rotational Anisotropy: Differentiating between Conformational Exchange and Anisotropic Diffusion in the B3 Domain of Protein G. J. Biomol. NMR 2003, 27, 261–275.

(83) Frishman, D.; Argos, P. Knowledge-Based Protein Secondary Structure Assignment. Proteins 1995, 23, 566–579. (84) Lindorff-Larsen, K.; Best, R. B.; DePristo, M. A.; Dobson, C. M.; Vendruscolo, M. Simultaneous Determination of Protein Structure and Dynamics. Nature 2005, 433, 128–132.

(73) Damberg, P.; Jarvet, J.; Gräslund, A. Limited Variations in 15N CSA Magnitudes and Orientations in Ubiquitin Are Revealed by Joint Analysis of Longitudinal and Transverse NMR Relaxation. J. Am. Chem. Soc. 2005, 127, 1995–2005.

(85) McGibbon, R. T.; Beauchamp, K. A.; Harrigan, M. P.; Klein, C.; Swails, J. M.; Hernández, C. X.; Schwantes, C. R.; Wang, L.P.; Lane, T. J.; Pande, V. S. MDTraj: A Modern Open Library for the Analysis of Molecular Dynamics Trajectories. Biophys. J. 2015, 109, 1528– 1532.

(74) Lange, A.; Casta˜ neda, C.; Hoeller, D.; Lancelin, J.-M.; Fushman, D.; Walker, O. Evidence for Cooperative and Domain-Specific Binding of the Signal Transducing Adaptor Molecule 2 (STAM2) to Lys63-Linked Diubiquitin. J. Biol. Chem. 2012, 287, 18687–18699.


17


(86) Feenstra, K. A.; Hess, B.; Berendsen, H. J. C. Improving Efficiency of Large Time-Scale Molecular Dynamics Simulations of Hydrogen-Rich Systems. J. Comput. Chem. 1999, 20, 786–798. (87) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577–8593. (88) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling Through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101. (89) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695–1697. (90) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182–7190. (91) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684–3690.


18

Page 18 of 19

Page 19 of 19

2 R1

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


1.5 expt MD 1

coarse-grained global motion

+

all-atom local motion

For Table of Contents Only


19

Ab-initio Prediction of NMR Spin-Relaxation Parameters from

Recommend Documents