Synergism between x-ray crystallography and NMR residual dipolar couplings in characterizing protein dynamics

The study of protein internal motions at the atomic level has been a topic of high biophysical interest for many decades. These motions are key to the entropic contribution to the free energy function that is of fundamental importance for characterizing protein–protein and protein–ligand interactions.^1,2 An improved view of protein dynamics may also shed new light on the age-old question of induced fit vs conformational capture.^3–5 

Advances in single particle cryo-electron microscopy have resulted in atomic resolution protein structures, with the heterogeneity in observed structures providing evidence for domain motions that can be of high functional importance.^6,7 However, at present the atomic details of protein motions and their time scales are more commonly analyzed by protein x-ray crystallography and NMR spectroscopy. In this perspective, we focus on recent innovations in these fields and discuss how their concerted use may lead to further improvements in characterizing the atomic details of molecular motions that drive protein function.

Temperature or B-factors contain information on the spatial distribution of dynamic fluctuations experienced by a protein,⁸ but their quantitative interpretation in structural terms is complex.^9–11 Using molecular dynamics to generate a “protein crystal” perturbed by known multimodal, anharmonic motions, Garcia et al. showed that experimental B-factors can substantially underestimate the true mean square displacements.¹² An analogous, in-depth analysis of simulated data showed that even at high crystallographic resolution, fitted B-factors of some well-resolved atoms can be underestimated by as much as sixfold.¹⁰ On the other hand, for proteins that diffract to very high crystallographic resolution, electron density often can be confidently interpreted in terms of distinct conformations that are populated in the crystalline state,^13,14 an observation sometimes referred to as polysterism.¹⁵ Such structural heterogeneities usually cluster in small regions, typically involving less than half a dozen residues, and the details and relative populations of these conformers vary with temperature.¹⁶ Residues that exhibit such structural heterogeneities in the crystalline state nearly always give rise to single resonances in solution NMR spectra, recorded at room temperature, indicative of conformational averaging between these substates that is rapid on the NMR timescale of milliseconds.^17,18 The functional importance of such dynamic interconversions has been well recognized in enzymatic mechanisms,^15,17,19 and is at the basis of Koshland's induced fit hypothesis which postulated that proteins can switch to alternate discrete conformational states upon ligand binding.³ 

For protein crystals that diffract to high resolution, i.e., better than ∼2-Å, ensemble refinement methods of increasing sophistication have been introduced over the past 30 years.^20–23 Rather than locally fitting a weighted average of two discrete conformers to the crystallographic data, these methods generate many protein conformers that are sampled from molecular dynamics (MD) trajectories, with these trajectories restrained to yield time-averaged agreement with the x-ray diffraction data.²⁰ Such ensemble refinement was shown to be prone to over-fitting of the experimental data, as exemplified by worse cross-validation statistics in terms of the crystallographic free R factor (R_free)²⁴ during early applications of this method.^25,26 However, the introduction of suitable adjustments in both the MD protocol, improved MD force fields, and deformable elastic net (DEN) restraints, mitigate over-fitting and can result in modest improvements in R_free.^21–23 An intrinsic advantage of ensemble refinement, over the fitting of multiple discrete local conformers, results from the correlation seen in dynamic excursions in different regions of a protein, as expected for example, during allosteric conformational switching. However, it is important to bear in mind that the correlations between motions depend upon the empirical force field used in the MD simulation, and do not result from the Bragg reflections per se. Scripts for carrying out crystallographic ensemble refinement using the very powerful and widely used Phenix software^27,28 have become available, making the approach accessible to non-expert users with sufficiently powerful computational resources.

