|Home | About | Journals | Submit | Contact Us | Français|
We performed density functional calculations of backbone 15N chemical shielding tensors in selected helical residues of protein G. Here we describe a computationally efficient methodology to include most of the important effects in the calculation of chemical shieldings of backbone 15N. We analyzed the role of long-range intra-protein electrostatic interactions by comparing models with different complexity in vacuum and in charge field. Our results show that the dipole moment of the α-helix can cause significant deshielding of 15N; therefore, it needs to be considered when calculating 15N chemical shielding. We found that it is important to include interactions with the side chains that are close in space when the charged form for ionizable side chains is adopted in the calculation. We also illustrate how the ionization state of these side chains can affect the chemical shielding tensor elements. Chemical shielding calculations using a 8-residue fragment model in vacuum and adopting the charged form of ionizable side chains yield a generally good agreement with experimental data.
Chemical shielding is determined by the local electronic environment of a nucleus of interest. In complex molecular systems, like proteins, the chemical shift tensors contain a wealth of potentially useful structural and dynamical information. There has been significant progress in using chemical shift information for characterization of protein structure (Cornilescu et al. 1999; Lipsitz and Tjandra 2003; Shen and Bax 2007; Shen et al. 2008; Shen et al. 2009) and dynamics (Hall and Fushman 2006). First-principles quantum chemical calculations have played an important role in understanding the relative importance of various effects, such as backbone and side chain conformation, hydrogen bonding, solvation, and protein electrostatics, on the chemical shielding of nuclei in proteins. These calculations have also provided parameters and conceptual ideas for the development of phenomenological models (Xu and Case 2001; Xu and Case 2002). Despite their contributions to our understanding of the chemical shielding phenomena and to model-building efforts, first-principles quantum chemical calculations have not become a routine practice in predicting chemical shifts in proteins because of several difficulties. First, chemical shifts could be influenced by a multitude of factors. In general, the backbone torsion angles, the side chain orientation, the hydrogen bonding and direct interactions with nearby side chains, the long-range intra-protein electrostatics, and solvent effect can all contribute. It is challenging to include all foreseeable effects and still be computationally efficient. Second, any experimental protein structure we use for calculating chemical shifts might require some geometry optimization in order to get an accurate prediction. This is particularly true for α-carbons because their calculated chemical shielding is very sensitive to bond lengths and bond angles (de Dios et al. 1993a). Third, proteins rapidly sample the available conformational space, hence no single snapshot (X-ray) or even a finite bundle of NMR structures of a given protein can fully represent the conformational ensemble (and dynamic ensemble averaging) for which the experimental chemical shifts are measured. Moreover, when using a solution structure for calculations, one has to consider not only the structural ensemble but also at what conditions (pH, temperature, etc) the structure was obtained. There can be fluctuations which may not only affect some structural parameters but also alter the ionization form of the side chains of Asp, Glu, Arg, and Lys (Vila and Scheraga 2007). Fourth, there are systematic errors associated with the level of theory used and the basis set chosen for calculating the electronic structure. In practice, an offset may exist that will need to be corrected if an absolute isotropic chemical shift prediction is the aim (Cai et al. 2008). An additional complication here is that this correction should differ between various residue types (Xu and Case 2001).
Despite these difficulties, theoretical chemical shift tensors are valuable in peptide and protein structure validation and refinement (Wylie et al. 2009) and therefore worth exploring. In this paper, we chose to focus on amide 15N chemical shielding tensor for several reasons: (i) it is influenced by numerous factors and, therefore, is arguably the most difficult case for computational prediction (de Dios et al. 1993b); (ii) perturbations in amide chemical shifts are often used for mapping protein-ligand interactions (Zuiderweg 2002); and (iii) amide 15N chemical shift tensor contributes to 15N spin relaxation rates that are widely used for analysis of protein dynamics (Fushman and Cowburn 2001; Hall and Fushman 2006). Thus, first-principles calculations may shed light onto the relationship between structural changes in proteins upon ligand binding and the accompanying changes in their chemical shifts. The quantum chemical calculations can also provide an understanding of the structural basis for site-to-site variability in 15N chemical shift tensors (Fushman et al. 1998; Hall and Fushman 2006).
Previous theoretical works (de Dios et al. 1993b; Le and Oldfield 1996; Xu and Case 2002; Cai et al. 2008) have established that 15N chemical shielding depends on the following factors with none seemingly dominating: ϕ, ψ, χ1, preceding side chain’s identity and conformation, hydrogen bonding partners, electrostatic interactions, and solvent effect. Whereas the short-range interaction with hydrogen-bonding partners can be included explicitly and exactly, model treatment has been required for long-range electrostatic interactions. Point charge representation can be an option here. For example, in addition to the main fragment, charges can be incorporated for the remaining atoms in the protein using some charge set. Thus, it was found that if the charge field perturbation (CFP) effects are included in the calculation of carbon chemical shielding in amino acids, the correlation between the theory and experiment is slightly improved (de Dios et al. 1994). The improvement is mainly for the sp2 carbonyl carbon, which is more sensitive to electrostatic field effects. The sp2-hybridized amide nitrogen shares this kind of sensitivity (Bader 2009) and showed some prediction improvement with inclusion of charge field perturbation in a dipeptide model (de Dios et al. 1993b). However, a later application to helical residues (Le and Oldfield 1996) indicated that a static charge field is inadequate to account for the long-range electrostatic field contribution to 15N shielding in an α-helix. There are probably two reasons for this: (i) each type of residue has a fixed set of approximated charges, which may not be accurate enough for the purpose of calculating 15N chemical shifts because the multiple-bond character of the peptide group makes the peptide nitrogen highly polarizable and hence very sensitive to the electric field that the charges generate; and (ii) the ionization state of Asp, Glu, Arg, and Lys is not well determined and can vary depending on the protein ensemble, temperature, and pH; this renders representation of the side chains of these residues as either neutral or charged somewhat arbitrary.
To address these issues we performed density functional calculations of backbone 15N chemical shielding for helical residues in protein G (GB3) and compared the results with experimental data from solid-state and solution NMR measurements. To assess the contributions from intra-protein electrostatic interactions, the density functional calculations were combined with charge field perturbations.
All calculations reported in this paper were performed using the GAUSSIAN03 suite of programs (Frisch et al. 2004). We used density-functional theory (DFT) with three-parameter Becke-Lee-Yang-Parr (B3LYP) exchange-correlation functional (Becke 1988; Lee et al. 1988; Becke 1993). The atom coordinates were taken from the best representative conformer (PDB code: 2OED) of the ensemble of GB3 solution NMR structures (Ulmer et al. 2003). The bond lengths and angles were taken directly from the experimental structure without geometry optimization.
Only ‘classical’ helical residues, A26 through Y33, of the GB3’s α-helix (that spans residues D22-N37) are examined in this paper. According to the definition used here, a classical helical residue i is hydrogen bonded through its NH group to residue i−4 and through the CO group to residue i+4. Thus the backbone nitrogen of residue i has a direct hydrogen bonding partner i−4 and an indirect hydrogen bonding partner i+3 (Fig. 1). To avoid terminal artifacts, four residues at either end of the α-helix were not examined here because they are missing either a direct or an indirect hydrogen-bonding partner.
Several models of different levels of complexity are investigated here (illustrated in Fig. 1b for the case of amide nitrogen of K28):
Model A: a simple dipeptide model containing only residues i and i−1. If not indicated directly, the ionizable side chains are assumed to be in their charged state. In some cases, as indicated in the text, the ionizable side chains were altered to be neutral.
Model B: a main fragment containing residues i and i−1 (as in Model A) and two additional fragments of the hydrogen bonding partners, both direct and indirect. This results in a “three-fragment” model. In some cases, as specified in the text, the side chains of the hydrogen bonding partners were altered to be neutral or modified to be −CH3.
Model C: a “long-chain” model containing a stretch of residues from the direct hydrogen bonding partner to the indirect hydrogen bonding partner as one main fragment. That is, for classical helical residues the “long chain” includes residues from i−4 to i+3. In some cases, as specified in the text, some side chains of the “long chain” were altered to be neutral or changed to −CH3, or the “long chain” was extended to include an additional residue at one end.
In our calculations, the main fragments in various models had the N-terminus capped by a formyl group (−COH) and the C-terminus capped by an amino group (−NH2), and the hydrogen bonding partners in Model B were modified to be CH3-CO-(NH-CH(R)-CO)-NH-CH3. When the charges from the atoms of the rest of the protein are included in the calculation, we refer to it as Charge Field Perturbation (CFP) calculation. If the charges are not included, we call it “vacuum” calculation. In CFP calculations, the point charges were all taken from AMBER charge set (Cornell et al. 1995) with the overall non-zero charge for ionizable residues. The calculations were performed using local dense basis sets (Chesnut et al. 1993). For the residue of interest (i), we applied a 6–311+G(2d,p) basis set for the Ni, Hi,Cαi, Ci−1, and Oi-1 atoms (shown in bold in Fig.1a). Where applicable, the 6–311+G(2d,p) basis set was also applied to similar atoms of its direct and indirect hydrogen bonding partners (shown in blue in Fig.1a), as well as the “second pair” (Xu and Case 2002) of hydrogen bonding partners (shown in red in Fig.1a). A 4-21G basis set was applied to the remaining atoms in the model. The charge field perturbation gauge-including atomic orbital method (Ditchfield 1974; Wolinski et al. 1990; de Dios and Oldfield 1993) was used.
It was suggested (Xu and Case 2002) that a 7–9 residue sequence should be used for chemical shift calculations in the α-helix conformation in order to eliminate the terminal artifacts. The terminal effects may be due to (i) the lack of hydrogen bonding partners for the residues at either end; or (ii) the underestimation of the interaction with the dipole moment of the helix.
In order to analyze the magnitude of the hydrogen bonding and helical dipole effects without side chains complicating the picture, all non-alanine residues other than residues i and i−1 were altered to be alanines (see Table 1), and we compared isotropic chemical shieldings from Model A to Model C vacuum calculations. As shown in Fig. 2a, the first pair of hydrogen bonding partners deshields 15N. The deshielding ranges from 4.75 ppm to 9.89 ppm depending on the strength of hydrogen bonding. The effect is the strongest for K31 because the peptide plane containing the amide group of K31 makes the strongest hydrogen bonds. In fact, the distance between K31’s amide proton and E27’s carbonyl oxygen is 1.83 Å, compared to the 1.89 to 2.86 Å range for the rest of the hydrogen bonds in the GB3 helix. In addition, the O-H distance for the corresponding indirect hydrogen bond is 1.89 Å. Model C deshields 15N further from Model B, by the amount ranging from 4.10 ppm to 6.59 ppm (Fig. 2a). This is due to the inclusion of atoms from residues i−3, i−2, and i+2 in Model C, which introduced new aligned NH and CO groups, namely NH groups from residues i−2 and i+2, and CO groups from residues i−3 and i+1 (see Fig. 1). These groups contribute to the helical dipole moment and bring about electrostatic field contribution to 15N chemical shielding, in agreement with previous findings (Le and Oldfield 1996; Xu and Case 2002). Our calculation showed an increase by 18.5 Debye on average in the dipole moment from “Model A neutral” to “Model B neutral −CH3”; and on average an additional increase by 10.2 Debye from “Model B neutral −CH3” to “Model C neutral −CH3” (Table S1). Since the “long-chain” model further deshields the 15N compared to the “three-fragment” model and this deshielding differs for different residues (Table 1), we conclude that the “long-chain” model includes some new nontrivial contributions, which are not present in the simpler models.
We note here that in the calculations with charged ionizable side chains (Fig. 2b), the contributions from the first pair of hydrogen bonds and from the additional helical dipole are similar to those in the calculations with neutral ionizable side chains (Fig. 2a). This indicates that the side chain charges do not interfere with hydrogen bonding effects and helical dipole effects.
Previously, when building the database for SHIFTS (Xu and Case 2002), all ionizable amino acids took neutral forms because it appeared that in gas-phase calculations, neutral side chain provided a better model for solution chemical shifts than the charged side chains did. More recently, it was suggested (Vila and Scheraga 2007) that the whole protein and its chemical environment should be taken into account when considering the proton binding/release equilibria and then the chemical shifts should be calculated as a weighted average. That study also found that for Cα the results obtained for Asp and Glu are almost systematically worse when using charged rather than neutral side chains, presumably because these two amino acids have much shorter side chains compared to Arg and Lys. When charged, these side chains may have a bigger influence on the Cα chemical shielding. Here we explore how the 15N chemical shift can be affected by the ionization form of the side chains.
As shown in Table 1, the side chain charge on the residue of interest has an apparent effect on its 15N chemical shielding (compare the results from “Model A” and “Model A neutral” calculations for E27 and K31). Deprotonation of the side chains introduces deshielding while protonation introduces shielding, similarly to what was observed for Cα (Vila and Scheraga 2007). For example, the negative charge on the side chain of E27 deshields its 15N by 11.74 ppm while the positive charge on the side chain of K31 shields its 15N by 6.34 ppm in Model A.
Unlike Cα, surrounding side chain charges on other residues can influence amide nitrogen chemical shielding. Take K31 as an example. Its “Model B −CH3” calculation (130.10 ppm) yields a more shielded value than its “Model B neutral −CH3” calculation (124.33 ppm) for 15N. This shielding effect of 5.77 ppm is due to the charge on K31. In the Model B calculation where K31’s hydrogen bonding partner E27 is also charged, the amide nitrogen of K31 has a chemical shielding of 125.80 ppm. This 4.3 ppm deshielding from the “Model B −CH3” calculation is mainly due to the charge on E27. This example illustrates that the charges on other side chains can make 15N chemical shielding calculation more interesting but also more challenging because it introduces another source of variation. In the case of neutral side chains, on the contrary, their identity does not play a significant role, such that they can be safely replaced by a −CH3 group. This is supported by the close resemblance of the results from “Model B neutral” and “Model B neutral −CH3” calculations and, likewise, of the results from “Model C neutral” and “Model C neutral −CH3” calculations.
We used two approaches when calculating 15N chemical shielding with Model C. The first approach makes all ionizable side chains neutral in the Model C vacuum calculation. The comparison between experimental solid-state chemical shifts and the calculated chemical shieldings is shown in Fig. 3a. The correlation is poor (|r|=0.52), and the regression line has an intercept of 181.00 ± 41.00 ppm and a slope of −0.49 ± 0.33. The slope is far from the ideal value of −1. No obvious outlier can be identified. The second approach represents all ionizable side chains in their charged states. This yields an obvious outlier, E27 (see Fig.4a). By excluding E27, a better correlation is achieved (|r|=0.93) and the linear regression yields a slope of −0.84 ± 0.15 and an intercept of 224.91 ± 18.83 ppm (Fig. 3b). A comparison with solution NMR data yields similar results: |r| = 0.91, slope = −0.98 ± 0.20, and intercept = 242.87 ± 24.70 ppm (Fig. 3c). In this approach it is important to include all charges that are close to the 15N site of interest (Bader 2009). Consider E27 as an example. In the Model C vacuum calculation that includes residues from A23 to F30, we obtained an isotropic chemical shielding of 113.98 ppm for the amide nitrogen. This calculation did not include K31, which makes a salt bridge with the side chain of E27. However, extending the long chain to include K31 in the Model C vacuum calculation yielded a more shielded value of 117.94 ppm, bringing it closer to the experimental value (Fig. 4a). This shielding effect of 3.96 ppm is due to the charge on K31’s side chain. In comparison, a proton unit charge placed on the ε-amine of K31 produced a 4.77 ppm shielding effect in the Model A vacuum calculation (Table S3). Interestingly, having E27 in its neutral form can significantly deshield its amide nitrogen, as suggested by the difference (11.74 ppm) between its Model A and “Model A neutral” calculations (Table 1). Indeed, making E27 neutral while keeping other ionizable side chains charged (in Model C) yields an isotropic chemical shielding of 124.25 ppm. Thus, protonation of the side chain can make E27 not an outlier.
In principle, there can be other possible combinations of ionization states of D22, D24, D27, K28, K31, and D36 that are involved here (Table 1). The two approaches shown here emphasize the complexity of calculating backbone 15N chemical shielding when dealing with uncertainties of the charge states of ionizable side chains. On the other hand, a detailed comparison with experimental NMR data might provide useful information on the ionization state of these residues.
Counterions were not explicitly included in our calculations, primarily because the amino acid composition and the structure of GB3’s α-helix are such that every charged side chain (except for D22 at the very N-terminus and D36 at the very C-terminus) has a salt-bridge partner (E27–K31, E24–K28), which could naturally balance its charge effect. Of the long-chain fragments used in Model C calculations, the one for K28 was neutral (as was the extended fragment for E27, see above), while the rest of the fragments had one or two unbalanced charges. Including the extra residues necessary for charge balancing would have further increased the size of the fragment. However, the effect of those unbalanced charges is generally small, as can be inferred from Table 1 (compare, for example, Model C with Model C −CH3), perhaps because of the rather long distances from most of these charged groups to the amide nitrogen of interest. To explore this effect further, we performed a dipeptide-model calculation in which we included point charges mimicking the side chains that could serve as “counterions” but were not present in the corresponding Model C fragments. The results (Table S3) show that the effect of these additional charges on 15N chemical shielding is generally small, except for E27, where including the side chain charge of K31 introduces a significant shift in the shielding, as discussed above.
Comparing the shielding tensor’s principal values between the calculations with “Model C neutral” and “Model C” (Table S2), the variation mainly exits in the least and most shielded components, σ11 and σ33, while the intermediate component, σ22, appears less susceptible to model selection. Residue E27 showed the biggest difference in all three principal components between these two calculations, with a 20 ppm difference in σ11. There is a generally good agreement between these two models in the orientation of the principal components of the chemical shielding tensor (Table S2).
Experimental data for all three principal components of GB3’s 15N chemical shift tensors are not available. However, these data exist for most residues of protein GB1 (Wylie et al. 2007), which is highly homologous to GB3. In fact, all helical residues of GB3 are present in GB1, except for E24 and A29, which are an alanine and a valine, respectively, in GB1. Therefore we compared our calculated principal values of the 15N chemical shielding tensor for A26, E27, K28, K31, and Y33 in GB3 with the corresponding values of the chemical shift tensor of these residues in GB1, obtained by solid-state NMR measurements (Fig. 5). We found that the individual principal values of the calculated chemical shielding tensor are in good agreement with the experimental data. Over the whole range of variation in the tensor’s principal values, the difference between the results of “Model C neutral” and “Model C” calculations is small (Fig. 5).
The difference between a vacuum calculation and a CFP calculation is a useful indicator of whether all necessary electrostatic effects are taken into account in the vacuum calculation. For example, when using a “three-fragment” model (Model B) there are large differences between the CFP calculation and the vacuum calculation (Table 1). This suggests that the fragment that is used in the vacuum calculation is too small. When using Model C which includes every residue from the direct hydrogen bonding partner to the indirect hydrogen bonding partner, not only the difference between vacuum and CFP calculations becomes smaller in general (Table 1 and Fig. 4b) but also the correlation between the vacuum calculation and the experiment becomes better, with a |r| value of 0.93, compared to 0.68 for Model B vacuum calculation and 0.54 for “Model B CFP” calculation (E27 was not included in this analysis).
From Table 1 and Fig. 4b, we observe that a simple CFP calculation for one model is close to a vacuum calculation for the next-complexity-level model. That is, Model A CFP calculation is very similar to Model B vacuum calculation and, likewise, Model B CFP is in general very close to Model C vacuum calculation. Although the CFP calculations include all point charges representing the atoms from the rest of the protein, only certain charges that are close enough to the fragment under consideration can influence the distribution of the electrons by perturbing the wave function and hence potentially affect the 15N chemical shielding. For example, Model B CFP calculation shows on average ~7 ppm deshielding of 15N compared to the Model B calculation, but including only point charges for the remaining helical residues from Model C yields most of the deshielding effect already (Table 1).
We performed density functional calculations of backbone 15N chemical shielding tensor for selected helical residues in protein G (GB3) and compared the isotropic chemical shielding and the principal values of the shielding tensor with experimental chemical shift data. We explored the effect of electrostatic interactions in a protein on the calculated 15N chemical shielding and found that:
Note that this study focused on an α-helix, where the spatial relationship between various atoms is different from that in, for example, extended conformation. Therefore these conclusions apply specifically to helical residues and may or may not be applicable to other secondary structures. Similar calculations (currently underway) for the other parts of the protein are expected to address this issue. Although this study covered different residues in the helix, it is inevitably limited by the amino acid composition of GB3. Comparison with experimental data for a broad range of proteins is required in order to reach a better understanding of the accuracy of the computational approaches used here and the ways to improve them.
Supported by NIH grant GM065334 to DF and by American Chemical Society Petroleum Research Fund (44481-G6) and Francqui Foundation to DSK. We thank Gabriel Cornilescu, Chad Rienstra, and Ben Wylie for making available detailed experimental chemical shift data for comparison.