Computational Tool for Interpreting NMR Spin

J. Phys. Chem. B 2009, 113, 13613–13625

13613

General Theoretical/Computational Tool for Interpreting NMR Spin Relaxation in Proteins Mirco Zerbetto,† Antonino Polimeno,*,† and Eva Meirovitch*,‡ Dipartimento di Scienze Chimiche, UniVersita` di PadoVa, PadoVa, Italy, and The Mina & EVerard Goodman Faculty of Life Sciences, Bar-Ilan UniVersity, Ramat-Gan, Israel ReceiVed: May 19, 2009; ReVised Manuscript ReceiVed: August 24, 2009

We developed in recent years the slowly relaxing local structure (SRLS) approach for analyzing NMR spin relaxation in proteins. SRLS is a two-body coupled rotator model which accounts rigorously for mode-coupling between the global motion of the protein and the local motion of the spin-bearing probe and allows for general properties of the second rank tensors involved. We showed that a general tool of data analysis requires both capabilities. Several important functionalities were missing in our previous implementations of SRLS in data fitting schemes, and in some important cases, the calculations were tedious. Here we present a general implementation which allows for asymmetric local and global diffusion tensors, distinct local ordering and local diffusion frames, and features a rhombic local potential which includes Wigner matrix element terms of ranks 2 and 4. A recently developed hydrodynamics-based approach for calculating global diffusion tensors has been incorporated into the data-fitting scheme. The computational efficiency of the latter has been increased significantly through object-oriented programming within the scope of the C++ programming language, and code parallelization. A convenient graphical user interface is provided. Currently autocorrelated 15N spin relaxation data can be analyzed effectively. Adaptation to any autocorrelated and cross-correlated relaxation analysis is straightforward. New physical insight is gleaned on largely preserved local structure in solution, even in chain segments which experience slow local motion. Prospects associated with improved dynamic models, and new applications made possible by the current implementation of SRLS, are delineated. 1. Introduction The dynamic properties of biological macromolecules are as critical for their function as their three-dimensional structure. One can elucidate protein and DNA/RNA dynamics by using either atomistic or mesoscopic approaches. The most advanced atomistic approach is molecular dynamics (MD), often used within the scope of advanced sampling techniques.1 This method can be applied currently to medium size proteins with the aqueous solvent accounted for explicitly, and simulation lengths extending up to microseconds.2 However, the modeling of water is still deficient to the extent that the global diffusion tensors are not determined satisfactorily.3 Also, the derivation of macroscopic properties (diffusion and ordering tensors and relevant features of local geometry) from the MD time correlation functions is problematic. Finally, force-field imperfections are still a matter of concern.3,4 The other type of approach is mesoscopic in nature. The manybody coupled rotator slowly relaxing local structure (SRLS) approach is such a method.5,6 It was applied as a two-body model in the context of electron spin resonance (ESR) of nitroxide probes dissolved in complex fluids7,8 and later to nitroxide-labeled biological macromolecules.9-11 Recently, we applied it to NMR spin relaxation in proteins,12,13 with the coupled rotators representing the locally reorienting spin-bearing probe (typically an 15N-1H bond14-16 or a deuterated methyl group17,18) and the globally reorienting protein. This application is the subject matter of the present study. * Corresponding authors. E-mail: [email protected] (A.P.); [email protected] (E.M.). † Universita` di Padova. ‡ Bar-Ilan University.

SRLS accounts rigorously for dynamical coupling between the local and global motions through a local ordering potential.5,6 Their time scale separation, and the strength of the coupling/ordering potential, is arbitrary. SRLS allows for general properties of the second rank tensors involved.6,9,13 Thereby general features of local geometries and asymmetric tensors can be treated. We found in previous work that SRLS generates physically insightful pictures of structural dynamics by NMR.12-18 However, several important functionalities were still missing in our previous implementations of SRLS in data fitting schemes,12,13 and in some important cases, the calculations were tedious. Here, we present a general and effective implementation which matches the sensitivity of typical experimental data sets. We illustrate this tool in the context of autocorrelated 15N spin relaxation of amide bonds, although it applies to any spin bearing moiety attached to a protein or a nucleic acid fragment, which experiences auto- or cross-correlated spin relaxation. The traditional method for analyzing NMR spin relaxation in proteins is model-free (MF).19-21 This method, which does not feature the important capabilities mentioned above, is a limiting case of SRLS.5,13,22 On the basis of comparison of results obtained by applying the full theory and its limiting case to the same experimental data, we found that MF has a small range of validity (see below). By starting with the MF limit and systematically decreasing tensor symmetry and increasing geometric complexity, we determined which parameter combination matches data sensitivity. This feature is particularly important in NMR spin relaxation analyses, where the spectral densities underlying the expressions for the experimental relaxation parameters are typically given by linear combinations of Lorentzian functions. If the theoretical spectral density is a good approximation to the experimental spectral density, then one obtains physically meaningful best-

10.1021/jp9046819 CCC: $40.75  2009 American Chemical Society Published on Web 09/23/2009

13614

J. Phys. Chem. B, Vol. 113, No. 41, 2009

fit parameters with good statistics by least-squares fitting the experimental spectral density to the theoretical spectral density. However, one can often fit the experimental data with good statistics using a linear combination of Lorentzian functions which represents a physically implausible situation. This leads to inaccurate best-fit parameters which have absorbed important unaccounted for factors. This process is called force-fitting. There are cases in which the addition of a variable to the fitting process (with the total number of variables being obviously equal to, or less than, the number of experimental data points) leads to an increase in the residuals. In such instances, the added variable is likely to be correlated with another variable. There are also cases in which unduly large uncertainties in the best-fit parameters are observed (similar problems have been encountered in MF analyses).23 Manifestations of this kind are usually termed “over-fitting”. We encountered both types of inappropriate data fitting processes in our work. Such calculations were not considered in deriving physical pictures. We found that a general tool of data analysis must be able to account for mode-coupling and realistic tensor symmetries.13 It must also account for the fact that in the presence of a local potential the Wigner functions are no longer eigenfunctions of the uncoupled local diffusion operator. A method that meets all of these requirements cannot yield analytical solutions, not even in the limit where statistical independence (or modedecoupling) between the local and global motions is valid within a good approximation. Currently the coupled-rotator methods are unique in fulfilling the requirements outlined above. In the field of protein dynamics by NMR the term “mode-coupling” is often referred to as a change in the shape of the protein brought about by local motions (e.g., ref 24). This differs from the meaning of this concept in traditional “cage”-type liquid dynamics theories, to which SRLS belongs. Let us consider first the “mode-decoupling” limit.25,26 A particle reorients in a solvent cage, the structure of which relaxes slowly as compared to the motion of the particle. This problem was solved in early work for simple tensorial properties of the system, in particular weak axial local ordering.25 The solution showed that the particle reorients freely except for a small (negative) correction term in the measurable spectral density.25 In the fully developed SRLS approach5,6 the cage and the particle became two rotators, with their rotational degrees of freedom “coupled” by a local potential. In the case of NMR spin relaxation in proteins, the two rotators are the protein and the probe, with the local potential representing the spatial restrictions imposed by the immediate (internal) protein surroundings at the site of the motion of the probe. May the tensorial properties be simplified in treating the restricted local motions which actually occur in proteins? This question can be answered by consulting the literature on spin relaxation in liquid crystals (LC), where the restricted motion of the “uncoupled” rotator is treated (the SRLS potential has the same form as the mean-field LC potential).27-31 These studies have shown that general tensorial properties are important. We found that this conclusion also applies to local motions in proteins.13,17,18 SRLS solves a two-body Smoluchowski equation.5 Similar to the LC scenario, 25 time correlation functions CKK′(t), with K ) 0, 1, 2 and K′ ) 0, 1, 2, are obtained (K is the order of the rank two local ordering/local diffusion tensor).6 Only 6 (9) distinct functions persist for axial (rhombic) magnetic tensors. Each function comprises an infinite number of terms, which depend on a small number of parameters (with formal MF analogues). In

