Quantum chemical ¹³C^α chemical shift calculations for protein NMR structure determination, refinement, and validation

Determination of a set of 10 protein conformations for YnzC (Set-bt).

A physics-based method, aimed at determining protein structures by using NOE-derived distance constraints together with observed ¹³C^α chemical shifts and those computed at the density functional level of theory, was applied here to determine 10 conformations (Set-bt) of the YnzC protein. The C-terminal ≈35-residues of the YnzC protein is intrinsically disordered in solution (15). Accordingly, calculations were carried out on 46-residue (YnzC[1–46]) and 52-residue (YnzC[1–52]) constructs of the 77-residue YnzC protein, including the first 46 and 52 residues, respectively, of the native sequence together with an additional 8-residue LEHHHHHH C-terminal affinity purification tag; the (His)₆ tag was omitted for quantum chemistry calculations, which were carried out only for the first 48 residues of YnzC[1–52]; i.e., including only ≈6 residues from the intrinsically-disordered C-terminal region. The determination of this structure was carried out in a “blind test” manner. Thus, the only information used (and provided) was the full set of NOE-derived distance constraints and the observed ¹³C^α chemical shifts. In other words, information about the atomic coordinates of structures previously solved by NMR spectroscopy (2JVD) (15) and x-ray crystallography (3BHP) (16), respectively, were not available to the investigators carrying out the quantum-chemical-based structure analysis at the time of the blind-test determination.

Data available for these quantum-chemical calculations included complete ¹³C^α chemical shift data along with 1,022 NOE-based distance constraints [supporting information (SI) Table S1 and SI Text]. Additional chemical shift data, although required for determining NOESY cross-peak assignments, were not used in the quantum-chemical structure determination and refinement process. The method for the determination of the 10 conformations (Set-bt) basically consists of three steps (see SI Text). The first of these steps is secondary-structure prediction based on ¹³C^α conformational shifts (CS), computed as described in Materials and Methods, that provides an initial abbreviated set of backbone torsional angles, namely, for those regions of secondary structure such as α-helix or β-sheet. Fig. S1 shows the corresponding distribution of ¹³C^α conformational shift (CS) values computed for the protein YnzC[1–48] using the observed ¹³C^α chemical shifts and statistical coil values from Wishart et al. (24). Two criteria were adopted for assigning the torsional constraints based on the conformational shifts. First, α-helical segments are those for which at least three consecutive residues possess CS values >1 ppm (green bars in Fig. S1), and second, for those residues that satisfied this condition, canonical backbone torsional-angle constraints were assigned, namely, φ = −60° ± 30° and ψ = −40° ± 30°. No torsional constraints were assigned for the remaining residues. Variations of the torsional angles within a tolerance range (±Λ) are not subject to energetic penalties; i.e., we use a flat-bottom pseudopotential-energy function. During the VTF procedure a tolerances range of Λ = 30° was adopted. However, a smaller tolerance range was adopted in the subsequent steps (see SI Text) allowing us to find a set of conformations that simultaneously satisfy a set of distance constraints derived from both the experimental NOEs and the ¹³C^α conformational shifts.

These torsional constraints were then used together with 1,022 NOE-derived distance constraints to generate 2,000 conformations by using the variable-target-function (VTF) procedure (25). From this ensemble, 10 conformations with lowest residual distance constraint violations (i.e., maximum distance violation <3.02 Å) were selected. For these 10 conformations, ¹³C^α chemical shifts were then computed by using the quantum-chemical DFT method, as explained in Materials and Methods. This calculation provided both backbone (φ and ψ) and all side-chain (χ) torsional constraints for all 48 residues in the sequence, as explained in SI Text.

The 1,022 NOE-derived distance constraints plus the updated backbone torsional constraints for all the residues were then used to generate a new set of conformations with the ECEPP/3 force field. Two iterations of the procedure (steps 2 and 3) follow the steps illustrated in Fig. S2. During the first iteration, 10 conformations with maximum distance violations <0.2 Å were selected among 175 generated by using both the 1,022 NOE-derived distances and the set of backbone and side-chain torsional constraints derived from the VTF procedure. In this iteration, a tolerance range ±Λ, with Λ = 20°, for the torsional constraints was adopted. The calculated ¹³C^α chemical shifts for these conformations provide a new set of backbone and side-chain torsional constraints for all 48 residues in the sequence. During this second iteration, 10 conformations with maximum distance violations <0.09 Å were selected from 330 generated by using both the 1,022 NOE-derived distances and the new set of backbone and side-chain torsional constraints derived from the previous iteration. In this second iteration a tolerance range Λ = 10°, rather than Λ = 20° used in the previous iteration, for the torsional constraints was adopted. The ¹³C^α chemical shifts for each residue of Set-bt enabled us to compute the rmsd between calculated and observed ¹³C^α values for each of these 10 conformations (see red bars in Fig. 1). A superposition of this final set of 10 conformations is shown in Fig. 2A