NMR spectroscopy has long been used for probing the rates and amplitudes of dynamic excursions from the time-averaged structure by measurement of nuclear spin relaxation rates. Motions that are faster than the rotational correlation time of the protein in solution are typically probed by measurement of longitudinal (R₁) and transverse (R₂) relaxation rates of ¹⁵N or ¹³C nuclei, combined with the magnitude of the heteronuclear nuclear Overhauser enhancement (NOE).^29,30 The rotational correlation time scales approximately linearly with the molecular volume or mass of the protein, ranging from 4.2 ns for ubiquitin (8.5 kDa) at room temperature, to ca 32 ns for hemoglobin (64.5 kDa), and therefore only rather fast motions can be quantitatively probed by such relaxation measurements. The type of motion associated with these relaxation parameters usually cannot be extracted from such data and is simply presented as an effective generalized order parameter, S², that reports on the angular amplitude of the excursions of corresponding bond vectors in a model-free manner.^31–33 Recent innovative work, relying on transient binding of a protein to slowly tumbling nanoparticles, can extend the timescale over which such dynamic excursions can be probed to values beyond the rotational correlation time of the free protein.^34,35 Unfortunately, the requirement that the time-weighted average rotational correlation time of the free and nanoparticle-bound protein remains shorter than ca 30 ns, the upper limit for which high-resolution amide resonances can be readily observed at adequate sensitivity when using TROSY methods, limits this method to relatively small proteins. Sidechain motions on time scales longer than 30 ns can be studied by the nanoparticle-assisted spin relaxation (NASR) method when using the more favorable methyl group relaxation properties to probe the system.³⁶ 

Internal motions on time scales much slower than the rotational tumbling, in the μs-ms regime, are often reflected in chemical exchange contributions, R_ex, to the transverse relaxation rates, R₂. R_ex measurement exploits the difference in chemical shifts of the conformationally exchanging substates. For the common situation where the exchange is faster than the difference in chemical shifts (expressed in Hz), the R_ex contribution to the transverse relaxation rate scales with the square of this shift difference, and inversely with the exchange rate.³⁷ In addition to providing precise values of the exchange rate, such measurements also can yield chemical shifts of the substates involved, which contain information regarding the structure of the substate. Although, in practice, it often proves difficult to extract the chemical shifts for more than a single unknown substate, the R_ex data can be invaluable for validating models of the exchange process that gives rise to the exchange contributions.^38–40 

If precise interproton distances are derived from highly quantitative NOE data, using so-called “exact NOEs” or eNOEs,⁴¹ regions in a protein structure can sometimes be identified that cannot simultaneously satisfy all precisely measured distances, indicative of switching between conformational substates.⁴² In such cases, multi-conformer ensembles can be derived that show much improved agreement with the NMR data, albeit using an N-fold increase in the degrees of freedom when fitting to an ensemble of N conformers. Like for crystallographic ensemble refinement, stringent cross-validation criteria are, therefore, needed to prevent over-fitting of the data.⁴³ As will be discussed below, residual dipolar couplings (RDCs) provide an unambiguous set of parameters to validate structural models. For example, RDCs that were not used to derive an NMR structure of a small GB3 domain agreed considerably better with a single-conformer eNOE-derived model than with the traditionally NOE-derived structure, validating the power of the eNOE method. However, no significant further improvement was obtained when the structure was refined as a multi-conformer model.⁴⁴ 

Multi-conformer NMR ensembles also have been generated by using the relaxation-derived backbone order parameters as restraints during MD trajectory calculations.⁴⁵ This dynamic ensemble refinement (DER) method was demonstrated for ubiquitin, showing considerable conformational heterogeneity throughout the structure, much larger than coordinate uncertainties in either the deposited NMR structure,⁴⁶ or calculated from x-ray data by using RAPPER.⁴⁷ However, it is important to note that the original NMR structures of ubiquitin represent best fits to the time-averaged coordinates, whereas the DER method aims to depict the ensemble of structures sampled over time. The DER ensemble showed remarkably good agreement with both scalar couplings and residual dipolar couplings that were not used as restraints during the MD calculation, supporting the methodology. However, we note that S² order parameters, derived from NMR relaxation data, are very sensitive to the N–H bond length used. Due to the low mass of hydrogen, N–H vectors are subject to ultrafast zero-point liberations of rather large amplitudes that decrease the effective ¹⁵N–¹H dipolar coupling. Therefore, N–H bond lengths should be artificially lengthened to 1.04 Å when interpreting ¹⁵N relaxation in terms of motion of the N and C backbone atoms.⁴⁸ The large amplitude dynamics seen in the DER ubiquitin ensemble (PDB entry 1XQQ), forced by imposing compliance with the ¹⁵N S² values that were derived for a 1.02-Å NH-bond length, resulted in conformers that transiently lack H-bonds (e.g., I30, Q40, L50, and L56), whereas very slow hydrogen exchange⁴⁹ is consistent with their uniform presence in all 190 PDB x-ray structures solved at a resolution ≤2.0 Å. This observation, therefore, suggests that, despite good validation statistics, the DER ensemble representation may be too dynamic.

