TRACT revisited: an algebraic solution for determining overall rotational correlation times from cross-correlated relaxation rates

Accurate rotational correlation times (τc) are critical for quantitative analysis of fast timescale NMR dynamics. As molecular weights increase, the classic derivation of τc using transverse and longitudinal relaxation rates becomes increasingly unsuitable due to the non-trivial contribution of remote dipole-dipole interactions to longitudinal relaxation. Derivations using cross-correlated relaxation experiments, such as TRACT, overcome these limitations but are erroneously calculated in 65% of the citing literature. Herein, we developed an algebraic solutions to the Goldman relationship that facilitate rapid, point-by-point calculations for straightforward identification of appropriate spectral regions where global tumbling is likely to be dominant. The rigid-body approximation of the Goldman relationship has been previously shown to underestimate TRACT-based rotational correlation time estimates. This motivated us to develop a second algebraic solution that employs a simplified model-free spectral density function including an order parameter term that could, in principle, be set to an average backbone S2 ≈ 0.9 to further improve the accuracy of τc estimation. These solutions enabled us to explore the boundaries of the Goldman relationship as a function of the H-N internuclear distance ( r ), difference of the two principal components of the axially-symmetric 15N CSA tensor (ΔδN), and angle of the CSA tensor relative to the N-H bond vector (θ). We hope our algebraic solutions and analytical strategies will increase the accuracy and application of the TRACT experiment.