Zerbetto et al. practice, a finite number of terms suffice for convergence of the SRLS solution. The CKK′(t) functions are linearly combined into the measurable time correlation function, C(t), in accordance to the local geometry (see below). The measurable MF time correlation function represents a two-term C00(t) function. This is actually shown in refs 32 and 33. For large time scale separation and weak local ordering, the form of C00(t) and the limiting form of SRLS25 are the same. For large time scale separation and strong local ordering, similarity with a limiting form of SRLS also exists.13,22 In ref 34, forms for the K ) 1 and K ) 2 components have also been developed. In ref 19, the physical parameters of the C00(t) function have been generalized in a mathematical manner. This clarifies the assessment that the MF spectral density has a small range of validity. Therefore, the theoretical spectral density is typically oversimplified in MF analyses. Similar considerations apply to the EMF spectral density,21 which is another simple limit of SRLS.13,26 The SRLS potential is diffusive in nature. It can be simplified to represent discrete jumps or wobble-in-a-cone. When the L ) 4, K ) 0, 2 terms are included, it can also treat diffusion in multiple wells, leading to jumps in the limit of very high barriers (see Figure 4 of ref 6). Diffusion in these potentials is significantly more general in nature35-37 than discrete jumps (e.g., see ref 38). To our knowledge, SRLS is the only currently available approach which has implemented “mode-coupling” as outlined above in the context of general tensorial properties and a versatile local potential. For 15N spin relaxation, one may vary at present 4-5 parameters, using typical (two-filed) experimental data sets. This level of parameterization of the SRLS model was found to be sufficient in most cases to provide consistent analyses of the available experimental data. Enhancements require incorporation of outside information, e.g., provided by MD simulations, quantum chemical calculations, and/or hydrodynamic modeling. Significant advances have been made lately in developing such integrated approaches.39-48 However, the generation of general and effective integrated tools is still a challenge. Hence it is timely, important, and useful to develop a general, effective, and convenient theoretical/computational tool based on SRLS for analyzing NMR spin relaxation in proteins. This objective has been accomplished herein in the form of the software package called C++OPPS (COupled Protein Probe Smoluchowski). The SRLS program has been enhanced to allow for rhombic local and global diffusion tensors, separate local ordering and local diffusion frames, and a local potential which includes Wigner matrix element terms up to rank four. The handling of the local geometry has been reformulated to make possible treating crosscorrelated relaxation besides autocorrelated relaxation. Importantly, the fitting scheme for SRLS has been integrated with an independent hydrodynamic program for calculating rhombic global diffusion tensors from 3D structures.49 Finally, efficiency has been improved significantly by using the C++ programming language, which among others allows for dynamic allocation of memory and object oriented programming, and by parallelizing the timeconsuming parts of the program. A convenient graphical user interface (GUI) is provided. Currently C++OPPS is designed to analyze 15N spin relaxation from 15N-1H amide bond vectors. It can be extended in a straightforward manner to additional autocorrelated relaxation scenarios and to cross-correlated relaxation. First results obtained with C++OPPS already provide new physical insight. Thus, we find that the local ordering is high not only for “rigid”

Interpreting NMR Spin Relaxation with SRLS

J. Phys. Chem. B, Vol. 113, No. 41, 2009 13615 M2F into OF. ΩD transforms OF into DF, and ΩOC transforms OF into CF; ΩC transforms DF into CF. In C++OPPS, we use ΩL and Ω as the set of stochastic variables;13 in earlier work, we used ΩL and ΩLO.12 All the other Euler angles are time-independent angles. Denoting the stochastic variables as X ) (ΩL,Ω), the Smoluchowski operator describing the time evolution of the probability density of the system is given by:

Γˆ (X) ) OJˆtr(Ω)OD2Peq(X)OJˆPeq-1(X) + [VJˆ(Ω) Vˆ J(ΩL)]trVD1Peq(X)[VJˆ(Ω) - VJˆ(ΩL)]Peq-1(X) (1) Figure 1. Representation of reference frames involved in the definition of the SRLS model.

N-H bonds, but also for N-H bonds which experience slow local motion. According to the local geometry determined, the slow motions are chain motions. Aspects of the currently used overdamped (Smoluchowski) picture which require improvement are identified, and solutions are suggested. Prospects associated with new applications, made possible by the SRLS implementation offered herein, are delineated. We distribute C++OPPS under the GNU Public License (GPL) v2.0. The software is available at the website http:// www.chimica/unipd.it/licc/software.html. 2. Theoretical Background The basics of the SRLS approach to NMR spin relaxation in proteins, as applied in this study, are summarized below. Further details are given in ref 13. SRLS solves a two-body Smoluchowski equation which describes the time evolution of the probability density of two coupled rotators.5,6 Dynamical coupling is brought into effect by a mean-field potential. In principle, the rotational diffusion rates of the coupled rotators are arbitrary. In practice, the unrestricted motion of the protein is slower than, or at most comparable to, the restricted motion of the probe. The coupling potential, exerted by the immediate (internal) protein surroundings at the site of the motion of the probe, represents the local ordering potential. The SRLS order parameters are defined in terms of this potential. Thus, the mesoscopic SRLS approach features dissipative parameters in the form of diffusion tensors and structural parameters in the form of potential coefficients, ordering tensors, and Euler angles which represent relative tensor orientations. The reference frames underlying SRLS are illustrated in Figure 1. LF is the inertial laboratory frame, M1F and M2F are, respectively, the protein and probe fixed frames where the diffusion tensors D1 (protein) and D2 (probe) are diagonal. VF is the protein-fixed local director frame, typically associated with the equilibrium orientation of the probe. OF is the probe-fixed local ordering frame. DF is the probe-fixed frame where the 15 N-1H dipolar tensor is diagonal. CF is the probe-fixed frame where the 15N chemical shift (CSA) tensor is diagonal. 15N-1H dipolar/15N CSA cross-correlated relaxation involves the DF and CF frames.12,50 Frame transformations are carried out as follows: ΩL, the stochastic angles which represent the diffusion of the protein with respect to the laboratory frame, transforms LF into VF; ΩLO transforms LF into OF and represents the diffusion of the probe with respect to the laboratory frame. Ω represents the diffusion of the probe with respect to the protein and transforms VF into OF. ΩV transforms M1F into VF, and ΩO transforms

where OD2 denotes the diffusion tensor of the probe in the OF frame, and VD1, the diffusion tensor of the protein in the VF frame. aJˆ(Ω) is the angular momentum operator expressed in frame a describing an infinitesimal rotation of the angles Ω. The equilibrium probability distribution, Peq(X), is given by

Peq(X) ) exp[-V(X)/kBT]/〈exp[-V(X)/kBT]〉

(2)

where kB is the Boltzmann constant and T is the absolute temperature. For a protein in isotropic solution, the equilibrium probability distribution is independent of ΩL. The local potential is expanded in the full basis set of the Wigner rotation elements. Preserving only the L ) 2 and L ) 4 terms, one obtains:

-V(Ω)/kBT ) c02D0,02(Ω) + c22[D0,-22(Ω) + D0,22(Ω)] + c04D0,04(Ω) + c24[D0,-24(Ω) + D0,24(Ω)] + c44[D0,-44(Ω) + D0,44(Ω)] (3) The SRLS solution consists of time correlation functions, CKK′(t), where K and K′ denote rank two tensor order. The spectral density functions, jKK′(ω), are Fourier-Laplace transformations of CKK′(t). The stochastic variables that enter the calculation of jKK′(ω) are the Euler angles ΩLO ) Ω + ΩL. One obtains the expressions: J J' jK,K'(ω) ) 〈DM,K (ΩLO)|(iω - Γˆ )-1 |DM',K' (ΩLO)Peq(ΩLO)〉 (4)

Since the medium is isotropic M ) M′ ) 0. Since all the tensors involved are symmetric and have rank two, and rotational reorientation is associated with second rank time correlation functions, only the terms with J ) J′ ) 2 are considered. The jKK′(ω) functions are calculated in the direct product basis set given by

|Λ〉 ) |λ1〉 X |λ2〉 ) |L1M1K1〉 X |L2M2K2〉

(5)

where:

|L1M1K1〉 ) |L2M2K2〉 ) and [L] ) 2L + 1.

[L1] 2

8π [L2]

L1 DM (ΩL) 1,K1

L2 DM (Ω) 2 2,K2

8π

(6)

13616


Zerbetto et al.

This basis set is infinite. However, the SRLS solution converges for a finite subset determined by a truncation parameter, L2,max. The dimension of the basis set is given by L2,max

N)5

∑ i)0

5 (2i + 1)2) (L2,max + 1)(2L2,max + 3 1)(2L2,max + 3) (7)

It is convenient to define the function 2CK,K' ) D0,K2 + D0,K'2 and calculate the symmetrized spectral densities as S jK,K' ) 〈CK,K'(ΩLO)|(iω - Γˆ )-1 |CK,K'(ΩLO)Peq(ΩLO)〉

(8) The jK,K′(ω) functions are then obtained as linear combinations of the symmetrized spectral densities: S S S jK,K'(ω) ) [2(1 + δK,K')jK,K' (ω) - jK,K (ω) - jK',K' (ω)]/2 (9)

Using the closure relation for the basis set |Λ〉, the integral in eq 8 can be rewritten in matrix form as S jK,K' ) vtr(iω1 - Γ)v

(10)

Figure 2. Schematic illustrating of the structure of the C++OPPS package. The red boxes represent the input data, the rounded blue boxes represent the main calculations, and the green hexagons represent the output of the main calculations. The magenta circle represents the nonlinear least-squares fitting procedure.

1

S jK,K',n (ω) )

where:

β12

R1 - iω (Γ)ij ) 〈Λi |Γˆ |Λj〉 (v)i ) 〈Λi |CK,K'(ΩLO)Peq(ΩLO)〉

R2 - iω -

(11)

We employ Lanczos tridiagonalization,51,52 which is an iterative algorithm creating at every step n a tridiagonal symmetric matrix, Tn. The latter is an approximation to the matrix representation of the Smoluchowski operator, Γˆ . The transformation between the approximate and exact matrices is carried out by the orthonormal matrix Qn:

β22 R3 - iω - ...

(15) where r is the n-dimensional vector containing the diagonal of Tn, and β is the (n - 1)-dimensional vector containing the subdiagonal of the tridiagonal matrix. The measurable spectral density for autocorrelated relaxation is given by 2

J (ω) ) µµ

Tn ) Qn ΓQn tr

(12)

∑