RDCs contain highly precise information on the time-averaged orientation of bond vectors in the molecular frame of a protein. Measurement of RDCs requires weak alignment of the protein relative to the external magnetic field, either by its own inherent magnetic susceptibility anisotropy⁵⁰ or by the application of an external alignment agent,⁵¹ most commonly a suspension of filamentous phage⁵² or anisotropically compressed hydrogel.^53,54

Dipolar coupling between two nuclear spins is described by a rank-2 tensorial interaction, and the averaged RDC for a bond vector that rapidly switches between multiple orientations in the molecular frame therefore will deviate from an RDC calculated for its time-averaged orientation.^55–57 Depending on the precision at which the experimental RDCs can be measured, such deviations only become detectable for root mean square angular excursions that are larger than ca 20°–30°. With the availability of a large set of very precisely measured RDCs and NOEs, ensemble refinement of ubiquitin yielded a small, cross-validated improvement for a two-conformer ensemble over a single conformer, but with local outliers that appeared to be artifactual.⁵⁸ 

Alternatively, if couplings can be measured under five orthogonal protein alignments (in rank-2 tensor space), these couplings allow the generation of NMR ensemble models that reflect the amplitudes and directions of bond vector fluctuations.^55,59,60 A landmark study that applied this technology to ubiquitin found evidence for large amplitude backbone dynamics on a microsecond timescale,⁵⁹ termed “recognition dynamics.”⁵⁹ The multi-alignment approach for extracting dynamics hinges on the inability to satisfy the five experimental RDCs, measured for a given bond vector under the five different alignment conditions, by a single vector orientation. However, each RDC also contains a random measurement error, and perfect fitting of five erroneous RDCs is possible due to the additional degrees of freedom that are available when representing the vector by an ensemble, even if its orientation is static. This problem of introducing false motion is aggravated when the alignment orientations are not sufficiently orthogonal to one another, in which case small measurement errors in the RDCs require a wider distribution of bond vectors, i.e., larger amplitudes of motion, to satisfy the data. In practice, application of the approach, therefore, requires RDCs that have been measured at very high precision, in particular when the five protein alignments are insufficiently orthogonal to one another. For ubiquitin, RDCs subsequently measured under a protein orientation that was poorly represented in the original study yielded better fits for most of the protein when using a single, static model rather than the ensemble representation,⁶¹ suggesting that microsecond motions pertain to a much smaller fraction of the protein than originally reported. We also note that, in practice, five linearly independent alignment orientations are difficult to achieve unless relying on paramagnetic tagging,⁶² which introduces other challenges.

Instead, as discussed below, even if measured in only a single alignment medium, the highly precise structural information contained in RDCs can be used to judge the validity and quality of both x-ray ensemble models and conventionally refined structures. Provided that sufficiently many x-ray structures are available for a given protein, RDCs can also be used to select an ensemble representation of such conventional low-energy structures that agrees better with RDC data than any individual structure.

Clearly, NMR relaxation data and x-ray crystal structures are highly complementary: the x-ray data yield precise coordinates for low-energy states, with the relaxation data reporting on the rates at which they interconvert if populated in solution and provided that this exchange occurs on a NMR-accessible timescale.⁶³ Quantitative information on the number of conformers and their relative populations is often difficult or impossible to extract from such relaxation data. However, as discussed below and recently demonstrated for the SARS-CoV-2 main protease (M^pro), RDCs hold strong potential to provide this missing information.⁶⁴ 