Introduction
A particle's rotational Brownian diffusion is characterized by the average time to rotate one radian, also known as the rotational correlation time (τ c ). It is related to the size and shape of a molecule, and in the case of a rigid, spherical particle, can be estimated from the Stokes-Einstein relation 1 . The rotational correlation time is frequently used in biophysics to gauge molecular aggregation and solvent viscosity. In protein NMR, rotational correlation time estimates are used to optimize interscan recycling delays, magnetization transfer delays in correlation experiments, and indirect dimension evolution times in multidimensional experiments 2 . Perhaps most significantly, τ c is the critical parameter for quantitative dynamics analyses, such as 'model-free' formalism, in which separation of overall and internal motion are required 3,4 .
The 1989 seminal work by Kay, Torchia and Bax 5 showed τ c can be estimated from 15 N longitudinal (R 1 ) and transverse (R 2 ) relaxation rates using equations established by Abragam 6 . These equations assume there is no rapid internal motion of the internuclear bond vector (i.e. motion faster than the rate of molecular tumbling), and that relaxation only results from i) dipole coupling with the covalently bonded nucleus and ii) the chemical shift anisotropy of the relaxing nucleus. In the original model free formalism, the amplitude and time parameter for this fast internal motion are referred to as the order parameter (S 2 ) and the effective correlation time (τ e ), respectively, where it is assumed that τ e ≪ τ c 3,4 . Note that τ m is frequently used in the literature to refer to overall molecular tumbling time, where 1/τ m = 1/τ e + 1/τ c . In this present work, unless stated otherwise, we assume fast motions on the time scale of τ e ≪ τ c do not exist and thus τ c is synonymous with τ m . In the absence of τ e , the S 2 order parameter essentially cancels out when R 2 is divided by R 1 5 . Thus, the R 2 /R 1 ratio can be used to estimate τ c by using the approximate formula in Eqn. 1, where v N is the 15 5 . This method of excluding spin systems from τ c estimates was subsequently automated by Clore et al 8 , who first suggested determination of the mean R 2 /R 1 ratio and exclusion of all values outside one standard deviation. Implicit in this model is that the tumbling is isotropic 5,8 ; Barbato et al 9 expanded application to anisotropic biomolecular systems. Anisotropic tumbling deviates R 2 and R 1 values, independent of chemical exchange or fast internal motions, but in opposing directions. Therefore, exchanging residues in which the R 2 /R 1 model does not apply can be detected by considering R 2 values that deviate more than one standard deviation from the mean, but without an associated decrease in R 1 relaxation 9 . Again, spin systems with 15 N{ 1 H}-NOE < 0.6 (indicating significant fast time scale motion) were excluded from τ c estimation. The remaining R 2 /R 1 values are then used to calculate τ c using Eqn 1.
Statistical selection of spin systems for τ c calculation is generally appropriate but does have the potential for mishandling by an inexperienced user. For example, a protein with significant regions undergoing chemical exchange will report high R 2 or R 2 /R 1 standard deviations from inclusion of unsuitable spin systems. It is also not hard to imagine systems with significant regions of high internal motions, but not complete random coil, which would skew the distribution of R 1 values and lead to the inclusion of spin systems where the R 2 /R 1 model is a poor approximation. Lastly, the R 2 /R 1 method does not account for the effects of dipolar couplings to remote (non-covalently bonded) protons on R 1 , which increase near exponentially as a function of molecular weight 10,11 ; although perdeuteration largely ameliorates this problem. Nonetheless, there are clear advantages to approaches that can estimate τ c from a measurable relaxation parameter that is insensitive to the effects of remote dipolar couplings and chemical exchange.
One such quantifiable phenomenon is transverse cross-correlated relaxation (CCR, η xy ). CCR results from the coordinated rotation of two nuclei in a magnetic field 12 and is primarily a function of dipole-dipole (DD) coupling, the chemical shift anisotropy (CSA) of the observed relaxing nucleus, and the molecular rotational correlation time. Measurements of CCR exploit the fact that the sign (±) of η xy depends on the spin states of the coupled nuclei. Cross-correlated relaxation only contributes to R 2 and R 1 in coupled systems because the opposing (opposite sign) contributions to relaxation cancel for decoupled spins 12 . 1 H-15 N and 1 H-13 C aromatic TROSY experiments exploit this property by only selecting signals from the spin state with relaxation interference 13,14 . Several methods have been developed to measure η xy in the 1 H-15 N spin system [15][16][17][18] where the most common approach is a pair of spectra that record the transverse relaxation rates of 15 N alpha (R α ) and beta (R β ) spin states. These rates are the sum of the auto-relaxation rate (R auto ), remote 1 H dipolar interactions (R D ), chemical exchange (R ex ), and η xy (Eqns. 2 and 3); where η xy is derived from subtraction of R α and R β 11,16,17 (Eqns. 4 and 5).
R β = R auto + R D + R ex + η xy R β − R α = 2η xy (4) η xy = R β − R α 2 (5) Using this method, the contributions to relaxation from remote protons and chemical exchange are cancelled and extraction of η xy is possible. Goldman 12 showed that η xy can be estimated given τ c , via the spectral density function, using Eqn. 6. η xy = pδ N 4J 0 + 3J ω N 3cos 2 θ − 1 (6) where: J ω = 2τ c 5 1 + τ c ω 2 (9) and: Goldman's relationship between η xy and τ c was first exploited experimentally by Lee and colleagues in the [ 15 N, 1 H]-TRACT (TROSY for rotational correlation times) pulse sequence 11 ; although, the manuscript does not explicitly detail how the Goldman relationship was solved. There are three important assumptions when applying the Goldman approach to TRACT data: i) the region analyzed is free of fast internal motions, which is rarely known a priori, ii) the system is a rigid body (i.e. S 2 ~ 1.0), and iii) the spin pair is isolated from remote DD/CSA relaxation interference. When pursuing our own TRACT analyses, we noted our τ c calculations were inconsistent with the original manuscript despite using identical physical and geometric constants.
Herein, we report an algebraic solution to the Goldman approach for straightforward calculation of τ c from measured R α and R β relaxation rate. Using this solution, we show that a numerical error in the original TRACT report has propagated into a high proportion of the citing literature. We use our algebraic solution to investigate complications from fast internal motions and propose analytical strategies to exclude unsuitable spin systems. The impact of order parameters motivated us to develop a second algebraic solution that includes S 2 as a parameter. We also noted that little attention has been paid to distributions of the difference of the two principal components of the axially-symmetric CSA tensor (Δδ N ), the CSA tensor angle relative to the internuclear bond vector (θ), and the internuclear distance (r). We show that the chosen value can have a non-negligible effect on τ c calculations but, as the relationship is near linear, a symmetrical random distribution around an average value would cancel out over many spin systems. Finally, we discuss how relaxation interference between remote 1 H DD and local CSA affects the observed CCR rate; a phenomenon that is independent of RD above but is generally negligible compared to other factors discussed in this paper.