[D*K,02(Ωµ)DK',02(Ωµ)]jK,K'(ω)

K,K')-2

(16) S (ω) Thus, at any given iteration step n, the approximant jK,K′,n S to jK,K′(ω) is given by

S jK,K',n (ω) ) (Qnv)tr(iω1n + Tn)-1(Qnv)

(13)

Based on the orthonormality of Qn, the last equation reduces

where Ωµ represents the Euler angles transforming the OF frame into the respective magnetic frame. This expression for Jµµ(ω) corresponds to the treatment of NMR spin relaxation in the Redfield limit by second-order perturbation theory.31 The angles Ωµ represent what we call the “local geometry”. We actually implemented in C++OPPS the following equation:

to 2

S jK,K',n (ω)

) (iω1n + Tn)1,1

-1

(14)

Jµν(ω) )

∑

[D*K,02(Ωµ)DK',02(Ων)]jK,K'(ω)

K,K')-2

(16a) i.e., the nth approximant to the symmetrized spectral density is the (1,1) element of the inverse of the matrix (iω1n + Tn). S (ω) can be expressed as a continued fraction: jK,K′,n

which yields the autocorrelated spectral density, JDD(ω), for µ ) ν ) D, and the 15N-1H dipolar/15N CSA cross-correlated


J. Phys. Chem. B, Vol. 113, No. 41, 2009 13617

spectral density, JDC(ω), for µ ) D and ν ) C. The involved (from the theoretical and practical point-of-view) calculation concerns the jK,K′(ω) functions; the difference between JDD(ω) and JDC(ω) only concerns the (time-independent) coefficients in the respective linear combinations. To calculate JCC(ω) one has to apply an additional frame transformation from the D frame to the C frame, involving the (timeindependent) angles ΩC (Figure 1). Note that in order to calculate JDC(ω) and JCC(ω) in general all three components, j00(ω), j11(ω), and j22(ω), are required. The model-free approach does not feature j11(ω) and j22(ω). Therefore, JCC(ω) and JDC(ω) can be calculated properly in MF only in the limit in which j11(ω) and j22(ω) may be ignored. This limit, which represents the extreme motional narrowing limit for the local motion, is specified in ref 53. Thus, the SRLS spectral densities jK,K′(ω) are the building blocks of a given dynamic model and the measurable spectral densities Jν(ω) are the building blocks of a specific geometric implementation of this model. For autocorrelated 15N spin relaxation, the measurable quantities are JDD(0), JDD(ωN), JDD(ωH), JDD(ωH + ωN), and JDD(ωH - ωN) for the dipolar 15N-1H interaction and JCC(0) and JCC(ωN) for 15N CSA interaction. Together with the magnetic interactions, these functions determine the experimentally measured relaxation rates given by54,55

1 ) d2[JDD(ωH - ωN) + 3JDD(-ωN) + 6JDD(ωH + ωN)] + c2JCC(-ωN) T1 1 d2 c2 ) [4JDD(0) + JDD(ωH - ωN) + 3JDD(-ωN) + 3JDD(-ωN) + 6JDD(ωH + ωN)] + [3JCC(-ωN) + 4JCC(0)] T2 2 6 γ H NOE ) 1 + d2 T1[6JDD(ωH + ωH) - JDD(ωH - ωH)] γN

(17)

where d2 ) 0.1γH2γN2p2〈rNH-3〉, c2 ) 2/15ωN2∆σ2 (e.g., ref 56), and the various constants have their usual connotation. 3. Software 3.A. Implementation. C++OPPS is a modular package of programs. The core programs, written in C++, perform the heavy computations, and the graphical user interface (GUI) is the front-end. In general, there are three categories of input: a 3D structure of the protein studied, a file comprising the physical parameters, and a file comprising the experimental data. On the basis of this information, C++OPPS can be used for the following purposes, which are represented schematically in Figure 2. Calculation of the Time Correlation Functions, CK,K′(t), and the Generic Spectral Densities, jK,K′(ω). This requires the SRLSrelated parameters residing in the physical input file (Figure 2). C++OPPS distinguishes between the local ordering and local diffusion frames; in the FORTRAN program,13 these frames were taken to be the same. This distinction is important when the main ordering axis is tilted relative to the principal axis of the dipolar tensor and the symmetry around the main ordering axis is nearly axial. In such cases, permutations of ordering axes labels to set the Z axis along the main ordering axis renders the ordering tensor nearly axially symmetric.57 Thereby, physical insight is gained, and the calculation is sped up significantly. On the other hand, since the principal axis of the dipolar frame is usually collinear with the cylindrical shape of the probe, there is no a priori reason to redefine the local diffusion tensor frame.

Calculation of the Relaxation Parameters. All the parameters comprised in the physical input file are required to calculate the generic jKK′(ω) functions. To assemble Jµµ(ω), or in general Jµν(ω), from the jKK′(ω) functions, the local-geometry-related information, which consists of ΩD and ΩC, is also required. To calculate the relaxation parameters from the Jµµ(ω) functions, the magnetic interactions are also required. The information associated with the global diffusion tensor can be (a) included in the physical input file if known from independent sources, (b) determined with DITE49 and transferred to the physical input file, or (c) fit together with the site-specific variables (Figure 2). Data Fitting. This calculation requires the information specified in the previous section and experimental data. The current fitting scheme features SRLS integrated with the MINPACK optimization package. Since SRLS is general in nature, various parameter combinations can be taken as a variable set. Obviously we require low (according to common criteria) χ2 for result acceptance. However, we also require physical viability of the best-fit parameters obtained, in the context of the dynamic probe used and its immediate protein surroundings. The sets of variable parameters are determined by systematically decreasing tensor symmetry and increasing geometric complexity. This process is pursued until further enhancements do not alter the results significantly and/or until symptoms of overfitting are observed. We found that for N-H

TABLE 1: Input Values for the Calculation of 15N T1, T2 and 15N-{1H} NOE with C++OPPSa D1/Hz

D2/Hz

1.1 × 10 1.1 × 107 1.1 × 107 7

5.19 × 10 5.19 × 107 4.56 × 108 7

F2

βD/deg

c02

c22

T1/ms

T2/ms

NOE

tC++/s

tF77/s

9.6 9.6 100.0

101.4 101.3 101.3

5.7 5.7 5.7

10.5 20.0 20.0

933.91 1281.7 1251.6

44.65 36.16 37.86

0.6581 0.7259 0.7775

112 125 108

1469 1956 1967

a D1 is the diffusion constant for isotropic global diffusion. D2 ) tr{D2}/3 is the “isotropic” rate constant of the axial local diffusion tensor, and F2 ) D2,||/D2,⊥ measures its axiality. c02 and c22 denote the coefficients of the L ) 2 axial and rhombic terms of the local potential, respectively. βD is the angle between the Z-axis of the axial dipolar frame and the actual main ordering axis of the rhombic local ordering frame. The angle γD was set equal to 90°,13 and the Euler angles ΩC to (0°, 17°, 0°).12 A magnetic field of 14.1 T and L2,max ) 30 were used. The columns denoted tC++/s and tF77 show the times elapsed. The calculations were carried out on an HP computer equipped with an Intel 2.7 GHz Dual Core CPU and 4 GB of RAM.

13618


Zerbetto et al.

TABLE 2: Relative Difference {[X(axial) - X(rhombic)]/ X(axial)} × 100 between X ) 15N T1, T2 and 15N-{1H} NOE Calculated for τm ) 1/6D1 ) 15 ns, a Time Scale Separation D1/D2 ) 0.01, ΩD ) (0°, 0°, 0°), ΩC ) (0°, 17°, 0°), and an Axial (c02 ) 8 and c22 ) 0) or a Rhombic (c02 ) 8 and c22 ) 4) Potentiala 15

N T1 N T2 15 N-{1H} NOE 15

11.7 T

14.1 T

-2.7 -7.4 30.0

-1.3 -7.5 40.7

TABLE 3: Relative Difference {[X(βV ) 0°) - X(βV ) 90°)]/X(βV ) 0°)} × 100 between X ) 15N T1, T2 and 15 N-{1H} NOE Calculated with τm ) 1/6D1 ) 15 ns, D1,⊥/D2 ) 0.01, an axial potential with c02 ) 8, ΩD ) (0°, 0°, 0°), ΩC ) (0°, 17°, 0°), and an Axial Global Diffusion Tensor with D1,||/D1,⊥ ) 1.2a

18.8 T 1.2 -7.6 44.2

15

N T1 N T2 15 N-{1H} NOE 15

11.7 T

14.1 T

18.8 T

7.4 -9.0 -2.8

7.1 -9.0 -3.6

6.1 -9.2 -4.1

a Calculations are shown for magnetic fields of 11.7, 14.1, and 18.8 T. The truncation parameter was set equal to 20.

a Calculations are shown for magnetic fields of 11.7, 14.1, and 18.8 T. The truncation parameter was set equal to 20.