Agreement between RDCs and protein atomic coordinates depends primarily on four factors: (1) The precision at which RDCs can be measured relative to the range spanned by the RDCs. (2) The precision at which the atomic coordinates are known, which is correlated with the crystallographic resolution and the R_free at which the structure was solved.^65,66 (3) The accuracy at which hydrogen atom positions can be added to an x-ray structure, which can be challenging for amide groups that sometimes deviate substantially from idealized geometry. (4) True differences between a protein structure in the crystalline and solution states. The latter are often attributed to intermolecular contacts and lattice-packing effects, but well-defined differences between local conformations can also be observed between high-resolution structures determined in the same space group. An example of the latter has been discussed for the orientation of the Thr-198 peptide group in the inter-domain linker adjacent to the P5 ligand-binding pocket of numerous M^pro x-ray structures, all in the C₁₂₁ space group.⁶⁷ 

Fits of M^pro RDCs to each of 350 x-ray structures available in the PDB showed good agreement for most of the α-helical C-terminal domain as well as part of the N-terminal β-sheet-rich domain, but considerably poorer fits for residues in the intervening half of the protein, including its active site and many of its substrate-binding pockets.⁶⁸ Despite the fact that this homodimeric enzyme is large by NMR standards (68 kDa), which adversely impacts the precision at which RDCs can be measured, the fit of ¹⁵N-¹H RDCs to structure was shown to be limited by the accuracy of the coordinates rather than RDC measurement error.⁶⁴ This protein, therefore, presents a good testing ground for evaluating the utility of RDCs to quantify internal dynamics. First, we focus on the application of RDCs to evaluate different x-ray ensemble refinement methods.

Recently, M^pro was studied by various crystallographic ensemble refinement procedures, with mixed results. Whereas the standard ensemble refinement procedure, without DEN restraints, applied to the temperature series of x-ray data (PDB entries 7MHF-7MHK) yielded ensembles (7MHL-7MHQ) that, with one exception (7MHL, 100 K), generated small improvements in R_free,⁶⁷ the dynamic distributions of these ensembles differed considerably.⁶⁴ 

Application of DEN restraints during ensemble refinement, starting from crystallographic data recorded at very high resolution (PDB entry 7K3T; 1.2 Å),⁶⁹ resulted in moderately (ca 10%) improved agreement with RDCs over the conventionally refined two-conformer model.⁶⁴ The results for M^pro obtained with the recently introduced ECHT ensemble refinement protocol²² have also been reported.⁶⁹ ECHT refinement breaks the disorder down to different size scales and provides a more intuitive view of the internal dynamics.

The ECHT-1 ensemble refinement protocol, applied to the PDB-REDO coordinates of the 7K3T M^pro structure, yielded a moderately improved correlation with the RDCs, in particular for the more dynamic fraction of residues. This improved correlation indicates that the directions of excursions seen in the ensemble agree with motions sensed by the RDCs, even though Q factors remain about 20%–40% worse than for highly refined small proteins such as ubiquitin,⁶¹ or the third IgG-binding domain of protein G.⁷⁰ Another indication that the x-ray ensemble refinement result remains imperfect is obtained when separately considering the most and least dynamic fraction of amides in the ensembles: if motions are correctly represented in the ensemble, separate RDC fits to the most and least disordered vectors in these ensembles should result in very similar D_a values. However, a 5%–10% increase in D_a value for the most disordered half of the residues indicated that the amplitudes of the dynamics assigned to their amides was considerably too large.