An algebraic solution to the Goldman relationship
While validating our own numerical solution to the Goldman relationship, we noted a 6.6% overestimation of η xy in Lee et al 11 . Specifically, Figures 3 and 4 in their paper present τ c = 21 ns and 24 ns, respectively, from which we calculate η xy = 27.1 Hz and 30.9 Hz, respectively, using Eqns. 6-9 and the quoted physical constants. This is inconsistent with the reported R α and R β which yield η xy rates of (64-13)/2 = 25.5 Hz and (80-22)/2 = 29 Hz, respectively, using Eqns. 4-5. Hypothesizing that this discrepancy may be the result of poor numerical minimization, we generated an exact solution to Eqn. 6 with respect to τ c , given B0 and η xy .
We start by expanding the spectral density function (J) in Eqn. 6 and substituting η xy with (R β − R α )/2 to give Eqn. 10, where ω N is the 15 N Larmor frequency (in radians per second): The right-hand side of Eqn. 10 is a constant once the relaxation rates have been measured. We therefore replace this side with the symbol 'c'.  (12) where, an increasing interest in this methodology 10-15 years after its' original publication. We focused our analysis on manuscripts with sufficient data to validate calculations using our algebraic solution (Eq. 12). Table 1 shows that out of 120 manuscripts, three referenced the original TRACT paper without performing any cross-correlated relaxation experiments, six were methodological reviews, and seventy-eight used the TRACT methodology but did not provide enough experimental information to confirm their calculations. Thirty-three papers (30%) reported sufficient data to authenticate 65 total calculations.
We defined an error measure by dividing the published τ c by our algebraically determined value. For example, the original paper reported τ c = 21 and 24 ns, while we calculated values of 19.8 and 22.5 ns, giving error ratios of 1.061 and 1.067, respectively. There are several features of this analysis worthy of report (Fig. 1A) We also noted that researchers are more likely to overestimate the rotational correlation time, especially for low τ c values. A scatter plot of τ c versus error ratio demonstrates an inverse trend with highly erroneous values when τ c < 2 ns (Fig. 1B). Plotting τ c errors by year of publication underscores the persistence of these miscalculations in contemporary literature (Fig. 1C).