bond dynamics and combined two-field 15N T1 or T2 and 15 N-{1H} NOE data (i.e., six data points), one may vary without encountering symptoms of overfitting 4-5 parameters. The fitting process is represented in Figure 2 by the magenta circle. One may use single-field data comprising 15N T1 or T2 and 15 N-{1H} NOE; clearly at most, three parameters can be varied in this case. Calculation of the Global Diffusion Tensor by DITE.49 This computation requires a 3D structure in Z-matrix (ZMT) or protein data bank (PDB) format from which the program Babel extracts the required input. The DITE approach can also determine global diffusion tensors in the presence of torsions about a small number of internal degrees of freedom. This important and unique functionality will be implemented in future work and will be used to devise strategies for calculating global diffusion tensors of multidomain proteins. 3.B. Efficiency. The matrix representation of the Smoluchowski operator (eq 11) can reach dimensions on the order of 104-105. This requires a substantial amount of RAM memory for matrix element storage. The required linear algebra operations are time-consuming. Hence the efficiency of the numerical calculation is important. The factors underlying the increased efficiency of C++OPPS relative to our older FORTRAN program13 have been discussed above. Table 1 illustrates the performance of these programs in the context of an example featuring isotropic global diffusion, axial local diffusion, and rhombic local potential. The data in the first row are typical of N-H bonds located within the mobile domain AMPbd of E. coli adenylate kinase (AKeco).13 The input parameters represent the best-fit parameters obtained with the FORTRAN-programbased fitting scheme (see data for residue 46 in Table 18 of ref 13). We could not analyze with data fitting the CORE domain using the FORTRAN program because the combination of high rhombic ordering and large time scale separation rendered the calculation impractical for repeated calculations inherent to fitting procedures. The input parameters in the third row of Table 1 represent this scenario within the scope of a single calculation.

on the local ordering (note that in MF the local ordering enters the calculation through the squared generalized order parameter, S2(ref 19)). The rhombicity of the local ordering, which affects the NOE to such a large extent, is quite limited. This can be appreciated by calculating the Cartesian ordering tensor components from c02 ) 8 and c22 ) 4 (see, for example, ref 13 for the relevant procedure). These components are given by Sxx ) -0.382, Syy ) -0.454, and Szz ) 0.836, yielding (Sxx - Syy)/Szz ) 0.09 on a scale extending from -1 to +1. Note that the extent of the local ordering is high, with Szz ) 0.836. Table 3 illustrates limited sensitivity of the analysis to small global diffusion axiality, given by F1 ) D1|/D1⊥ ) 1.2, as one would expect. This is in stark contrast with the large effect F1 ) 1.18 has on the MF analysis of dihydrofolate reductase (DHFR).58 Thus, 50% of the residues of this protein require substantial conformational exchange contributions, Rex, when isotropic D1 is used instead of axial D1 with F1 ) 1.18.58 If, however, D1 is allowed to be axial, then the Rex contributions disappear and the unaccounted for rhombicity of the S tensor (see next paragraph), which has a large effect on the analysis (Table 2), is absorbed by an apparent axiality of the D1 tensor. Evidence that Rex absorbs the unaccounted for asymmetry of the local ordering appears in ref 50 where Ribonuclease H (RNase H) and AKeco have been studied in this context. It is also shown in that study that using axial potentials instead of (appropriate) rhombic potentials, and an axial global diffusion tensor instead of an (appropriate) isotropic global diffusion tensor, implies inaccurate order parameters obtained with data fitting. The findings of ref 50 are based on extensive predictive calculations, and back-calculations of experimental data, carried out in the context of a concerted analysis of the autocorrelated relaxation parameters, 15N T1, T2 and 15N-1H NOE and the transverse 15N-1H dipolar/15N CSA cross-correlated relaxation rate, ηxy. The fact that the experimental value of 1/T2 (ηxy) depends (does not depend) on Rex is a key element in the analysis. The information outlined above for DHFR, comprised in ref 50 and illustrated in Tables 2 and 3, highlights the vulnerability of 15N spin relaxation analysis to force-fitting when a dominant featuresthe rhombicity of the local potentialsis not accounted for. The quantities which are particularly prone to absorb this feature include the Rex contribution and the axiality of the global diffusion tensor (which is the only tensor allowed to have symmetry lower that the highest possible symmetry). Therefore it is important to develop strategies for ascertaining that Rex * 0, and/or the axiality of the tensor D1, are real properties. The fact that 15N spin relaxation is sensitive to the rhombicity of the local ordering was shown by the 3D Gaussian axial fluctuations (GAF) model59 applied to rigid N-H bonds of

4. Illustrative Calculations with Significant Consequences 4.A. Potential Asymmetry/Rhombicity versus Global Diffusion Axiality. Table 2 illustrates the high sensitivity of the analysis to the rhombicity of the local potential. It can be seen that the 15N-{1H} NOE is affected substantially by allowing for rhombic potentials/rhombic ordering instead of imposing axial potentials/axial ordering. The larger effect on the 15 N-{1H} NOE as compared to 15N T1 and T2 is most likely due to the fact that the NOE represents a ratio of two relaxation rates,54 each depending intricately (through the jKK′(ω) functions and their coefficients in the expressions for JDD(ω) and JCC(ω))



TABLE 4: SRLS:Five Dominant Eigenmodes and Corresponding Eigenvalues, 1/τn, (in Units of D2) Contributing to C00(t) for c02 ) 8, 10, and 20 and Time Scale Separations D1/D2 ) 0.02, 0.01, and 0.001.a MF limit:19 local motional eigenvalues given by 3c02,29 global motional eigenvalues given by 6D1/D2, and global and local motional eigenmodes given by (S20)2 and [1 - (S20)2], respectively D1/D2 ) 0.02 SRLS

D1/D2 ) 0.01 MF limit

SRLS

D1/D2 ) 0.001 MF limit

c0 ) 8 0.060 (0.754) 24.00 (0.246)

SRLS

MF limit

2

0.1178 (0.762) 21.87 (0.198) 36.34 (0.017) 40.65 (0.009) 37.49 (0.004)

0.120 (0.754) 24.00 (0.246)

0.0594 (0.758) 21.62 (0.203) 35.97 (0.016) 40.27 (0.009) 37.05 (0.005)

0.1177 (0.808) 28.20 (0.170) 47.89 (0.010) 52.56 (0.006) 67.55 (0.002)

0.120 (0.803) 30.00 (0.197)

0.0594 (0.806) 27.89 (0.172) 47.38 (0.010) 52.06 (0.006) 67.95 (0.002)

0.1177 (0.903) 59.08 (0.091) 111.2 (0.002) 115.4 (0.002) 0.1177 (0.001) a

0.120 (0.901) 60.00 (0.099)

0.0594 (0.903) 58.46 (0.093) 110.1 (0.002) 114.3 (0.002) 158.1 (0.0003)

c02 ) 10 0.060 (0.803) 30.00 (0.197)

c02 ) 20 0.060 (0.901) 60.00 (0.099)

0.006 (75.3) 21.39 (0.206) 35.60 (0.015) 39.92 (0.010) 36.67 (0.007)

0.006 (0.754) 24.00 (0.246)

0.006 (0.803) 27.60 (0.175) 46.91 (0.010) 51.61 (0.007) 67.28 (0.002)

0.006 (0.803) 30.00 (0.197)

0.006 (0.900) 57.90 (0.094) 109.1 (0.003) 113.3 (0.002) 0.005994 (0.000488)

0.006 (0.901) 60.00 (0.099)

The truncation parameter used was L2,max ) 20.

Ubiquitin,60 by SRLS applied to rigid N-H bonds of RNase50 and by extensive application of SRLS to AKeco.13-16,50 A residual dipolar coupling (RDC) study reported on anisotropic probability distribution functions of N-H orientations61 and MD simulations revealed asymmetric N-H fluctuations.62 Hence, one expects the experimental local potential, i.e., the experimental ordering tensor, to be asymmetric. Local motional effects can also impact the accuracy of the global diffusion tensor. The determination of the latter is designed to utilize experimental 15N T1/T2 data from rigid N-H bonds unaffected by local motions. Typically, the complete 15N T1/T2 data set is filtered according to MF criteria.56,63 We found that N-H bonds considered rigid per MF might contain local motional effects absorbed by the value of S2 (ref 13; see also, ref 9). Also, D1 analysis is very sensitive to the nature of the filtered 15N T1/T2 data set if axiality is forced upon an actually isotropic tensor.64 The strategy we have been using to minimize local motional effects in determining D1 from 15N T1/T2 ratios is outlined in ref 13. On the basis of the evidence outlined above, it can be concluded that in quite a few cases the axiality of the global diffusion tensor has been over-rated in MF analyses. Clearly, the global diffusion tensor might be genuinely axial, affecting the analysis significantly. Our comment relates to cases in which the axiality of the tensor D1 is apparent, stemming typically from unaccounted for rhombicity of the local ordering. 4.B. Large Time Scale Separation and High Ordering Are Prerequisites for MF Validity in the “Relatively-RigidProbe” Regime. The time correlation function ˜

J J CKK'(t) ) 〈DM,K Peq1/2 |e-Γt |DM,K' Peq1/2〉

(18)

can be recast into the form: ∞

CKK'(t) )

∑ c*K,ncK',ne-t/τ

n)0

∞

n

)

(n) -t/τ e ∑ AKK'

n

n)0

(19)