The significance of using ¹³C^α-derived torsional constraints for all residues.

To determine whether ¹³C^α-derived torsional constraints for all of the residues is a necessary condition to derive a good-quality set of conformations, such as Set-bt, the following test was carried out. The first iteration of our method (from here on referred to as run_A) was repeated twice. The additional two new runs differ from run_A only in the number, distribution, and type of torsional constraints used. In other words, all three runs start from the same initial conformation and with the same set of distance constraints, namely 1,022 NOE-derived distances. These different conditions were used to assess the specific contributions to the structure refinement of the specific types of torsional constraints based on quantum-chemical calculations. As a scoring function to judge the quality of the generated ensembles, for each set of conformations, we computed the maximum distance violation (using all of the 1,022 NOE-derived distances) and the backbone rmsd with respect to the average structure of the solution NMR structure (i.e., the 2JVD set) (15). Thus, the first test (from here on referred to as run_B) was carried out by removing all ¹³C^α-derived torsional constraints used in run_A, except those for the residues in the α-helix segments, namely for residues 4–18 and 24–39, which are indicated by green bars in Fig. S1. This experiment allowed us to assess the role of ¹³C^α-based torsional constraints in the turn region between the helices in determining the relative orientations of these helices. The second test (from here on referred to as run_C) was carried out by replacing all the chemical-shift constraints used in run_B by the set of backbone-torsional constraints used in the initial VTF step, namely, φ = −60° ± 30° and ψ = −40° ± 30°. This experiment allowed us to assess the value of the constraints generated by the DFT process in defining the accuracy of interhelical packing.

As discussed in the previous section, run_A yielded a set of 10 conformations with a maximum distance violation <0.2 Å and backbone rmsd of 1.7 Å relative to the average coordinates of the solution NMR structure ensemble 2JVD. The set of 10 conformations derived from run_B possessed a maximum distance violation <0.2 Å and a backbone rmsd relative to the average NMR structure coordinates of 2.6 Å. From this, we conclude that the torsional constraints on the loop region between the two helices contribute, but do not dominate, in defining the relative orientations of the two helices. On the other hand, the set of 10 conformations derived from run_C possessed a maximum distance violation <0.4 Å and a backbone rmsd relative to the NMR structure 2JVD of 1.9 Å. Accordingly, the constraints provided by the DFT calculations modestly improve the accuracy of the interhelical orientations, probably by providing somewhat more accurate side-chain packing. It should be noted that, although run_C provided a reasonably good structure, run_A could be further refined by iterative computation with DFT, to provide even better agreement with the combined chemical shift and NOE data. These results demonstrate that run_A, which makes use of the ¹³C^α-derived torsional constraints for all the residues in the sequence, provides a more accurate set of conformations than run_B or run_C in terms of the scoring functions used.

Validation of the 2JVD set of 20 conformations.

A superposition of the 20 conformations of 2JVD is shown in Fig. 2B. Computation of the ¹³C^α chemical shifts for each amino acid residue of the 20 conformations of the 2JVD set was carried out as described in Method Used to Compute the ¹³C^α Chemical Shifts in Materials and Methods. Blue bars in Fig. 1 show the mean rmsd between the observed and computed ¹³C^α chemical shifts of residues 1–48 for each of these 20 conformations. The black horizontal line (1.52 ppm) indicates the computed conformationally averaged ca-rmsd (1). Results obtained from each of the 10 conformations of Set-bt (red bars) are also shown in Fig. 1. The rmsd between calculated and observed ¹³C^α chemical shifts for each of the 20 conformations of 2JVD is higher, except for conformation no. 4 (for which these two conformations are indistinguishable), than those obtained from the Set-bt. This is not surprising because agreement between calculated and observed ¹³C^α chemical shifts was used as a constraint during the process of protein structure determination of the latter set. However, the ca-rmsd computed from 2JVD (1.52 ppm, black horizontal line in Fig. 1) is slightly better than the one computed from Set-bt (1.62 ppm, red horizontal line in Fig. 1). This is a consequence of the higher conformational dispersion in the C-terminal segment (residues 42–48) of the 2JVD structural ensemble as compared with Set-bt. This is clearly illustrated by comparing Fig. 2 A and B. The influence of this higher conformational dispersion at the C terminus on the ca-rmsd value is discussed below.