Evaluation and analysis using experimental data
While the CCR experiments eliminate complications from chemical exchange and remote dipolar couplings (Eqns. 2 and 3), the existence of fast internal motions cannot be established from these data alone. Spin systems possessing non-negligible τ e would artificially reduce τ c values using Goldman's relationship (Eqns. 6 and 12) 12 . This concern is especially relevant to the TRACT approach because spectra are often only collected in the directly-acquired 1 H dimension, and then analyzed by integration over a chosen 1 H N region (typically 1 H N δ > 8 ppm) to improve S/N. This leads to significant signal overlap, especially in the high molecular weight target proteins for which these experiments were designed 10,11 . Overlap itself can be mitigated by acquiring 2D versions with indirect 15 N evolution (although at a significant time expense), but confirmation of fast timescale motions still requires 15 N{ 1 H}-NOE data which are especially problematic in perdeuterated, high molecular weight systems 21 .
Our algebraic solution enables rapid point-by-point τ c calculation, which we used to explore the signal overlap problem using [U-15 N, 2 H]-OmpX prepared in DPC micelles (Fig. 2). The 1D 1 H N TROSY (i.e. N α spin state) spectrum at a relaxation delay = 1 ms is well dispersed with high intensity from 8.5 -8.0 ppm indicative of many overlapped spin systems and/or rapid local motions ( Fig. 2A). As expected, regions with weak signal intensity (e.g. 9.8 ppm ≥ 1 H δ ≤ 8.0 ppm) give wildly variable τ c estimates (Fig. 2B). There is sufficient signal intensity between 9.5 -8.8 ppm to calculate reasonable τ c values of 30-50 ns.
A significant drop in τ c is observed as data approaches the center of the amide region (~8.5 ppm; Fig. 2C), which reflects OmpX's unstructured loops 3 and 5 19 that resonate at random coil chemical shifts 22 . Together, this demonstrates that arbitrary selection of a region for integration, without knowledge of underlying dynamic processes, is problematic. Even estimations from relatively invariable regions, such as 9.5 -8.8 ppm, still possess high variance on a point-by-point basis. We applying a 200 point (5% of 4096 point amide region) sliding window as an optimal compromise for both sensitivity (via integration) and variability (region selection). The 200 points in the window are integrated and a τ c prediction generated (Fig. 2C). As illustrated in Figure 2, the sliding window approach narrows the standard deviation for straight-forward identification of regions with consistent τ c values. For example, calculations from 9.5 -9.17 ppm, with and without a sliding window, result in τ c of 37.92 ± 4.25 ns and 37.17 ± 0.59 ns, respectively (Fig. 2C,D). While commonly-applied assumptions about dispersed signals and global tumbling are useful when analyzing overlapped signals, point-by-point calculations coupled with the sliding window approach enable data-driven verification of consistent properties across the region of interest.

Rigid-body approximation as source of systematic error
In its current form, the Goldman relation (Eqn. 6) and our algebraic solution (Eqn. 12) do not account for bond motions. As well-ordered protein regions generally possess backbone 15 N H 0.85 ≤ S 2 ≤ 0.95 23,24 , we hypothesize the rigid-body assumption results in an underestimate of rotational correlation times estimated from TRACT experiments. We modified the spectral density function in Goldman's relationship (Eqn. 6) to include an order parameter (Eqn. 14) and solved for τ c as a function of B 0 , R α , R β , and S 2 (Eqn. 15). Note, we use a simplified form of the model-free spectral density function 3 that excludes fast timescale motions (τ e ) which TRACT data alone is insufficient to estimate. Therefore, Eqns. 14 and 15 are only applicable to spin systems where τ e has been established to be negligible through the process above.
J ω = 2τ c O We then modelled the effect of S 2 order parameters on τ c estimation using Eqn. 15 (Fig.  3). Relative to the rigid body assumption, typical backbone order parameters (0.85 ≤ S 2 ≤ 0.95) increase the rotational correlation time 5-25% depending on the CCR rate; this error increases to ≥45% at an S 2 ≈ 0.7. Plotting errors in τ c as a function of η xy for various order parameters reveals a pronounced biphasic character centered ≈ 4 Hz that becomes largely invariant at η xy > 10 ns (Fig. 3B). This simulation was performed with B0 = 800 MHz where η xy = 4 Hz gives a τ c of ~ 4 ns (assuming S 2 = 1.0). There is insufficient data to estimate S 2 from TRACT data alone; however, as shown, the rigid body assumption leads to significant underestimation of the rotational correlation time. This problem has been previously discussed by Wand and coworkers 25 who found TRACT underestimated τ c by an ≈ 20% for five different proteins, with errors ranging from 15-35%, when compared to more rigorous methods. Here we show that, in principle, these errors could be substantially reduced by inclusion of an order parameter to the Goldman relationship.