where 1/τn are the eigenvalues of the Smoluchowski operator and the coefficients cK,n are given by cK,n )(uKV*)n. V denotes the matrix of eigenvectors of the Smoluchowski operator and uK denotes the vector of projections of the physical observable J on the basis functions, i.e., (uK)j ) 〈Λj|DM,K Peq1/2〉. We focus on the parameter range typical of N-H bonds located in canonical elements of secondary structure. Although the actual local potential is high and rhombic, we use for illustration potentials which are high and axial. We show in Table 4 SRLS and MF results for potentials given by c02 equal to 8, 10, and 20 and time scale separations, D1/D2, of 0.02, 0.01, and 0.001. The SRLS time correlation functions are formally analogous to the MF time correlation function of ref 19. Table 4 comprises the five SRLS eigenmodes with the largest fractional contributions (numbers in parentheses) aside their corresponding eigenvalues, 1/τn. The latter are given in units of D2. The correlation times, τn, can be calculated in units of seconds by dividing the inverse of the eigenvalue by D2 ) 1.1107 × 107/(D1/D2) Hz. The factor 1.1107 × 107 Hz is equal to 1/(6 × 15.006 ns). For example, the main local motional eigenvalue for D1/D2 ) 0.02 and c02 ) 8 is equal to 82.33 × 10-12 s. The MF data have been calculated as follows. Assuming diffusive local motion, the local motional eigenvalue is given by the “renormalized” eigenvalue, 1/τren;13,40 in units of D2, one has 1/τren ) 3c02. The global motional eigenvalue is given by 6D1/D2. The global motional eigenmode is given by the squared order parameter, (S02)2, corresponding to the c02 value used. Consequently, the local motional eigenmode is given by [1 (S02)2]. τren in seconds is calculated as outlined above for τn. The time correlation functions, C00(t), are shown in Figures 3, 4, and 5 for D1/D2 ) 0.02, 0.01, and 0.001, with c02 set equal to 8, 10, and 20, respectively. In Figure 6, we show C00(t) for D1/D2 ) 0.001 and c02 ) 8, 10, and 20. The X-axes extend up to about three times the longest local motional correlation time. SRLS and MF generate similar results for the global motion. Thus, the plateau values of the C00(t) functions in Figures 3-6 are close to the global motional eigenmodes of (S02)2 ) 0.754

13620


Figure 3. Logarithmic plot of the time correlation function C00(t) calculated for c02 ) 8 and time scale separations D1/D2 of 0.001 (SRLS solid line, MF circles), 0.01 (SRLS dashed line, MF triangles), and 0.02 (SRLS dotted line, MF squares).



for c02 ) 8, (S02)2 ) 0.803 for c02 ) 10, and (S02)2 ) 0.901 for c02 ) 20. The final decay to zero (not shown) agrees with the global motional eigenvalue of ∼15 ns, which appears in Table 4 in units of D2. The data associated with the local motion are not the same. The main SRLS local motional eigenvalue and the MF local

Zerbetto et al.

Figure 6. Logarithmic plot of time correlation function C00(t) calculated for D1/D2 ) 0.001 and c02 ) 8 (SRLS solid line, MF circles), 10 (SRLS dashed line, MF triangles), and 20 (SRLS dotted line, MF squares).

motional eigenvalue differ from 3.6% for c02 ) 20 and D1/D2 ) 0.001, up to 9.7% for c02 ) 8 and D1/D2 ) 0.02. The main local motional eigenmode and the MF local motional eigenmode differ from 5.6% for c02 ) 20 and D1/D2 ) 0.001, up to by 24% for c02 ) 8 and D1/D2 ) 0.02. The extra SRLS local motional eigenmodes make individual contributions which are below 5% of the main SRLS local motional eigenmode. However, the corresponding eigenvalues, 1/τn, are on the order of the main local motional eigenvalue (Table 4). On the whole, the combined effect of the extra local motional eigenmodes on the shape of the SRLS time correlation functions at short times is significant, rendering the short-time decays quite different from their MF counterparts (Figures 3-6; note that the time correlation functions are plotted on a natural logarithm scale). Clearly both large time scale separation and strong potentials should be in place to obtain accurate MF results in the relativelyrigid-probe regime. The physical considerations underlying these requirements are as follows. Large time scale separation is the fundamental assumption that underlies statistical independence between the global and local motions, on which the MF rationale is based.13,19 For very high potentials, C00(t) is approximated well by the global motional eigenmode and a single local motional eigenmode with eigenvalue given by the “renormalized” eigenvalue of ref 29. Simple eigenfunctions, given in refs 29 and 25 agree with the MF formula.19 As the coupling potential is reduced, the correlation function for the relative motion (i.e., for the corresponding Wigner matrix elements) becomes more complex, involving several eigenmodes of this motion.13 This is illustrated in Table 4 and Figures 3-6. When D1/D2 is larger than 0.02, the coefficient c02 smaller than 8, and the asymmetry of the local potential, as well as realistic local geometry, are accounted for, then the discrepancy between SRLS and MF increases further substantially.13 Fitting actual data with oversimplified time correlation functions leads to force-fitting (good statistics but inaccurate best-fit parameters, which have absorbed unaccounted for factors).13 4.C. Contributions Associated with K ) 1 and 2. When the angle βD is not zero, then all three generic spectral densities, jKK(ω), K ) 0, 1, 2, contribute to JDD(ω). Yet, the MF spectral density of ref 19 is given by JDD(ω) ) JCC(ω) ) j00(ω). We use for illustration the input values D1 ) 1.1 × 107 Hz, D1/D2 ) 0.8, c02 ) 3, ΩD ) (0°, 20°, 0°), and ΩC ) (0°, 17°, 0°). The dominant contributions to the time correlation functions, CKK(t), are shown in Table 5, and the generic spectral densities, jKK(ω), are depicted in Figure 7.



TABLE 5: Dominant Eigenmodes (in Parentheses) and Corresponding Eigenvalues (in Units of D2) Contributing to the Time Correlation Functions, CKK(t), K ) 0, 1, 2, for c02 ) 3 and D1/D2 ) 0.8 1 2 3 4 5

C00(t)

C11(t)

C22(t)

0.243 (0.737) 6.08 (0.223) 2.72 (0.030) 3.04 (0.007) 4.27 (0.004)

0.203 (0.793) 5.67 (0.175) 2.64 (0.017) 2.92 (0.007) 3.45 (0.005)

0.140 (0.944) 4.77 (0.040) 3.45 (0.008) 2.49 (0.005) 2.62 (0.003)

calculation is the same as the number of variables in the analogous MF calculation. Hence, the difference between SRLS and MF consists of the manner in which C(t) (hence J(ω)) is calculated. Corresponding SRLS and MF parameters will have the same meaning if the MF spectral densities, given below for convenience, are valid. The original MF spectral density is given by19

J(ω) ) S2τm /(1 + τm2ω2) + (1 - S2)τ|e /(1 + τ|e2ω2), (20) where 1/τ|e ) 1/τm + 1/τe, with τm denoting the correlation time for global motion, τe, the effective correlation time for local motion, and S2, the squared generalized order parameter. This formula requires that τm . τe. Hence, within a good approximation, τ|e may be replaced by τe. The parameter S2 is formally analogous to (S02)2 in SRLS. For high local ordering τe agrees13,22 with the renormalized correlation time for local motion given by τren ) 2 τ/c02, with τ ) 1/(6D2).29 The parameter τren is the local motional correlation time29 obtained by solving a simple diffusion equation valid in the large time scale separation limit for strong axial potentials.28 The EMF spectral density is given by21

Figure 7. Plot of j00(ω) (full line), j11(ω) (dashed line), and j22(ω) (dotted line) spectral densities corresponding to the data shown in Table 5.

For βD ) 20°, the coefficients of j00(ω), j11(ω), and j22(ω) in the expression for JDD(ω) are given by the reduced Wigner matrix elements (d002)2 ) 0.680, 2(d012)2 ) 0.310, and 2(d022)2 ) 0.010, respectively. The contributions of the KK ) (1,1) and (2,2) components (Figure 7) to JDD(ω) are clearly important. These results should be viewed as merely illustrative because they have been obtained assuming axial local potentials,14-16 while actual experimental data feature rhombic local potentials.13 When rhombic potentials were allowed, the angle βD assumed the physically meaningful value of approximately 100°, which is close to the canonical angle of 101.3° between the N-H bond and the Ci-1R - CiR axis.60 The latter has been determined to be the main local ordering/local diffusion axis. In this case, one has (d002)2 ) 0.196, 2(d012)2 ) 0.111, and 2(d022)2 ) 0.693, with j22(ω) dominating JDD(ω).

J(ω) ) Sf2[Ss2τm /(1 + ω2τm2) + (1 - Ss2)τs| /(1 + ω2τs|2)] + (1 - Sf2)τ|f /(1 + ω2τ|f2) (21)

14,15

5. Typical Data Fitting Scenarios Let us denote N-H bonds located in well-structured regions of the protein, notably elements of secondary structure, as “rigid”, and those located in mobile domains, loops, end-chain segments not interacting with the protein surface, etc., as “flexible”. In many cases the MF analysis of rigid N-H bonds (and in some cases of flexible N-H bonds) requires the inclusion of Rex contributions (added to the expression for 1/T2 given in eq 17). SRLS/C++OPPS fitting of data from these three types of N-H bonds is illustrated below. We selected the protein E. coli adenylate kinase (AKeco) as an example because it features all three N-H bond categories.14,15 The 15N relaxation parameters were acquired at 14.1 and 18.8 T and 303 K (BMRB, accession no. 5746). Further details are given in refs 13, 14, and 64. The magnetic interaction parameters used in the calculations presented below are rNH(eff) ) 1.015 Å,65 ∆σ ) -169 ppm,66 and βD ) 17°.67 An Lmax value of 24 was always found to suffice. We also show the results of the corresponding MF analyses. In each case, the number of variables in a given SRLS