Analysis of the x-ray crystal structure of YnzC in terms of the ¹³C^α chemical shifts.

Although efforts to obtain diffraction-quality crystals of the same construct used in these NMR studies, YnzC[1–46], were not successful, diffraction-quality crystals were obtained for a slightly longer construct YnzC[1–52]. Both the YnzC[1–46] and YnzC[1–52] constructs are monomeric in solution (see Fig. S4), but YnzC[1–52] crystallized as a trimeric structure. A ribbon diagram of the 2.0-Å x-ray structure (16) (PDB ID code 3BPH, with three chains in the asymmetric unit) is shown in Fig. 2C. The three chains of this trimer, although essentially identical in conformation, exhibit some notable differences. In particular, in chains A and C, helix α2 extends through residue 52, whereas in chain B, helix α2 extends only to residue 48. The per-residue all-heavy atom rmsds (Å) between chains A–B and A–C, respectively, are shown in Fig. S3, illustrating the distribution of small structural differences among the three chains of the trimer observed in this crystal structure.

These subtle structural differences among the three chains of the 3BHP trimer result in small differences in calculated ¹³C^α chemical shift values. This is illustrated by the black, yellow and green bars, respectively, in Fig. 3. These calculations demonstrate that these ¹³C^α chemical shifts are very sensitive to the small conformational differences among the three chains of trimeric YnzC[1–52] observed in the crystal structure.

Although YnzC[1–52] is trimeric in the crystal structure, biophysical studies, including gel filtration, static light scattering, and ¹⁵N relaxation measurements, demonstrate that both YnzC[1–52] and YnzC[1–46] (used for NMR studies) are monomeric in solution (see Fig. S4); the full-length YnzC [1–77] is also monomeric in solution (15). Comparison of coordinates between the monomeric solution NMR structure and each monomer chain in the trimeric crystal structures reveals that residues 39–46 become more ordered (helical) in the crystal structure. This further extension of helix α2 in the crystalline environment is stabilized by intermolecular interactions within the trimeric structure (Fig. 2). However, for residues in the range 2–38, the NMR averaged (over all 20 models) and the x-ray averaged (over all three chains) structures are essentially identical (i.e., backbone rmsd between averaged coordinates of 0.57 Å).

Roles of ensemble averaging and protein flexibility.

The amino acid sequence of the YnzC[1–52] (3BHP) (16), YnzC[1–46] (2JVD) (15), and Set-bt structures are identical only for the first 46 residues. Hence, a straightforward comparison between these conformations requires that the rmsd and ca-rmsd for these sets of structures must be reevaluated by considering only the first 46 residues. Comparisons among these three sets are shown in Fig. 3. Some important conclusions can be derived from this analysis. First, the ca-rmsd of both 2JVD and Set-bt, black (1.54 ppm) and red (1.38 ppm) horizontal lines, respectively, in Fig. 3, are lower than the rmsd of any of the single chains of the crystal structure 3BHP (2.82, 2.06, and 2.00 ppm, for chains A, B, and C, respectively). This result is in line with a previous analysis of ubiquitin (1) showing that an ensemble of NMR-derived protein conformations provides better agreement, in terms of the ca-rmsd, than some single x-ray structures do. Fig. 3 also shows that 4 of 20 models of 2JVD (blue bars in Fig. 3) and all 10 conformations (red bars in Fig. 3) of Set-bt give lower rmsds than any of the three chains of 3BHP.