Potential sources of systematic error
The internuclear distance (r), difference of the two principal components of the axially symmetric CSA tensor (Δδ N ), and angle of the CSA tensor relative to the N-H bond vector (θ) are three additional parameters assumed constant in the equations above. These values are typically applied uniformly across the protein in 15 N relaxation analyses although they're well documented to be dependent on local structure 26,27 . Further, there are multiple commonly-employed values used throughout the literature, which are, themselves, interdependent and likely sources of systematic error 26,27  To explore the effect of each parameter on the Goldman relationship, we calculated rotational correlation times as a function of each parameter individually and plotted the percentage error relative to τ = 1.02 Å, θ = 17° and Δδ N = −160 ppm (Fig. 4). While τ and θ deviations generally effect rotational correlation time estimates within a ±5% error, variations of Δδ N can result in errors up to ±15% (Fig. 4).
It is important to note, however, that these errors approximate a mathematically-odd function centered around the elected value. That is, integration over many spin systems would cancel out small deviations around the average value for these three parameters, regardless of the parameter's chosen magnitude. Although, this characteristic would have little benefit in situations when the geometric constant is obviously inappropriate, such as using a helical CSA tensor to evaluate a beta-stranded protein.
Finally, we consider the contribution of remote dipole-dipole interference with local CSA on measured transverse cross-correlated relaxation rates. These interactions are distinct from the R D contributions in Eqns. 2 and 3 and do not cancel out when subtracting R β from R α . Generally, remote protons are not found within 2.2 Å of a given 1 H-15 N spin system, regardless of secondary structure, due to steric exclusion 34 . Given the τ −3 dependence of the dipole effect, protons at a distance of 2.2 Å contribute approximately 10 times less than protons at 1.02 Å to dipole relaxation. Liu and Prestegard have previously simulated the contributions of remote dipole/local CSA interference on measured CCR rates for the yARF1 protein 18 . They demonstrate that the average error in the CCR rate is 0.75% with an upper limit of ~3.5%{Liu, 2008 #17. In the perdeuterated case, which is essential for the study of proteins > 30 kDa, this effect would be even further minimized. Moreover, if the structure is known, the method described by Lui and Prestegard could be used to determine the error contribution.

Conclusion
All NMR methods for estimating the rotational correlation time depend on a number of assumptions concerning how molecules behave in solution. The most comprehensive, and time-consuming, approach involves the collection of multidimensional T 2 , T 1 and heteronuclear 15 N{ 1 H}-NOE experiments. Establishing which spin systems have suitable behavior is not trivial and difficult to automate; nonetheless, this strategy is independent of order parameters. This advantage is significant and enables site-specific τ c estimation. Two major drawbacks are the difficulties associated with applying these experiments to high molecular weight systems, and the significant effect of remote dipolar couplings on measured longitudinal relaxation rates. In an attempt to circumvent complications from chemical exchange and remote protons, Lee and coworkers developed the TRACT experiment to estimate rotational correlation times from cross-correlated relaxation rates using the Goldman relationship.
Herein, we developed two algebraic solutions to the Goldman relationship for accurate calculations assuming the rigid-body approximation or a specific order parameter. These solutions enabled us to explore the boundaries of the Goldman relationship without relying on numerical minimization, which is computationally slow and potentially inaccurate. However, as we have discussed in this paper, accurate analysis of TRACT data also requires careful consideration. First, there is no way to directly detect spin systems with fast internal motions that would undervalue τ c estimates. Second, experiments are frequently collected in a one-dimensional mode that is quite fast, but sensitive to signal overlap. The algebraic solutions facilitate rapid, point-by-point calculations for straightforward identification of appropriate spectral regions where global tumbling is likely to be dominant. Combining this approach with a sliding window simulates the advantages of integration while minimizing the inclusion of inappropriate spin systems. We also demonstrate that the rigid-body approximation can substantially underestimate TRACT-based rotational correlation time estimates. Our algebraic solution incorporates a simplified model-free spectral density function with order parameter that could, in principle, be set to an average backbone S 2 ≈ 0.9 to further improve the accuracy of τ c estimation. This has not been considered previously. Deviations in τ, θ and Δδ N contribute modest errors to τ c estimation, although these would be expected to cancel out over a large number of spin systems. We hope our algebraic solutions and analytical strategies will increase the accuracy and application of the TRACT experiment.