where 1/τf| ) 1/τf + 1/τm and 1/τs| ) 1/τs + 1/τm, with τf (τs) denoting the effective correlation time for fast (slow) local motion. Similar to eq 20, this equation also requires large time scale separation between the global motion and all of the local motions. In practice it is used in cases where τs and τm are comparable. Sf2 and Ss2 denote the squared generalized order parameters associated with the local motions. The parameters τf and τs are formally analogous to τ| ) 1/(6D2,|) and τ⊥ ) 1/(6D2,⊥), respectively. The parameters Sf2 and Ss2 can be expressed in terms of (S02)2 and (S22)2.13,22 Often the reduced EMF formula, given by eq 21 with the third term omitted based on the assumption that τf f 0, is used. Example 1: Rigid Residue with Axial Local Ordering. Data from the 15N-1H bond of residue 197 of the CORE domain of AKeco were analyzed with the MF program DYNAMICS.68 MF model 2, where S2 and τe are varied, led to the best-fit values of S2 ) 0.84 (corresponding to c02 ) 12.4, calculated assuming that S2 f (S02)2 ) 〈D002〉) and τe ) 12.7 ps (corresponding to τ ) 78 ps; see the expression for τren given above), with χ2 ) 2. In analogy, we allowed c02 and τ to vary in the SRLS calculation. This led to c02 ) 11.6 ((S02)2 ) 0.83) and τ ) 69 ps, with χ2 ) 0.6. These results are shown in Table 6, rows 1 and 2. The SRLS calculation in this case took ∼19 s to be completed. The differences between the SRLS and MF results are 1.2% for the squared order parameters, 13% for the local motional correlation times, and 6.9% for the potential coefficient c02. Although χ2 ) 0.6 in the SRLS calculation and χ2 ) 2 in the MF calculation, both calculations are considered appropriate since both χ2 values lie below 5.99, which is the percentile value for χ2 distribution for 4 degrees of freedom (six data points and two variables) for a commonly used 5% threshold (Table 39 of ref 69). The differences versus MF stem from (1) accounting in SRLS for the frame transformation between the 15N-1H dipolar and

13622


Zerbetto et al.

TABLE 6: Illustrative SRLS Data-Fitting Calculationsa residue 197 197 209 209 209 46 46 46

method MF SRLS MF SRLS SRLS MF SRLS SRLS

ax ax rh ax rh

c02 ((S02)2)b 12.4 (0.84) 11.6 (0.83) 9.0 (0.78) 10.1 (0.81) -15.9 8.9 (0.778) 3.6 (0.448) -6.8

τ⊥/ns

τ|/ns c

0.078 0.069 0.042c 0.031 14.9 0.91 7.12 4.0

βD/deg

Rex/s-1

c

0.078 0.069 0.042c 0.006 0.013 0.0 0.004 0.021

4.4 3.6 12.2d 21.4

1.8 4.0

χ2 2.0 0.6 8.8 5.8 2.4 5.1 0.9 2

a Results of applying SRLS to the combined 14.1 and 18.8 T, 302 K, data of selected residues of E. coli adenylate kinase. The isotropic global diffusion correlation time was 14.9 ns. The underlined figure in row 5 indicates that τ⊥ was fixed at 14.9 ns in the respective calculation. The parameter values of rNH(eff) ) 1.015 Å,65 ∆σ ) -169 ppm,66 and βC ) 17° 67 were used. The local ordering and local diffusion frame were taken the same. The calculations were carried out on an HP computer equipped with an Intel 2.7 GHz Dual Core CPU and 4 GB of RAM. b The SRLS calculations vary c02, whereas the MF calculations vary S2. The complementary quantity has been calculated in each axial potential case. c MF data fitting yields τe, which in the limit of high axial local potentials agrees with the renormalized correlation time, τren ) 2t/c20.29 The data which appear in the table are τ values obtained from τe using this formula. These values should be compared with τ from SRLS. d The angle βD is derived from the MF parameters Sf2 according to Sf2 ) (1.5 cos2 βD - 0.5)2.12

15

N CSA frames13 (Figure 1) and (2) possible deviations in MF from the single-decay approximation for the local motion.13 As pointed out above, in the presence of local motions carrying out this frame transformation requires (besides C00(t)) the time correlation functions C11(t) and C22(t), which do not exist in MF; this is not always appreciated.70 By the same token, the MF formula for transverse 15N-1H dipolar/15N CSA crosscorrelated relaxation is only applicable when the local motion is practically in the extreme motional narrowing limit.53 In contrast, SRLS can treat cross-correlated relaxation in general over the entire parameter range relevant for proteins. It is therefore recommended to use SRLS instead of MF also in axial potential scenarios. Example 2: Rigid Residue with Rhombic Local Ordering. Residue 209 is also a rigid residue. In this case, the MF calculation did not pass the goodness-of-fit (GOF) criteria of the program DYNAMICS.68 The best results generated by this program, obtained with model 3 MF, are S2 ) 0.78 (c02 ) 9), τe ) 9.3 ps (τ ) 41.7 ps), and Rex ) 4.35 s-1. The χ2 value is 8.8, which is higher than the relevant threshold of 7.81 (Table 39 of ref 69). By using as variables c02, D2 (isotropic), and Rex, we obtained with SRLS c02 ) 10.1 ((S02)2 ) 0.805), τ ) 30.9 ps, and Rex ) 3.56 s-1. The χ2 value is 5.8, which is below the relevant threshold. These results are shown in Table 6, rows 3 and 4. As pointed out above, Rex ) can absorb the rhombicity of the local potential. With this in mind, we set βD ) 101.3° and RD ) 90° (as appropriate for Ci-1R - CRi considered to be the main local ordering/local diffusion axis)13-16,59 and allowed c22 to vary instead of Rex. To obtain good statistics and effective convergence, we had to set D2,⊥ ) D1 and allow D2,| to vary. The results of this calculation are shown in Table 6, row 5. The potential coefficients are c02 ) -15.9 and c22 ) -3.4. The Cartesian tensor components calculated from these coefficients (e.g., see ref 13) are Sxx ) -0.4010, Syy ) +0.8736, and Szz ) -0.4726. Viewed in the context of RD ) 90° and βD ) 101.3°, the fact that Syy ) +0.8736 means that relatively high ordering prevails along the Ci-1R - CiR axis. The anisotropy of the local ordering is evaluated by (Sxx - Syy)/Szz ) 0.082. The value of D2,| ) 1.3 × 1010 s-1 corresponds to τ ) 12.8 ps. This represents fast fluctuations of the N-H bond. Fixing D2,⊥ to be equal to D1 ) 1.117 × 107 s-1 (14.9 ns) means that backbone motions, which are expected (based on geometric considerations) to be associated with D2,⊥, will not be detected. Within the scope of combined two-field data, the rigid N-H bond of residue 209 does not sense backbone motions. However,

it makes possible determining in the context of realistic local geometry the magnitude and symmetry of the local ordering, the form of the potential in terms of which the order parameters are defined, and the rate of the local N-H fluctuations. Note that the search for a set of variables which differs from the simple scenario of axial potential and main local ordering/ local diffusion axis along the N-H bond was motivated by physical considerations. If only statistical criteria were considered, we would have accepted the results shown in row 4 of Table 6. At present, the local geometry is fixed at RD ) 90° and βD ) 101.3°, and D2,⊥ ) D1 is fixed. The analysis of combined three-field data, concerted analysis of temperaturedependent data, or combined analysis of several probes with their equilibrium orientation lying within the peptide plane, to be pursued in future work, might allow varying D2,⊥. The time required to complete the calculation illustrated in Table 6, row 5, was approximately 1 h. The local potential is high (c02 ) -15.9), the local ordering is high (Syy ) +0.8736), and the time scale separation between D2| and D1 is large (0.00086). The Lmax value required was 24. Higher potentials do not require much larger Lmax values, and the time scale separation has already reached a limiting value for which a robust fitting calculation should stop. Therefore this example might be considered to be representative of a long fitting calculation. Example 3: Flexible Residue with Rhombic Local Ordering. Residue 46 of AKeco is located in the mobile domain AMPbd. The program DYNAMICS68 selected model 7, but the calculation did not pass the GOF criteria. The best results were S2 ) 0.778 (c02 ) 8.9), Sf2 ) 0.87 (corresponding to βD ) 12.2°, according to Sf2 ) (1.5 cos2 βD - 0.5)2; see ref 12), τs ) 0.91 ns, τf ) 0.0 ps, Rex ) 1.8 s-1, and χ2 ) 5.1. Using SRLS with axial potentials and assuming that D2,| . D2,⊥, in analogy with τs . τf in MF, we obtained14 c02 ) 3.6 ((S02)2 ) 0.448), βD ) 21.4°, τ⊥ ) 1/(6D2,⊥) ) 7.12 ns, and τ| ) 1/(6D2,|) ) 0.004 ns. These results are shown in Table 6, rows 6 and 7. As reported previously,14 the SRLS and MF results differ significantly mainly because mode-coupling is not accounted for. The SRLS results are also problematic because c02 ) 3.6 represents too weak a potential inside a folded protein, and a 21.4° tilt from the N-H bond does not identify a structural element which can serve as a main local ordering/local diffusion axis.13 Row 8 of Table 6 shows the results obtained with C++OPPS by allowing the local potential to be rhombic and the local diffusion tensor, D2, to be axially symmetric. The angles RD and βD were set equal to 90° and 101.3°, respectively. The bestfit values of the potential coefficients are c02 ) -6.8 and c22 )