When an artificial “superensemble” was generated from all 350 M^pro PDB x-ray structures (distributed over eight space groups, different ligation states, and ranging from 1.2 to 2.85-Å in resolution), considerably better agreement with RDCs was obtained than for any individual structure or any of the crystallographic ensemble refinements (Fig. 1).⁶⁴ Moreover, fits of RDCs to the 50% of residues with the highest disorder (H50%) and lowest disorder (L50%) in this superensemble yielded nearly the same D_a value,⁶⁴ indicating that the amplitudes of excursions in such an ensemble quantitatively agreed with solution RDCs. For the L50% fraction of residues, the improvement in the ensemble fit over the fit to a single model, ⟨Xray^*⟩, where each vector adopts the ensemble-averaged vector orientation, is very small (Table I). This result suggests that for these least dynamic residues the improved ensemble RDC fit over the fit to any individual structure is dominated by decreased uncertainty in vector orientations after averaging over 350 different models. In contrast, for the dynamic H50% fraction of residues, the ensemble fit improves by 14% over fitting the single, ensemble-averaged model (Table I), indicating that it is the ensemble distribution, not simply the averaging of uncertainties in orientation, that is responsible for the better fit. For the ensemble fit to the 350 x-ray structures, the Q factor for the most dynamic, H50%, residues is only slightly worse than for the L50% half (Table I). This result indicates that this “random collection” of x-ray structures, the vast majority of which were determined by molecular replacement methods but also representing different ligation states and eight space groups, represents a remarkably good quantitative representation of the protein in its unligated solution state that was used for the RDC measurements.

Ubiquitin has also been studied extensively by both NMR and x-ray crystallography. Coordinates for more than 350 ubiquitin chains are available in the PDB, distributed across 31 space groups. Below, we focus on 190 ubiquitin chains that were solved at a resolution ≤2.0 Å. Many of these concern poly-ubiquitin constructs, frequently complexed to other proteins. A large set of highly precise RDCs in solution, measured under four different alignment orientations, is available for free ubiquitin and can be used to evaluate the backbone dynamics of this protein, which often has been used as a benchmark for dynamics analysis by NMR relaxation methods.^{34,45,59,71–73}

A plot of the RDC Q-factor against crystallographic resolution confirms that, on average, higher resolution crystal structures agree better with RDCs [Fig. 2(a)]. However, many of the structures solved at above 2-Å resolutions also yield good fits. As was observed for M^pro, when considering these chains as a dynamic ensemble with equal weights for all members, the fit to the RDCs improves to well beyond what is obtained for any of the individual structures [dashed line in Fig. 2(a)]. A best fit superposition of these chains shows a fairly tight bundle, with a backbone C^α coordinate root mean square deviation (rmsd) of 0.65 Å for residues Q2-L71 [Fig. 2(b)].

The high-order parameters seen for most of the protein do not exclude the presence of transient large amplitude motions, such as associated with flipping of aromatic rings or transient unfolding, whose rates are simply too slow to contribute measurably to the ensemble-averaged S² values in relaxation analysis. Similarly, the corresponding high energy states cannot be crystallized and are therefore not observed in the PDB.

When multiple x-ray structures for the same protein are present in the PDB, these constitute a rather arbitrary subset of the conformational space sampled by the protein, with each conformer impacted by intermolecular contacts in the crystalline lattice or in molecular complexes that they may be part of. Even if they span the conformational space sampled by the protein in solution, it is highly unlikely that the population of conformers in solution quantitatively corresponds to those deposited in the PDB. Therefore, the above RDC fitting procedure, carried out by assuming equal weights for each of the PDB depositions, is likely to deviate from the true situation and therefore to be sub-optimal. On the other hand, optimizing the weights for each of N PDB depositions introduces N-1 adjustable parameters in the RDC fit and can result in overfitting. This problem can be mitigated by cross validation, leaving out a small fraction of the RDCs used in the fit, while the remainder is used to optimize the weights and keeping fixed all previously optimized parameters related to alignment tensor orientation and magnitude. Computationally, such cross-validated optimization can become burdensome and improved procedures for ensemble optimization require further development. Efforts at optimizing the weights of the 350 M^pro PDB structures, based on a simulated annealing protocol, improved the RDC fit by about 5% over weighing all structures equally.⁶⁴ 

Another question one may ask is: What is the minimum number of PDB structures needed to fit the RDCs satisfactorily? This latter question becomes of particular interest when focusing on small, local regions with increased dynamics. Can the RDCs of such local regions be satisfied by using a linear combination of just a few structural models? If so, this may yield intuitive insights regarding such motions.