Heteronuclear ¹⁵N-¹H NOE data obtained for YnzC[1–46] demonstrate that residues 41–46, along with the C-terminal eight-residue affinity purification tag, are quite flexible on the subnanosecond time scale (see Fig. S5). In considering only the first 46 residues, computation of the ca-rmsd for Set-bt (1.38 ppm, indicated by a red horizontal line in Fig. 3) is lower than the value obtained by considering the 48 residues (1.62 ppm, indicated by the red horizontal line in Fig. 1); this effect of dynamic regions on ca-rmsd is less pronounced in calculations carried out with the 20 conformations of 2JVD, i.e., considering only the first 46 residues gives a ca-rmsd of 1.54 ppm (black horizontal line in Fig. 3), whereas considering the first 48 residues gives a ca-rmsd of 1.52 ppm (black horizontal line in Fig. 1). These results are attributable to larger differences between observed and computed ¹³C^α chemical shifts for the two C-terminal residues of Set-bt than for the 2JVD set. In fact, the absolute values of the average error, and their SD (see Eqs. 1. and 2, in SI Text, for residues 47 and 48 are 0.36 ± 3.12 ppm and 1.54 ± 3.25 ppm for the 2JVD set, and 6.27 ± 0.39 ppm and 0.19 ± 0.17 for the Set-bt set, respectively, with the standard deviation reflecting the dispersion of the conformations in each ensemble. These results illustrate how highly flexible regions, such as the C-terminal residues of these YnzC constructs, are better represented by wide, rather than tight, sets of conformations. From this point of view, the flexible C-terminal region of YnzC[1–46] is better represented by the 2JVD set than by the Set-bt.

The ca-rmsd is lower when the extent of conformational sampling is more complete. For example, for the Set-bt, analysis of the dependence of the ca-rmsd on n, with 41 ≤ n ≤ 48 being the total number of residues considered for the computation of the ca-rmsd, reveals a constant ca-rmsd value of ≈1.4 ppm for 41 ≤ n ≤ 46, a higher ca-rmsd value (≈1.6 ppm) for n = 47, and no change after the inclusion of residue 48. Among all 48 residues, residue 47 exhibits both the highest absolute value of the average error, namely 6.27 ppm, and residue 47 has the lowest heteronuclear ¹⁵N-¹H NOE values (see Fig. S5). It appears that for these two residues, the refinement protocol has resulted in an undersampling of the solution space contributing to the ca-rmsd.

In the protocol used here, the ¹³C-chemical shift-derived torsional constraints on the backbone and side-chain torsional angles used a tolerance of 10° in the last step of the procedure. Although this appears to be appropriate in well defined regions of the structure, in more poorly defined regions of the structure, these constraints are too tight and result in undersampling of the conformational space contributing to the observed chemical shift data. As we have pointed out previously, special care must be taken in using ¹³C chemical shift constraints in flexible regions of the protein structure (26). In retrospect, experimental data indicating structural flexibility, such as heteronuclear ¹⁵N-¹H NOE data (see Fig. S5), could be used to identify regions of the structure where ¹³C constraints can be used only with looser tolerances. In this blind test case, however, these data were not provided nor used in the determination of the Set-bt.

A ¹³C^α-reference-dependent analysis.

Reference inaccuracy problems (1), which affect the ca-rmsd and the rmsd values, can be detected by computing the frequency distribution of the errors between observed and computed ¹³C^α chemical shifts. Such analysis (for the Set-bt) indicates that the resulting distribution can be modeled by a Normal (or Gaussian) function with a characteristic mean value (x_o = 0.10 ppm) and standard deviation (σ = 1.08 ppm). The resulting mean value (x_o = 0.10 ppm) is very close to the ideal one (x_o = 0.0 ppm) indicating that there is no need for further reference corrections. Moreover, it also assures that 99.7% of these errors lie within 3σ (≈3.2 ppm).

A ¹³C^α-reference-independent analysis.

The Pearson coefficient R (27) provides a measure of the quality of the agreement in terms of the shielding (rather than the chemical shift) and, hence, is not affected by reference accuracy. The R values obtained for Set-bt, 2JVD, and chains A, B, and C of 3BHP are: 0.932, 0.917, 0.815, 0.881, and 0.865, respectively. These results confirm the conclusions derived from the ¹³C^α chemical shift analysis, namely through the comparison of rmsd and ca-rmsd shown in Fig. 3; the conformations of Set-bt are in best agreement with the ¹³C^α data, the conformations of set 2JVD are in good agreement, and the individual conformers of the x-ray crystal structure is in somewhat poorer agreement with these data.