Interpreting NMR Spin Relaxation with SRLS -4.40, and the corresponding Cartesian tensor components are Sxx ) -0.4261, Syy ) +0.8758, and Szz ) -0.4497. The anisotropy is (Sxx - Syy)/Szz ) 0.027. D2,⊥ ) 4.15 × 107 s-1 corresponds to τ⊥ ) 4.0 ns, and D2,| ) 8.11 × 109 s-1, to τ| ) 20.6 ps. The physical picture is as follows. According to the local geometry τ⊥ ) 4.0 ns is associated with motion about the Ci-1R - CiR axis which is 3.7 times faster than the global tumbling of 14.9 ns. This is very likely to represent internal domain reorientation. The other local motional component detected, occurring with a rate of D2,| ) 8.11 × 109 s-1, represents fast fluctuations of the N-H bond. The value of τ⊥ is 3.7 times larger for the flexible N-H bond of residue 46 than for the rigid N-H bond of residue 209 (for which τ⊥ ) τm). The value of τ| is 1.6 times slower for the “flexible” residue. The “flexible” N-H bond of residue 46, and the rigid N-H bond of residue 209, are ordered locally around the Ci-1R CiR axis practically to the same extent: we obtained Syy ) +0.8758 for the former and Syy ) +0.8736 for the latter. On the other hand, the anisotropy of the local ordering, defined as (Sxx - Syy)/Szz, is nearly 3 times higher for the rigid site than for the flexible site. This is interesting new information. There are controversial views in the field of protein NMR on whether in solution proteins are well-structured or whether they prevail as conformational ensembles.71 The information emerging from the SRLS/C++OPPS analysis, which indicates that the protein might be altogether well-structured in solution, might help resolve these controversies. The 3D GAF model59,60 can also quantify the magnitude and anisotropy of the local ordering. However, it requires very fast local motions, the availability of MD trajectories, and it does not provide local potentials (which can be used to calculate thermodynamic quantities). Additional Experimental Relaxation Parameters. Extending the current version of C++OPPS to treat 15N-1H dipolar/ 15 N CSA cross-correlated relaxation, and other types of auto-, and cross-correlated relaxation scenarios, is straightforward (see Theoretical Background). Note that a given set of the jKK′(ω) functions obtained with 15N spin relaxation analysis can also be used to analyze any autocorrelated and cross-correlated relaxation data from probes with equilibrium orientations (local directors) attached rigidly to the peptide plane, where the equilibrium orientation of the N-H bond is assumed to lie. One only has to assemble the appropriate Jµν(ω) functions, write down the appropriate expressions for the experimentally measured relaxation parameters, and use the appropriate magnetic interactions. C++OPPS-Based 15N Spin Relaxation Analysis of an Entire Protein. Let us consider a large (by NMR standards) protein comprising 300 residues and assume that 10% of the residues are flexible. For the latter, one may use the parameter combination of example 3. This amounts to a total of 10.5 h. For the rigid residues, the parameter combination of example 2 is suggested as a paradigm. However, let us consider 1.5 h per residue for 270 residues, to account for the examination of a second parameter combination for 50% of the rigid residues; this amounts to 405 h. Altogether, approximately 17 days are required. This estimate is based on calculations carried out on a HP computer equipped with an Intel 2.7 GHz Dual Core CPU and 4GB RAM. On a Quad-Core i7 Extreme CPU with a 3.2 GHz clock speed 1600 MHz 8 MB cache and 24 GB of 1300 MHz CL6 RAM, the analysis of the 300 residue protein selected as an example will be completed in 4-5 days. A computer with these specifications is affordable. Utilization of the parallelized

J. Phys. Chem. B, Vol. 113, No. 41, 2009 13623 version of C++OPPS in the context of a computer cluster will further reduce the running time significantly. Problems Encountered in Some Cases and Prospects. The present data-fitting scheme features the publicly available MINPACK minimization package. This minimizer has not been adapted/optimized, and other minimizers have not yet been implemented/examined. With rigid N-H bonds, we encountered in some cases problems associated with the exit criteria of the MINPACK minimizer. In the context of SRLS/ESR, the Levenberg-Marquardt minimizer has been adapted/optimized successfully.72 We might be able to overcome the problem noted by optimizing the MINPACK minimizer, or another minimizer, for SRLS/NMR. Such efforts are underway. It is easy to fit the rigid N-H bond data with high axial potentials, fast isotropic local diffusion, and frequent inclusion of Rex contributions. This is similar to the results obtained with MF analyses. When anisotropic features are allowed for (rhombic potentials, axial local diffusion, and/or tilted local ordering/local diffusion and magnetic frames) the fitting process leads often to unphysical results. It is known that in the presence oflargelocalpotentialsoneshouldusetheFokker-Planck-Kramers (FPK) equation with both orientation and angular momentum explicitly included.5,36 This will allow the probe to engage in torsional oscillations in the potential well, expected on physical grounds, which in the overdamped Smoluchowski treatment are relaxed instantaneously. The full two-body FPK model is treated explicitly in ref 5. Efforts to implement this model for 15N spin relaxation in proteins are underway. The additional parameters required are moments of inertia which can be derived from 3D structures. 6. Conclusions SRLS is an advanced mesoscopic approach which accounts for the indispensable features required to extract properly from the experimental spin relaxation data the inherent information on structural dynamics of proteins. C++OPPS is a general and effective implementation of SRLS within the scope of a data fitting scheme. It can also be used for separate calculations of SRLS time correlation functions, generic spectral densities, measurable spectral densities, and NMR spin relaxation parameters. C++OPPS has already provided new physical insight. The underlying SRLS model features as physical parameters diffusion tensors for the global and local motions, a local ordering tensor, and Euler angles associated with the various frame transformations. It is practical to analyze the 15N spin relaxation data from a large protein using a modern portable computer within this physical framework. The following specific new physical information has been obtained, or can be obtained potentially. (1) In a geometric context where the main local ordering/local diffusion axis is parallel to the Ci-1R - CiR axis, rigid and flexible N-H bonds experience similar ordering strength, but different ordering asymmetry, which is higher for rigid N-H bonds. (2) Experimental data from rigid N-H bonds, which require a Rex contribution when axial local ordering is imposed, might be force-fitted. It is suggested to try and fit these data with rhombic ordering; the physical picture might be quite different. (3) Combining DITE, which calculates anisotropic global diffusion tensors from structure with SRLS within the scope of C++OPPS, provides a number of option for determining the global diffusion tensor, D1. One can calculate D1 with DITE, determine D1 with data fitting, or vary D1,| with D1,|/ D1,⊥ kept fixed. Comparison between the results obtained with these different strategies is expected to provide new information.

13624