Application of this procedure to the M^pro linker, A191-T199, bordering the P5 substrate binding pocket, shows that a satisfactory fit to the data cannot be obtained for any single x-ray structure [Fig. 3(a)]. A fit to all 350 conformers, evenly weighted, agrees better with the RDC data, but not as well as an ensemble of three loop conformers with populations of 20%, 16%, and 64% [Fig. 3(b)]. These three conformers, which represent the minimum needed to satisfy the RDCs to within their experimental uncertainties, exhibit fairly large differences in local structure, providing insights into the type of dynamic excursions at this functional linker site [Fig. 3(b)].

The high degree of complementarity between RDCs and x-ray structures holds strong potential to provide improved, quantitative insights into structural dynamics. To date, the measurement of RDCs for characterizing dynamics has remained a niche, and their analysis in terms of motions remains challenging in the absence of crystallographic data. However, RDC data can rapidly and unambiguously be measured and assigned, and in conjunction with crystallographic data offer new opportunities to characterize protein motions in greater detail.

Recent improvements in crystallographic ensemble refinement, in particular the use of deformable-elastic-net-based restraints, appears to make the procedure more robust.²² On the other hand, x-ray ensemble refinements^67,69 that aimed to describe the internal motions in SARS-CoV-2 M^pro only provided a moderate improvement in agreement with RDCs,⁶⁴ measured by NMR for the apo-enzyme in solution. That analysis indicated that the motional amplitudes of the most dynamic regions in the ensemble refinements were substantially overestimated and inconsistent with the RDCs. These findings, therefore, suggest that RDCs can be used as valuable restraints to further optimize the crystallographic ensemble refinement method. However, analysis of M^pro also highlighted areas that will remain impossible to capture by crystallographic ensemble refinement. For example, substantial structural differences adjacent to the P5 binding pocket of M^pro were seen among different crystal structures, far outside uncertainties related to their electron densities.⁶⁷ Correspondingly, RDCs measured for T198 in this P5 pocket strongly disagree with the highest resolution x-ray structures and one of their ensemble refinements, but agree well with a multi-conformer structure that mostly populates an alternate state seen in several of the other high-resolution structures.

We observed that for two proteins that have been very extensively studied by both x-ray crystallography and NMR, M^pro and ubiquitin, improved RDC agreement is obtained for pseudo-ensembles of all PDB x-ray structures over fits to any individual structure. Modest improvement in RDC fits for regions of well-defined secondary structure with below average C^α coordinate RMSD appears to be dominated by a more accurate orientation of the bond vectors for the ensemble average over any individual structure, rather than by dynamic effects. However, the substantially larger improvement in pseudo-ensemble RDC agreement seen for the loop regions,⁶⁴ even when compared to ensemble-averaged vector orientations (Table I), indicates that for these more dynamic regions the RDCs contain valuable information on the distribution of conformations sampled by the protein. The improved RDC fits obtained for large ensembles of M^pro and ubiquitin x-ray structures suggest that these ensembles of x-ray structures provide fairly realistic descriptions of the conformations sampled in solution. NMR relaxation data then can provide precise time scales for the rates at which such structures interchange, whereas RDCs may aid in adjusting the relative populations.

For both ubiquitin and M^pro, the large diversity of PDB structures results from the work of dozens of laboratories over multiple years. It appears likely, however, that a much smaller number of structures may suffice to represent the sampled conformational space. At that point, RDCs can be used to judge how representative such a “micro-ensemble” is for its solution structure, and weight factors for its members may be optimized to improve quantitative agreement with the data. Such work will represent an important step toward the ultimate goal of pinpointing the functionally relevant conformers and learning how they interconvert.

We thank Stephan Grzesiek, Dennis Torchia, and Charles Schwieters for helpful discussions. This work was supported by the Intramural Research Program of the National Institute of Diabetes and Digestive and Kidney Diseases (Grant No. DK029046).

The authors have no conflicts to disclose.

Yang Shen: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Software (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Ad Bax: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Software (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.