(4) DITE allows for internal torsions. This capability, the implementation of which is straightforward, enables the determination of global diffusion tensors in the presence of internal domain reorientation which can be modeled as torsional motion. (5) SRLS/C++OPPS can treat diffusion in double-well potentials. In the limit of high barriers, this leads to jumps. Diffusion in such potentials is a very general way for treating internal motions in proteins. (6) 15N-1H residual dipolar couplings (RDCs) affected by picosecond to nanosecond motions are determined by the alignment tensor of the medium and by the same order parameters for local motion and local geometry as in 15N spin relaxation. Hence, concerted RDC/spin relaxation analysis, with obvious benefits, is possible. RDC and MF analyses are also integrated but one has to overcome with specific formalisms the fact that RDC requires rhombic global diffusion and rhombic local ordering tensors, whereas MF typically features (in this context) isotropic global diffusion and (implicitly) axial local ordering tensors. We propose straightforward integration of RDC and 15N spin relaxation analysis within the scope of SRLS; this is expected to provide new information. C++OPPS has to be enhanced to make this possible. (7) Identifying the necessity of treating 15N spin relaxation from rigid N-H bonds with a two-body FPK model instead of a two-body Smoluchowski model is an important new insight. Large parts of the protein structure are rigid. The overall picture of protein dynamics by NMR spin relaxation might change following the utilization of the FPK equation. Acknowledgment. We are grateful to Prof. Jack H. Freed for suggesting the implementation of the Fokker-Planck-Kramers equation to NMR spin relaxation in proteins. M.Z and A. P. acknowledge support provided by the Ministero dell’Istruzione, Universita` e Ricerca (MIUR), grant PRIN 2006 (2006033728), and by the Consorzio Interuniversitario per la Scienza e la Tecnologia dei Materiali (INSTM), grant PROMO 2007. E.M. acknowledges support from the Israel Science Foundation, grant no. 347/07; the Israel-US Binational Science Foundation, grant no. 2006050; the German-Israeli Science Foundation for Scientific Research and Development, grant no. 928-190.0/2006; and the Damadian Center for Magnetic Resonance Research at Bar-Ilan University. References and Notes (1) Arora, K.; Brooks, C. L., III Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 18496. (2) Maragakis, P.; Lindorf-Larssen, 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., III; Shaw, D. E. J. Phys. Chem. B 2008, 112, 6155. (3) Wong, V.; Case, D. A. J. Phys. Chem. B 2008, 112, 6013. (4) Trbovic, N.; Kim, B.; Friesner, R. A.; Palmer, A. G., III Proteins 2007, 71, 684. (5) Polimeno, A.; Freed, J. H. AdV. Chem. Phys. 1993, 83, 89. (6) Polimeno, A.; Freed, J. H. J. Phys. Chem. 1995, 99, 10995. (7) Sastry, V. S. S.; Polimeno, A.; Crepeau, R. H.; Freed, J. H. J. Chem. Phys. 1996, 107, 5753. (8) Earle, K. A.; Moscicki, J. K.; Polimeno, A.; Freed, J. H. J. Chem. Phys. 1997, 106, 9996. (9) Liang, Z.; Freed, J. H. J. Phys. Chem. B 1999, 103, 6384. (10) Liang, Z.; Lou, Y.; Freed, J. H.; Columbus, L.; Hubbell, W. L. J. Phys. Chem. B 2004, 108, 17649. (11) Liang, Z.; Freed, J. H.; Keyes, R. S.; Bobst, A. M. J. Phys. Chem. B 2000, 104, 5372. (12) Tugarinov, V.; Liang, Z.; Shapiro, Y. E.; Freed, J. H.; Meirovitch, E. J. Am. Chem. Soc. 2001, 123, 3055. (13) Meirovitch, E.; Shapiro, Y. E.; Polimeno, A.; Freed, J. H. J. Phys. Chem. A 2006, 110, 8366. (14) Tugarinov, V.; Shapiro, Y. E.; Liang, Z.; Freed, J. H.; Meirovitch, E. J. Mol. Biol. 2002, 315, 171. (15) Shapiro, Y. E.; Kahana, E.; Tugarinov, V.; Liang, Z.; Freed, J. H.; Meirovitch, E. Biochemistry 2002, 41, 6271.

Zerbetto et al. (16) Shapiro, Y. E.; Meirovitch, E. J. Phys. Chem. B 2006, 110, 11519. (17) Meirovitch, E.; Polimeno, A.; Freed, J. H. J. Phys. Chem. B 2006, 110, 20615. (18) Meirovitch, E.; Shapiro, Y. E.; Polimeno, A.; Freed, J. H. J. Phys. Chem. B 2007, 111, 12865. (19) Lipari, G.; Szabo, A. J. Am. Chem. Soc. 1982, 104, 4546. (20) Lipari, G.; Szabo, A. J. Am. Chem. Soc. 1982, 104, 4559. (21) Clore, G. M.; Szabo, A.; Bax, A.; Kay, L. E.; Driscoll, P. C.; Gronenborn, A. M. J. Am. Chem. Soc. 1990, 112, 4989. (22) Meirovitch, E.; Shapiro, Y. E.; Liang, Z.; Freed, J. H. J. Phys. Chem. B 2003, 107, 9898. (23) Choy, W.-Y.; Kay, L. E. J. Biomol. NMR 2003, 25, 325–333. (24) Prompers, J. J.; Bru¨schweiler, R. J. Am. Chem. Soc. 2001, 123, 7305. (25) Freed, J. H. J. Chem. Phys. 1977, 66, 4183. (26) Lin, W. J.; Freed, J. H. J. Phys. Chem. 1979, 83, 379–401. (27) Nordio, P. L.; P. Busolin, P. J. Chem. Phys. 1971, 55, 5485. (28) Polnaszek, C. F.; Bruno, G. V.; Freed, J. H. J. Chem. Phys. 1973, 58, 3185–3199. (29) Polnaszek, C. F.; Freed, J. H. J. Phys. Chem. 1975, 79, 2283. (30) Emsley, J. W.; Lindon, J. C. NMR Spectroscopy Using Liquid Crystal SolVents; Pergamon Press: Oxford, 1975. (31) Freed, H. J.; Nayeem, A.; Rananavare, S. B. The Molecular Dynamics of Liquid Crystals; Luckhurst, G. R., Veracini, C. A., Eds.; Kluwer Academic Publishers: The Netherlands, 1994; Chapter 12, pp 271-312. (32) Szabo, A. J. Chem. Phys. 1980, 72, 4620. (33) Lipari, G.; Szabo, A. Biophys. J. 1980, 30, 489. (34) Lipari, G.; Szabo, A. J. Chem. Phys. 1981, 75, 2971. (35) Polimeno, A.; Nordio, P. L.; Moro, G. J. Chem. Phys. Lett. 1987, 144, 357. (36) Moro, G. J.; Polimeno, A. Chem. Phys. 1989, 131, 281. (37) Moro, G. J.; Polimeno, A. Chem. Phys. Lett. 1992, 189, 133. (38) Wong, V.; Case, D. A.; Szabo, A Proc. Natl. Acad. Sci. U.S.A. 2009, 106, 11016. (39) Tombolato, F.; Ferrarini, A.; Freed, J. H. J. Phys. Chem. B 2006, 110, 26248. (40) Tombolato, F.; Ferrarini, A.; Freed, J. H. J. Phys. Chem. B 2006, 110, 26260. (41) Barone, V.; Polimeno, A. Phys. Chem. Chem. Phys. 2006, 8, 4609. (42) Barone, V.; Brustolon, M.; Cimino, P.; Polimeno, A.; Zerbetto, M.; Zoleo, A. J. Am. Chem. Soc. 2006, 128, 15865. (43) Zerbetto, M.; Carlotto, S.; Polimeno, A.; Corvaja, C.; Franco, L.; Toniolo, C.; Formaggio, F.; Barone, V.; Cimino, P. J. Phys. Chem. B 2007, 111, 2668. (44) Carlotto, S.; Cimino, P.; Zerbetto, M.; Franco, L.; Corvaja, C.; Crisma, M.; Formaggio, F.; Toniolo, C.; Polimeno, A.; Barone, V. J. Am. Chem. Soc. 2007, 129, 11248. (45) Barone, V.; Polimeno, A. Chem. Soc. ReV. 2007, 36, 1724. (46) Sezer, D.; Freed, J. H.; Roux, B. J. Phys. Chem. B 2008, 112, 5755. (47) Sezer, D.; Freed, J. H.; Roux, B. J. Chem. Phys. 2008, 128, 165106. (48) Sezer, D.; Freed, J. H.; Roux, B. J. Phys. Chem. B 2008, 112, 11014. (49) Barone, V.; Zerbetto, M.; Polimeno, A. J. Comput. Chem. 2009, 30, 2. (50) Meirovitch, E.; Shapiro, Yu.E.; Tugarinov, V.; Liang, Z.; Freed, J. H. J. Phys. Chem. B 2003, 107, 9883–9897. (51) Moro, G. J.; Freed, J. H. J. Chem. Phys. 1981, 74, 3757. (52) Moro, G. J.; Freed, J. H. Large-Scale EigenValue Problems, Mathematical Studies Series; Cullum, J. A., Willough, R. K., Eds.; Elsevier: New York, 1986. (53) Tjandra, N.; Szabo, A.; Bax, A. J. Am. Chem. Soc. 1996, 118, 6986. (54) Abragam, A. The Principles of Nuclear Magnetism; Oxford University Press (Clarendon): London, 1961. (55) Peng, J. W.; Wagner, G. Methods in Enzymology; James, T. L., Oppenheimer, N. J., Eds.; Academic Press: New York, 1994; Vol. 239, pp 563-595. (56) Kay, L. E.; Torchia, D. A.; Bax, A. Biochemistry 1989, 28, 8972– 8979. (57) Barnes, J. P.; Freed, J. H Biophys. J. 1998, 75, 2532. (58) Osborne, M. J.; Wright, P. E. J. Biomol. NMR 2001, 19, 209. (59) Bremi, T.; Bru¨schweiler, R. J. Am. Chem. Soc. 1997, 119, 6672. (60) Lienin, S. F.; Bremi, T.; Brutscher, B.; Bruschweiler, R.; Ernst, R. R. J. Am. Chem. Soc. 1998, 120, 9870. (61) Ulmer, T. S.; Ramirez, B. E.; Delaglio, F.; Bax, A. J. Am. Chem. Soc. 2003, 125, 9179. (62) Fushman, D.; Ohlenschlager, O.; Ruterjans, H. J. Biomol. Struct. Dynamics 1994, 11, 1377. (63) Pawley, N. H.; Wang, C.; Koide, S; Nicholson, L. K. J. Biomol. NMR 2001, 20, 149–165. (64) Shapiro, Yu. E.; Kahana, E.; Meirovitch, E. J. Phys. Chem. B 2009, 113, 12050. (65) Yao, L.; Vogeli, B.; Ying, J.; Bax, A. J. Am. Chem. Soc. 2008, 130, 16518.

Interpreting NMR Spin Relaxation with SRLS (66) Damberg, P.; Jarvet, J.; Graslund, A. J. Am. Chem. Soc. 2004, 127, 1995–2005. (67) Fushman, D.; Cowburn, D. J. Biomol. NMR. 1999, 13, 139. (68) Fushman, D.; Cahill, S.; Cowburn, D. J. Mol. Biol. 1997, 266, 173. (69) Spiegel, M. R.; Liu, J. Mathematical Handbook of Formulas and Tables, Schaum’s Outline Series, 2nd ed.; McGraw-Hill: New York, 1999; p 266.

J. Phys. Chem. B, Vol. 113, No. 41, 2009 13625 (70) Chatfield, D. C.; Szabo, A.; Brooks, B. R. J. Am. Chem. Soc. 1998, 120, 5301. (71) Clore, G. M.; Schwieters, C. D. J. Mol. Biol. 2006, 355, 879. (72) Budil, D. E.; Lee, S; Saxena, S.; Freed, J. H. J. Magn. Res. A 1996, 120, 155.

JP9046819

Computational Tool for Interpreting NMR Spin

Recommend Documents