Computational Protocol to Evaluate Side-Chain Vicinal Spin–Spin

Article pubs.acs.org/JCTC

Cite This: J. Chem. Theory Comput. 2019, 15, 4252−4263

Computational Protocol to Evaluate Side-Chain Vicinal Spin−Spin Coupling Constants and Karplus Equation in Amino Acids: Alanine Dipeptide Model J. San Fabián, S. Omar, and J. M. García de la Vega* Departamento de Química Física Aplicada, Facultad de Ciencias, Universidad Autónoma de Madrid, 28049 Madrid, Spain

Downloaded via UNIV OF SOUTHERN INDIANA on July 23, 2019 at 18:23:32 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: A computational protocol has been applied to the alanine dipeptide model in order to study the side-chain conformation, the calculated spin−spin coupling constants involved in the side-chain χ1 angle, and theoretical extended Karplus equations developed for amino acids. Two structures within the backbone secondary conformation are used to predict coupling constants which in addition are employed to analyze the effect of those structures on the resulting Karplus extended equations. The number of Fourier coefficients to be included within the Karplus equations is critically analyzed. Wave function and density functional methods as well as different basis sets are compared to find appropriate Karplus coefficients in the most efficient way. The influence of exchange and correlation functionals and the solvent effect on the calculated couplings are considered.

1. INTRODUCTION The structure of proteins, which are the central building blocks of life, plays a key role in chemical and biological processes. Knowledge of three-dimensional structures is essential to fully understand the protein interaction with other molecules and therefore its functional mechanism and biological function.1−3 Nuclear magnetic resonance (NMR) is a popular technique for studying biological systems in near physiological conditions4,5 and is crucial for both the assessment of ligand bindings and conformational studies.6 This popularity mainly results from the development of a multitude of experiments designed to obtain line-widths, scalar couplings, residual dipolar couplings, and Nuclear Overhauser Effects (NOEs). Among them, scalar couplings, together with NOEs, are the two most important NMR tools and the most used for the configurational and conformational analysis of small organic molecules.7−10 NMR spin−spin coupling constants (SSCCs) of proteins are valuable experimental parameters to investigate conformational properties.11−16 In particular, vicinal coupling constants 3J, those between three bonds, show a dependence on the dihedral angle θ between the coupled nuclei X and Y. This dependence between 3J and θ was established almost 60 years ago,17,18 3

J XY = C0 + C1 cos(θ ) + C2 cos(2θ )

conformation in proteins using relationships between vicinal SSCCs and torsion angles χ1 has been developed, although to a lesser extent.24−26 The theoretical determination of Karplus equations, together with the empirical results, can help the conformational analysis of biomolecules from different points of view. Nowadays, theoretical results using wave function (WF) or density functional theory (DFT) methods predict correlations between spectroscopic parameters and molecular structure that approximate those derived from empirical evaluations. In addition, theoretical calculations can be applied to specific biochemical models with an environment similar to that of interest. For instance, a full set of Karplus coefficients can be obtained from the SSCCs involved in the side-chain bond for each amino acid (AA) that will include explicitly the substituent effects27 and therefore minimize the errors. Besides, the conformational space, usually limited in experimental systems, can be fully covered in theoretical approximations. The Karplus equations which, due to empirical limitations, are usually restricted to three coefficients, eq 1, could be extended from theoretical results. Moreover, errors derived from the limitations of experimental Karplus coefficients can be evaluated from theoretical results, and estimations in the uncertainty can be established. In molecular fragments X−M−N−Y, vicinal SSCC between two nuclei X and Y, 3JXY, can be related to the torsional angle θ formed by the planes X−M−N and M−N−Y. A general

(1)

Several studies on the relationship between vicinal SSCCs and the dihedral backbone angles in polypeptides have been carried out by the groups of Bax,19,20 Case,21,22 and Rü terjans.23 Also, the determination of the side-chain © 2019 American Chemical Society

Received: February 11, 2019 Published: May 29, 2019 4252

DOI: 10.1021/acs.jctc.9b00131 J. Chem. Theory Comput. 2019, 15, 4252−4263

Article

Journal of Chemical Theory and Computation

2.1. AA Molecular Geometry. Most of the AA usually have either the α (α-helix with ϕ ≈ −64° and ψ ≈ −44°) or β conformation (β-sheet with ϕ ≈ −121° and ψ ≈ 128°). These two main conformations are the result of combining AAs into peptides and proteins that present, in addition to local steric interactions, through space interaction among the AA residues. When a small dipeptide model, i.e., one containing two peptide linkages (−CO−NH−), is used, the noncovalent interactions from further residues and from the surrounding media do not exist. Due to the lack of these noncovalent interactions, a full geometry optimization of this dipeptide fragment leads to conformations different from those indicated above (α or β), and therefore, the resulting geometry is less attractive for protein property studies. Thus, in this work the geometry optimizations will be restricted to those secondary structures (α-helices and β-sheets). Effects on the side-chain SSCCs for these two conformations are compared and analyzed. Besides, AA model geometries driving χ1 angle must be obtained, that is, the additional parameter χ1 = N′−Cα−Cβ− Hβ1 will be constrained between 0° and 360° at intervals of 30°. Due to the symmetry of the alanine (Ala) side-chain CβH3 group, only two (0 and 60°) or four conformations (0, 30, 60, and 90°) are necessary to obtain 6 or 12 values for each SSCC, respectively. As indicated below, the choice of a specific hydrogen Hβ1 introduces some nonequivalences or asymmetries in the fragment within the side-chain group of protons. 2.2. SSCC Calculations. Once the geometries are obtained, the vicinal SSCCs around the Cα−Cβ bond are calculated. For the particular case of Ala, see Figure 1, three vicinal SSCCs can be calculated around the Cα−Cβ bond: 3 α β 3 JH ,H , JN′,Hβ, and 3JC′,Hβ.

formulation between vicinal couplings and torsional angles is based on Fourier series, P 3

J XY = C0 +

Q

∑ Cn cos(nθ) + ∑ Sn sin(nθ) n=1

n=1

(2)

28−30

From theoretical calculations, it was shown that eq 2 with terms up to P = 3 and Q = 2 includes the most important contributions. The well-known classical Karplus (eq 1) uses the terms up to P = 2 and Q = 0, neglecting the sine terms that are responsible for the asymmetry of vicinal couplings around θ = 180°. Moreover, although the coefficients C3 and S1 are usually small, the coefficient S2 is not negligible and its accurate determination could improve to resolve ambiguities on the molecular conformation.31 The Fourier coefficients, Cn and Sn, have been derived in many ways. Generally, the experimental SSCCs are fitted to torsional angles obtained from high-resolution X-ray crystallography and for specific substituent patterns.19,20,32−34 Alternatively, Schmidt et al. obtained both Fourier coefficients and torsional angles from a redundant set of experimental SSCCs using a self-consistent method.23,35 A less explored approach consists of the employment of theoretical SSCCs and torsional angles to obtain Karplus equations for biomolecules.21 Once these Fourier series, also called extended Karplus equations, have been obtained and parametrized, they can be used to predict, for instance, the side-chain conformation. Therefore, initially we need to obtain two important sets of parameters: vicinal SSCCs and their corresponding dihedral angles. In addition to the torsional or dihedral angle, vicinal SSCCs 3 JXY depend on other factors as the local geometry (bond lengths and angles) and the substituents attached to the X− M−N−Y fragment.18 In a first approximation, we consider that those effects are implicitly included, at least partially, into the Fourier coefficients when they are obtained directly for specific AA model fragments.

2. COMPUTATIONAL PROTOCOL The computational protocol for evaluating SSCCs and Karplus equations in AAs (see Scheme 1) is summarized as follows. Scheme 1. Computational Protocol Figure 1. Molecular model of Ala dipeptide. Definition of ϕ (C−N′− Cα−C′), ψ (N′−Cα−C′−N), and χ1 (N′−Cα−Cβ−Hβ1) are also indicated. Newman projection of the χ1 angle.

Initially, the proper or more accurate way to calculate these couplings should be ab initio WF methods that have proved to give accurate results in small molecules. In this work, these ab initio calculations will be carried out at the limit of our computational resources. However, due to the size of the AA fragments and the large amount of calculated SSCC required, it is more appropriate to perform computations using DFT methods. Consider the future application of the results to larger AA residues. Besides, the accuracy of a specific functional for the calculation of SSCCs depends largely on the basis sets, the exchange-correlation term, and the type of SSCC.36−39 Therefore, DFT calculations must be tested in specific cases to find the best basis set and functional for each type of SSCCs. The lack of experimental SSCCs to cover the 4253


Article


and DFT results, basis sets effects, choosing a functional showing the best quality/cost, etc. Two main topics are studied in this work: (1) the ability of theoretical calculations (WF and DFT) to predict SSCCs and to obtain Karplus equations good enough to be used in conformational protein studies. (2) The use of Karplus equations to estimate the problems, inaccuracies, and the intrinsic deviations of the Karplus equation because of its limited number of terms, the secondary effects as the substituent or local geometry effects, the backbone conformation, etc. With this study, we expect to shed light on the possible problems of Karplus equations and on the competition between WF and DFT methods, basis sets, and Fourier series truncation for determining the geometries of peptides. The rest of the paper is organized as follows: in the next section 3, the main features of the computational details are described. In section 4, results for the Ala model using the proposed computational protocol for WF and DFT methods are described, discussed, and compared with the available experimental data. Finally, some general conclusions are drawn.

full torsional space obliges us to compare DFT results with those obtained using WF calculations. 2.3. Fourier Coefficients. Both SSCCs and dihedral angles are used to obtain the Fourier coefficients for the three types of couplings. In Ala, because of the free rotation of the methyl group, the experimental couplings can only be related with C0 or C0 + C3 Fourier coefficients (see below). From six values of 3JXY and their respective dihedral angles, it is possible to obtain up to six Fourier coefficients C60, C61, C62, C63, S61, and S62 (set66), where the superscripts show the total number of SSCCs available or used to obtain the Fourier coefficients. This set will be derived from either the costly ab initio WF calculations where only six values are calculated or from DFT calculations when comparisons are required. On the other hand, if 12 values of SSCCs are available, it is possible to 12 12 12 12 12 obtain up to 12 coefficients C12 0 , C1 , C2 , ..., C6 , S1 , S2 , ..., and 12 S12 (set12 ). Obviously, in this case, we can also use six of the 5 SSCCs to obtain the set66. One can also take only the first six coefficients in the set1212, this set can be named set612. It should be noted that set66 and set612 are slightly different and they are related by C06 = C012 + C612 , C16 = C112 + C512 ,

S16 = S112 − S512 ,

C26 = C212 + C412 ,

S26 = S212 − S412 ,

C36 = C312

3. COMPUTATIONAL DETAILS In order to study the side-chain conformation and SSCCs in Ala, we need to choose a molecular model large enough for incorporating all the properties under study and small enough to carry out the computational calculations with a rational cost and within the available resources. The dipeptide model chosen (N-formyl L-alanine), shown in Figure 1, is appropriate considering that the SSCCs are quite local properties and that the α- and β-substituent effects are included. Moreover, the many intrinsic local geometry effects are also considered since the SSCCs were calculated using partially optimized geometries. This model comprises eight hydrogen and eight first row atoms. It represents a simplification of the so-called dipeptide model40 after the methyl terminal groups are replaced by hydrogen atoms. These models have been frequently used to study and analyze the properties of AA side-chains.40 The dipeptide model including two methyl terminal groups40 has also been used in this work to compare it with our main model (see below). The full and partial optimized geometries have been carried out at the B3LYP/6-31G(d,p) level41−45 using the Gaussian program.46 This level of theory has been used satisfactorily to obtain geometries in previous AA studies.31,47,48 SSCCs have been calculated using WF and DFT calculations. Three WF methods have been used: the properly called SOPPA (secondorder polarization propagator approximation) that uses single and double excitation coefficients from the Møller−Plesset (MP) methodology,49−51 SOPPA(CC2) which is a SOPPA method that uses the single and double amplitudes from CC2,52 and SOPPA(CCSD) that uses single and double amplitudes from CCSD (coupled cluster singles and doubles).53 The last method presents the best description of the electron correlation.52 Although it is the most expensive, it will be taken as the best reference in this work. Dalton program54 has been used to perform the calculations with the three WF methods. Within DFT, the B3LYP functional41,42 will be used to check a series of initial effects and geometrical traits. Finally, several additional functionals will be tested. Those functionals have been used previously for SSCCs.38,39,55,56 In addition, some

(3)

To analyze the vast amount of results, sets of 6 and 12 SSCCs or Fourier coefficients, the following root mean squared deviation (rmsd) is used as a statistical parameter to compare pairs of results. For instance, a set 1 and a set 2 are compared using 2π

rmsd =

∫0 (J set1(ϕ) − J set 2 (ϕ))2 dϕ 2π

(4)

Within this definition, rmsd between two Fourier series represented by eq 2 is 1 ijjj 2 j∑ (ΔCn) + 2 jj n = 1 k P

rmsd =

(ΔC0)2 +

yz

∑ (ΔSn)2 zzzzz Q

n=1

{

(5)

with ΔCn = − and ΔSn = − To compare different methods and/or basis sets, it is more convenient to assemble the relative rmsd for the three studied SSCCs. Thus, we define the average weighted rmsd (awrmsd) values in %, where the relative factors or weights are the average couplings taken as the respective |C0| values. Cset1 n

Cset2 n

Sset1 n

Sset2 n .

ÅÄÅ ÑÉÑ Å ÑÑ ij rmsd yz ij rmsd yz ÑÑ 1 ÅÅÅÅijj rmsd yzz j z j z Ñ zz jj zz jj zz awrmsd = ÅÅjjj + + j |C | z 3 j |C | z3 ÑÑÑÑ × 100 3 ÅÅÅk |C0| z{3 0 0 k { k { J HH J CH J NH Ñ ÅÅÇ ÑÑÖ

(6)

The values taken for C0 are those calculated at the SOPPA(CCSD)/aug-cc-pVTZ-J, that is, 6.70, 3.92, and −3.14 Hz for 3JHH, 3JCH, and 3JNH, respectively. In this work on protein side-chain SSCCs, we explore different problems found when trying to settle appropriate Karplus equations: the validity of dipeptide model, geometry constrains, the chosen reference when artificially driving χ1 angle, the truncation of Fourier series, differences between WF 4254


Article

Journal of Chemical Theory and Computation Table 1. Basis Sets Functions (Uncontracted/Contracted) For the Individual Atoms and for the Ala basis set 6-31+G*-J 6-31++G**-J 6-311+G*-J 6-311++G**-J aug-cc-pVTZ-J

C, N, O 37/21 37/21 41/26 41/26 55/46

[14s6p1d/7s3p1d] [14s6p1d/7s3p1d] [15s7p1d/9s4p1d] [15s7p1d/9s4p1d] [15s6p3d1f/9s5p3d1f]

H

Ala model

7/4 [7s/4s] 11/8 [8s1p/5s1p] 8/6 [8s/6s] 12/10 [9s1p/7s1p] 24/20 [10s3p1d/6s3p1d]

352/200 384/232 392/256 424/288 632/528

solvent effect PCM to the B3LYP/6-31G(d,p) level; and (iii) calculated SSCCs include the model PCM solvent to the geometries optimized in the previous point (ii). Two characteristic solvents were considered: water and dimethyl sulfoxide (DMSO), with dielectric constants of 78.36 and 46.83, respectively. For Ala residues, the experimental side-chain SSCCs are averaged values due to the free rotation of the methyl group. At room temperature ⟨3JXY⟩ ≈ C0 and only the first Fourier coefficient can be used to test the calculated values.67 From 17 residues of Ala, Pérez et al.35 obtained empirical values between 5.58 and 6.30 Hz (average 5.88 ± 0.20 Hz) for 3JHH, between 3.66 and 4.27 Hz (average 3.96 ± 0.15 Hz) for 3JCH and between −2.17 and −2.90 Hz (average −2.61 ± 0.24 Hz) for 3JNH (only 16 residues). All the results presented in this work refer to 1H, 13C, and 15 N isotopes. Experimentally, the isotope most used for nitrogen is 15N. Therefore, the 3JNH SSCCs and the Fourier coefficients for the isotope 14N can be obtained dividing by −1.40242, the quotient between the gyromagnetic ratios. As the sign is not usually determined experimentally, some experimental results are provided for 15N as positive and the SSCCs or the Fourier coefficients derived from them show the sign changed with respect to the theoretical results.68

functionals with a large amount of Hartree−Fock (HF) exchange38,39 have also been tested as well as the effect of a fraction HF exchange on the SSCCs.42,57,58 Basis sets used in these calculation are the 6-31+G*-J, 6-31+ +G**-J, 6-311+G*-J, 6-311++G**L-J,59 and the larger aug-ccpVTZ-J.60 These basis sets were especially developed for calculations of the coupling constants.59,60 WF methods have advantages and drawbacks over present DFT methods. The main benefit of WFs is that the exact solution can be approached systematically. That is, in general, the results of a property improve when electron correlation is considered and a larger basis set is used. In SSCC calculation this last statement is not always true because a small basis set that describes properly the density close to the coupled nuclei using the well-known tight functions can give better values than a large basis set that does not take this point into account. Therefore, it is recommended to use a basis set specially optimized for SSCCs as those indicated above. Other interesting basis sets are those developed by Jensen;61−63 however, and although those basis sets allow one to study the convergence with basis set size, they grow too quickly and its use is limited for medium-sized molecules. Taking into account these points, we assume that the WF results will improve hierarchically from SOPPA < SOPPA(CC2) < SOPPA(CCSD) and as the basis set size increases: 6-31+G*-J < 631++G**-J < 6-311+G*-J < 6-31++G**-J < aug-cc-pVTZ-J. The main disadvantage of WF methods is the high computational cost that limits both the method and the basis set length as the size of molecule increases. Table 1 shows the size of the basis set functions for Ala. The largest one has 528 contracted basis functions making the WF calculations quite expensive. This problem can be solved using DFT methods that are much more economical than those of WF theories. However, their main drawback is that they cannot be systematically improved. DFT is an exact theory, but it must be approximated with nonsystematic routes since the exact functional is unknown. This fact introduces some conceptual problems. In addition, the SSCCs calculated by DFT are very dependent on the type of coupling, the used functional, and basis set in a way that is not always consistent. That is, a “worse” functional for instance in a lower step in the Jacob’s ladder and a small basis set can yield reasonable SSCCs that are usually attributed, lacking another explanation, to “error compensation”. In this work, solvent effects on the SSCCs are treated using the polarizable continuous model (PCM)64−66 implemented in the Gaussian program.46 The PCM model considers the molecule inside a cavity with the shape of the molecule and in a dielectric medium. Polarization effects are introduced through charges on the cavity surface. Here, the effects are estimated using the B972 functional and within three different approaches of the model: (i) the PCM model is applied to calculate the SSCCs using the geometries partially optimized in the gas phase at level B3LYP/6-31G(d,p); (ii) SSCCs are calculated on the partially optimized geometries including the

4. RESULTS AND DISCUSSION 4.1. Molecular Model and Geometry. We need appropriate geometries to calculate SSCCs that resemble, as much as possible, those shown in proteins. That is, how to optimize peptide geometries to obtain the Karplus coefficients that would be valid to apply them into usual proteins conformations. First, we have to cut the peptide chain to a size that allows one to calculate the SSCCs accuracy. In other words, the calculated SSCC values must include, as much as possible, substituent and local geometry effects. Considering that SSCC is a local property, we estimate that the dipeptide fragment shown in Figure 1 should be appropriate to study the side-chain SSCCs. Obviously, a full optimization of our model of Ala, which lacks long distance noncovalent interactions, gives a predicted full optimized geometry that is far from the usual secondary structures of α-helix or β-sheet. Although it is not the purpose of this work to make an in-depth study of the optimized geometries, many studies have been conducted recently,16,69,70 we want to show some results for comparison. Figure 2 shows the ϕ−ψ conformational space calculated at the B3LYP/6-31G(d,p) level with 144 grid points fitted by a double Fourier series. The resulting potential energy surface (PES) is topologically similar to that previously obtained.69,70 From this PES we fully optimize three local minima shown in Table 2 and marked approximately in Figure 2. Unfortunately, as was already detected,69,70 the α helical conformation is not a local minimum within this reduced model. As a reasonable starting point, we can consider folding proteins in the two most common backbone conformations: α4255


Article


Figure 3. Largest dihedral angle deviations involved in the Cα−Cβ bond against the tetrahedral angles.

geometrical parameters. In the Ala, the side-chain substituent is a CH3 group, and obviously when studying the angular dependence of side-chain SSCCs, the election of one of the three protons should give, by symmetry, identical results in a torsional property. Nevertheless, the election of one hydrogen to define the χ1 angle introduces an artificial breaking in the symmetry and the other two hydrogen atoms present slightly different behavior. Consequently, slightly different Fourier coefficients for each hydrogen are obtained. A possibility to address this problem is to treat the methyl group as a rigid symmetric rotor with local C3v symmetry (see Table S2). However, the results will be not realistic73 because the neighbor group to which the methyl group is attached, the CαN′C′Hα group, does not have 3-fold symmetry. On the other hand, for other AA residues, all except glycine and alanine, it is necessary to fix the χ1 angle (N′−Cα−Cβ−Xγ), choosing an Xγ substituent. This choice will introduce similar artifacts into the Fourier coefficients, and furthermore, in such cases the rotating group has no local C3v symmetry. That is, the artifacts created when one chose a torsional χ1 angle can be considered intrinsic to the theoretical study of the angular dependence of SSCCs and must be evaluated. 4.2. Fourier Series. In this work, the effect on the Fourier coefficients for using a specific γ-substituent is evaluated. Table S3 shows the Fourier coefficients obtained at the B3LYP/augcc-pVTZ-J level for the α-conformer using different sets of SSCCs. Table 3 sums up the main results showing the rmsd and awrmsd values. Let us only discuss the 3JHH SSCCs which are the biggest. For the other two SSCCs, 3JNH and 3JCH, the conclusions are similar. First, the results of 12 calculations, those obtained when driving χ1 from 0 to 330° at 30° increments, are shown. Here, the hydrogens of the methyl group are considered different (sets s1, s2, and s3). Although by symmetry they should be equivalent, some differences between the smaller coefficients, specially C1, C5, S3, and S5 are observed. rmsd between these three sets are 0.34 (s1−s2), 0.35 (s1−s3), and 0.22 Hz (s2−s3). Set s4 has been obtained from four calculations (χ1 = 0, 30, 60, and 90°). From these four calculations and considering equivalent the three β-hydrogens, 12 SSCCs are obtained which will give us 12 Fourier coefficients. rmsd between these coefficients and those of s1, s2, and s3 are 0.31, 0.26, and 0.25 Hz, respectively. We can also obtain six Fourier coefficients from six calculations (sets s5, s6, and s7) or from two calculations (set s8). In the first case, six 3JHH SSCCs are obtained for each proton in the sidechain group. In the second one, six SSCCs are obtained considering equivalent the three β-hydrogens. It is important to note the Fourier coefficients for these four sets are similar. Moreover, we observe that the C61 coefficient obtained in these

Figure 2. Potential energy surface for Ala calculated at the B3LYP/631G(d,p) level. Three local minima (A, B, and C) fully optimized are represented (see Table 2).

Table 2. ϕ and ψ Values (in Degree) For Three Local Minima Calculated at B3LYP/6-31G(d,p) for the Ala Model conformer

ϕa

ψa

total energy (au)

A B C

−81.39 (−81.38) −160.46 (−157.61) 72.21 (74.89)

72.66 (78.46) 168.27 (162.02) −57.67 (−70.08)

−417.237437 −417.235378 −417.233714

a

Values in parentheses were calculated at the RMP2/6-311G** in ref 69.

helix or β-sheet. Thus, the ϕ and ψ angles will be constrained in values close to those of α (−64° and −44°, respectively) and β (−121° and 128°, respectively) structures.71,72 Additionally, the side-chain group is rotated around the χ1 angle (Figure 1). This χ1 angle was rotated at 30° intervals to cover the full rotation. Thus, in the general case, 12 geometries are obtained. The effect of taking a 60° interval and the consideration of the symmetry of the rotating group is also studied. The number of optimizations can be reduced due to the symmetry of the side-chain group in Ala. In protein studies, the dihedral angles θ between the substituents around the Cα−Cβ bond are related to χ1 by the equation θ = χ1 + Δθ, where the phase shift Δθ is usually taken as 0, 120, or −120°.35 In our optimized geometries, we observe some systematic deviations between the dihedral angles optimized and those calculated using the previous relationship. These deviations, θcalc − (χ1 + Δθ), are negative and average up to 5 degrees for the dihedral angles Hβi−Cβ−Cα−C′ (i = 1, 2, and 3) for the α-conformer (see Figure 3 and Table S1). We propose, specially for these dihedral angles, to use the Δθ of −5, 115, and −125° (Table S1) instead of the usual 0, 120, and −120° (Table S1). For the β-conformer, the deviations are smaller because of its less strained conformation and the tetrahedral Δθ values can be used as well as for the remaining dihedral angles in the α-conformer. Although these corrections of five degrees on the dihedral angles will give small deviations on the predicted SSCCs derived from the Karplus equations, the consideration of these Δθ values will avoid additional errors. Partial optimized geometries were obtained constraining the ϕ, ψ, and χ1 torsional angles and relaxing the remaining 4256


Article

Journal of Chemical Theory and Computation Table 3. rmsd (Hz, eq 5) and awrmsd (%, eq 6) between the Indicated Sets and s8 Coefficients (see Table S3)a 3

JHH JCH 3 JNH awrmsd 3

s1

s2

s3

s4

s5

s6

s7

Karplusb

0.19 0.10 0.14 3.3

0.27 0.13 0.13 3.8

0.18 0.19 0.11 3.6

0.25 0.14 0.09 3.3

0.01 0.01 0.01 0.3

0.01 0.00 0.00 0.1

0.01 0.01 0.01 0.2

1.10 1.02 0.10 15.2

a Sets s1 to s3 are derived from 12 SSCCs obtained from 12 calculations (χ1 = 0, 30, 60, ..., 300, and 330); s4 is derived from 12 values obtained from 4 calculations (χ1 = 0, 30, 60, and 90); sets s5 to s7 are derived from 6 values obtained from 6 calculations (χ1 = 0, 60, 120, 180, 240, and 300); and s8 is derived from 6 values obtained from 2 calculations (χ1 = 0 and 60). bs8 coefficients neglecting C3, S1, and S2 ones, i.e., considering the classical Karplus equation.

12 sets is close to C12 1 + C5 obtained in the sets s1 to s4. 6 12 6 12 12 Relationships S1 ≈ S1 − S12 5 and S2 ≈ S2 − S4 are also obtained (see eq 3). In order to save computational resources and considering that the differences between the results of two and four calculations are small, in the rest of the paper we only consider results from two calculations. rmsd between set s8 and s1 to s4 are 0.19, 0.27, 0.18, and 0.25 Hz, respectively (see Table 3). Therefore, 0.25 Hz can be considered the error produced when instead of 12 only 6 Fourier coefficients are used for predicting 3JHH in Ala. Similarly, maximum errors of 0.19 and 0.09 Hz are obtained for 3JNH and 3JCH. awrmsd values for these sets amount to 3 or 4% (Table 3). Empirically Karplus equations consider only three coefficients (C0 to C2). Last column of Table 3 shows the rmsd and awrmsd between the set s8 and a set where the non-Karplus coefficients (C3, S1, and S2) are eliminated. Because of the large values of S2, rmsd for 3JHH and 3JCH are large and also the awrmsd amount is 15%. The dipeptide model used in this work is the smallest that includes α- and β-substituent effects27 within the vicinal SSCCs. To test the differences with a larger model, we carried out the calculations of a model that includes two additional terminal methyl groups40 at the B3LYP/aug-cc-pVTZ-J level. The rmsd and awrmsd between the Fourier coefficients obtained using the model with and without the terminal methyl groups are summarized in Table 4. The full set of

is appropriate for reducing the computational cost. However, it should be noted that large interaction effects, as hydrogen bonds, are not included within the dipeptide model. 4.3. Empirical Values. Experimental vicinal SSCCs to methyl groups 3Jmethyl (with X = H, C, and N), which is our XH case in Ala, are average-weighed on the conformer density distribution.67 In the rotational isomeric state approximation,67 the experimental SSCCs are related by 3 methyl J XH

3

0.09 0.38

0.06 0.16

JHH

α-conf β-conf

3

JCH

JNH

awrmsd

0.03 0.05

1.3 3.8

(7)

At high temperature, where the density distribution is the same for all possible conformers, the relation reduces to67 3 methyl J XH

= C0

(8)

Table 5 compares the C0 and C0 − C3 coefficients calculated using the SOPPA(CCSD)/aug-cc-pVTZ-J level and the empirical ones. The average empirical values correspond to those obtained for Ala residues in Desulfovibrio vulgaris flavodoxin.35 Theoretical results agree with the empirical ones, albeit, some deviations are observed. Theoretical results, in magnitude, are overestimated for 3JHH (between 6 and 19%), 3JNH (between 2 and 24%), and for 3JCH in the α conformer (between 1 and 4%), whereas 3JCH in the β conformer is slightly underestimated (between 3 and 8%). See Table S5 for additional data. In this work, we chose the two most usual conformations for the peptide backbone: α-helix and β-sheet. Results for both conformations are compared in Table 6. Significant differences are observed. Thus, at the SOPPA(CCSD)/aug-cc-pVTZ-J level, rmsd between α-helix and β-sheet SSCCs amount to 0.73, 0.54, and 0.44 Hz for 3JHH, 3JCH, and 3JNH, respectively. Since these values are not insignificant and in order to use Karplus equations for side-chain conformational studies, it should be more accurate to use, when possible, Fourier coefficients obtained for a main-chain conformer close to that studied. Considering the C0 Fourier coefficient, which corresponds to the average SSCC value, those rmsd represent 11, 14, and 20% of error, respectively. The three larger

Table 4. rmsd (Hz) and awrmsd (%) between the Fourier Coefficients Obtained Using the Ala Molecular Model with and without Methyl Terminal Groupsa 3

= C0 − C3 + C6

a

Calculated at the B3LYP/aug-cc-pVTZ-J level.

calculated Fourier coefficients are detailed in Table S4. Only the rmsd values of 0.4 and 0.16 Hz for 3JHH and 3JCH in the βconformer are highlighted due to a reduction of the Hα−Cα− N′ angle caused by the additional methyl groups. The remaining rmsd are smaller than 0.1 Hz. The model proposed

Table 5. C0 and C0 − C3 Coefficients (Hz) Calculated with SOPPA(CCSD)/aug-cc-pVTZ-J and Empirical Values35a 3

3

JHH

α-conf β-conf empiricalb

3

JCH

JNH

C0

C0 − C3

C0

C0 − C3

C0

C0 − C3

6.70 (14%) 6.22 (6%) 5.88 (3%)

6.99 (19%) 6.49 (10%)

3.92 (1%) 3.65 (8%) 3.96 (4%)

4.10 (4%) 3.84 (3%)

−3.14 (20%) −2.67 (2%) −2.61 (9%)

−3.23 (24%) −2.79 (7%)

a

Relative deviations against empirical values are given between parentheses. parentheses.35 4257

b

Empirical average values and relative uncertainty between DOI: 10.1021/acs.jctc.9b00131 J. Chem. Theory Comput. 2019, 15, 4252−4263

Article


Table 6. Fourier Coefficients (Hz) Calculated at the SOPPA(CCSD)/aug-cc-pVTZ-J Level for the α-(First Row) and βConformers (Second Row) SSCC 3

JHH

3

JCH

3

JNH

C0

C1

C2

C3

S1

S2

rmsda

6.70 6.22 3.92 3.65 −3.14 −2.67

−0.59 −0.63 −1.11 −1.48 0.70 0.90

6.01 5.70 3.53 3.39 −3.06 −2.59

−0.30 −0.27 −0.18 −0.19 0.09 0.12

−0.07 0.10 0.19 0.05 0.03 −0.08

1.27 0.56 −1.25 −0.74 −0.08 −0.22

0.73 0.54 0.61

rmsd (Hz) between the SSCCs of α- and β-conformer for the indicated SSCC. awrmsd amounts to 15%.

a

coefficients, C0, C2, and S2 are smaller for β-conformer than for the α-one. 4.4. WF and Basis Set Analysis. Probably our best calculated results are those of SOPPA(CCSD)/aug-cc-pVTZ-J. Theoretically, the SOPPA(CCSD) method includes more electron correlation than the other two approximations SOPPA and SOPPA(CC2). However, SOPPA(CCSD) is more expensive from a computational point of view than the other two approaches. Although the results of the three methods are usually similar (see below), the SOPPA(CCSD) seems to predict better SSCCs.39,56 rmsd and awrmsd results for the three SOPPA methods are presented in Figure 4 for the α-conformer. Values for α- and βFigure 5. rmsd (Hz) and awrmsd (%) between the results calculated with the indicated basis sets and those calculated with the aug-ccpVTZ-J one. Results obtained with the SOPPA(CCSD) method for α-conformer are shown (see Table S7 for numerical and additional data).

+G**-J basis set instead of the aug-cc-pVTZ-J one will yield results with deviations smaller than 0.2 Hz for 3JHH and 0.1 Hz for 3JCH and 3JNH. It should be noted the important reduction in the number of basis set functions (Table 1). The global awrmsd values are 5, 4, 3, and 2% for 6-31+G*-J, 6-31++G**J, 6-311+G*-J, and 6-311++G**-J results, respectively, with respect to aug-cc-pVTZ-J ones. In conclusion, the use of SOPPA(CC2) instead of SOPPA(CCSD) and the basis set 6311++G**-J (or even the 6-311+G*-J) instead of the aug-ccpVTZ-J will be appropriate to reduce the computational effort when necessary. 4.5. DFT Analysis. DFT results will be analyzed in three steps: first, the effectiveness of the basis sets with the hybrid B3LYP functional is checked; second, a set of functionals will be analyzed to detect the best performance for the three SSCCs studied; and third, the effect of the HF exchange on the SSCCs is examined. Table 7 shows rmsd between the B3LYP results and those obtained with SOPPA(CCSD). Only the results of αconformer are shown, those of the β-sheet being similar. Two types of comparisons are presented: first, B3LYP and SOPPA(CCSD) results are compared using the same basis set for both methods. It is observed that B3LYP and SOPPA(CCSD) present significant rmsd, between 0.8 and 1.1 Hz for 3 JHH, and that the best agreement (0.82, 0.39, and 0.32 Hz for 3 JHH, 3JCH, and 3JNH, respectively) is obtained for the smaller basis set (6-31+G*-J). rmsd values increase as the basis set size increases and specially when polarization and diffuse functions are added. Coupled cluster methods are very demanding on the basis set size. The results shown in Table 7 reflect deficiencies in the DFT calculations because as the basis set size increases the CCSD method more appropriately includes

Figure 4. rmsd (Hz) and awrmsd (%) between the results calculated with the indicated methods. Results calculated with the aug-cc-pVTZJ basis set for the α-conformer are shown (see Table S6 for numerical and additional data).

conformers are similar. See Table S6 for additional and numerical values. For the α-conformer, rmsd values between SOPPA(CCSD) and SOPPA(CC2) are 0.28, 0.12, and 0.06 Hz for 3JHH, 3JCH and 3JNH, respectively. The respective values between SOPPA(CCSD) and SOPPA increase to 0.37, 0.22, and 0.18 Hz. With the use of the global awrmsd parameter, the difference between SOPPA and SOPPA(CC2) and between SOPPA(CC2) and SOPPA(CCSD) is 3% and that between SOPPA and SOPPA(CCSD) is 6%. Figure 5 shows rmsd and awrmsd between results obtained with different basis sets and with the SOPPA(CCSD) approach. See Table S7 for additional and numerical values. rmsd between the aug-cc-pVTZ-J results and those calculated with the 6-31+G*-J, 6-31++G**-J, 6-311++G**-J, and 6-311+ +G**-J basis sets are, respectively, 0.47, 0.46, 0.19, and 0.18 Hz for 3JHH, 0.19, 0.12, 0.13, and 0.06 Hz for 3JCH, and 0.08, 0.03, 0.06, and 0.01 Hz for 3JNH. That is, using the 6-311+ 4258


Article


Table 7. rmsd (Hz) and awrmsd (%) between B3LYP and SOPPA(CCSD) Results Using the Indicated Basis Set for the αConformer 6-31+G*-J

6-31++G**-J

6-311+G*-J

6-311++G**-J

aug-cc-pVTZ-J

0.82 0.86

1.05 0.97

0.86 0.91

1.04 1.01

1.10 1.10

0.39 0.55

0.50 0.58

0.41 0.52

0.49 0.53

0.51 0.51

0.46 0.56

0.53 0.57

0.47 0.55

0.54 0.55

0.54 0.54

12 14

15 15

13 15

15 15

16 15

3

JHH B3LYP vs SOPPA(CCSD)a B3LYPb vs SOPPA(CCSD)/aug-cc-pVTZ 3 JCH B3LYP vs SOPPA(CCSD)a B3LYPb vs SOPPA(CCSD)/aug-cc-pVTZ 3 JNH B3LYP vs SOPPA(CCSD)a B3LYPb vs SOPPA(CCSD)/aug-cc-pVTZ awrmsd B3LYP vs SOPPA(CCSD)a B3LYPb vs SOPPA(CCSD)/aug-cc-pVTZ a

Both B3LYP and SOPPA(CCSD) using the same basis set, that indicated in the column labels. bB3LYP uses the basis set indicated in the column name, while SOPPA(CCSD) uses the aug-cc-pVTZ-J basis set.

Figure 6. awrmsd (%) between the results obtained with the indicated functional and those of SOPPA(CCSD).

exchange. Although the results depend on the type of SSCCs and to a lesser extent on the conformer studied, the best results, those with lower awrmsd, are those of the B972 functional. For 3JHH the best rmsd are those of the B972 functional (0.11 and 0.12 Hz for α and β-conformer, respectively). For 3JCH, they are those B97D functional (0.18 and 0.23 Hz), while for 3JNH, it is that of S55VWN5 (0.21 Hz) for α and B972 (0.13 Hz) for β-conformer. It should be noted that the functional S55VWN538 corresponds to SVWN5 functional where 55% of DF exchange is replaced by HF exchange. Considering Figure 6, the effect of the HF exchange on these vicinal SSCCs is not clear. For instance, the PBEPBE gives better results than the PBE1PBE which replaces 20% of DFT exchange by HF exchange. On the other hand, some functionals as S66P86 (66% of HF exchange), tHCTHhyb, S55VWN5, and the group of B97 functionals, all of them including HF exchange in same way, present mixed results. Some recent works38,55,56 showed that the HF exchange is important to produce, at least, direct (one bond) and geminal SSCCs and that the best results were obtained with high fractions of HF exchange. With an HF exchange of 47, 55, and 66%, the functionals S47LYP, S55VWN5, and S66P86 were derived. Two of these functionals yield awrmsd smaller than 10%. However, the S47LYP gives worse results than those of SLYP, with no HF exchange (see the Supporting Information). In order to analyze the HF exchange effect on the vicinal SSCCs, the results of four functionals which replace the DF exchange by HF exchange are represented in Figure 7. Table

the electronic correlation, and it is now when the differences with the B3LYP are greater. That is, DFT results resemble more those of CCSD when the basis set is small. In a second group of comparisons, SOPPA(CCSD)/aug-cc-pVTZ-J coefficients are compared with B3LYP with the five used basis sets. In other words, B3LYP results with the five used basis sets are contrasted with a fixed reference. Hence, it is observed that the basis set does not affect much the DFT results. rmsd ranges between 0.86 and 1.10 Hz for 3JHH, between 0.51 and 0.55 Hz for 3JCH, and between 0.38 and 0.41 Hz for 3JNH. Once we use tight functions to describe properly the density close to the nuclei, the size of the basis set is less important within the DF methods than with the WF ones. Considering a fixed reference, DFT results for 3JCH and 3JNH slightly improve with the basis set size. This is not observed in 3JHH. Next, we try to find out if the standard functionals can correctly reproduce these vicinal SSCCs. Figure 6 shows the rmsd and awrmsd obtained with some functionals in comparison with the SOPPA(CCSD) values (see Table S8, for additional results). The largest basis set, aug-cc-pVTZ-J, was used in these calculations. We have noted that functionals developed from the seminal B9774 give some of the best rmsd. For this reason, we test and include in Figure 6 many of these functionals (B98, B971, B972, wB97, wB97X, B97D, and wB97XD). It should be noted that these functionals add, but do not replace, some HF exchange, usually 21%, in addition to the 100% of density functional (DF) exchange. As will be indicated below, this gives results similar to those obtained with the replacement of a higher DF exchange by HF 4259


Article


Figure 7. awrmsd (%) between the results obtained with the indicated functionals replacing the percentage of DFT exchange by HF exchange (HFX %) and those obtained with SOPPA(CCSD). Results obtained using the aug-cc-pVTZ-J basis set for the α-conformer are shown.

Figure 8. awrmsd between the reference SOPPA(CCSD) results and those obtained with the B972 functional adding the indicated percentage of HF exchange (HFX %). Basis set aug-cc-pVTZ-J was used.

water. The results for both dielectric media are similar (see Tables S11 and S12 for detailed results of DMSO), therefore, only a summary of those calculated using water as a solvent are presented. As indicated above, the effect of the solvent on the SSCC may be due to a polarization of the electronic charge in the solvation, or to the geometry changes induced by solvation.75,76 Both effects are shown in Table 8, containing

S9 presents additional and numerical results obtained with these four DFs. The effect on the SSCCs depends on the type of coupling and the need for HF exchange decreases as 3JHH > 3 JCH > 3JNH. The SLYP functional presents its best results when between 0 and 20% of HF exchange is considered. However, for the other three functionals the lowest awrmsd are obtained when large HF exchange is considered. It should be noted that the results in Figure 7 are obtained replacing the DF exchange by the HF exchange. This is not a common procedure used in the B97 group of functionals that adds the HF exchange (usually 21%) without eliminating the corresponding DF exchange fraction. This way of proceeding requires a lower amount of HF exchange (see below) than when it is replaced. Thus, for the B972 functional, the lowest awrmsd is obtained by replacing about 60% HF exchange. Results presented in Figure 8 show the possibility to change the fraction of HF exchange which is added to the B972 functional. These results suggest that although the amount HF exchange depends on the type of SSCCs, the original 21% considered in the B97 group of functionals is a good value to calculate at least the vicinal SSCCs studied in this work. Table S10 presents additional and numerical results obtained adding the HF exchange to the B972 functional. 4.6. Solvent Effect. Considering the solvent effects on the calculated SSCCs is highly recommended,9,10,75 despite being quite insensitive to these effects. The solvent effect studied in this section has been computed at B972/aug-cc-pVTZ-J level because this functional it is the best balanced with respect to the accuracy of the SSCC studied and the calculation time, as we have remarked previously. The calculated SSCCs on the gas phase geometry at B3LYP/6-31G(d,p) level are taken as reference. Two dielectric media have been used, DMSO and

Table 8. rmsd (Hz) and awrmsd (%) between the Fourier Coefficients Calculated Considering the Implicit Solvent Effect (PCM, Water) and Those of B972/aug-cc-pVTZ-J (Gas Phase)a method b

PCM-GP GP-PCMc PCM-PCMd

3

JHH

0.06/0.06 0.15/0.03 0.17/0.08

3

JCH

0.18/0.03 0.05/0.03 0.22/0.07

3

JNH

awrmsd

0.04/0.19 0.06/0.07 0.08/0.18

2.3/2.6 1.9/1.2 3.6/2.9

α/β conformer results are indicated. bPCM coupling constants using the gas-phase geometry. cGas phase coupling constants using PCM geometries. dPCM coupling constants using PCM geometries.

a

the deviations (rmsd and awrmsd) using the gas phase as a reference for three different calculations: a) those obtained with the original gas phase geometry, that is, the medium only affects the SSCC; b) using partially optimized geometry with medium effect (water solvent) but not during SSCC calculations; c) using the mean effect for geometry and for SSCC calculations. From Table 8, we note that the water solvent effects as calculated using PCM in this work are small, always presenting an awrmsd lower than 4%. The same results are obtained for DMSO solvent, showing that the optimized structures 4260


Article

Journal of Chemical Theory and Computation obtained for gas phase and PCM present negligible changes in bond distances and angles for this alanine dipetide model, particularly for the β conformer. These results can be confirmed analyzing the Fourier coefficients (see Table S12). The C0 coefficients, that can be compared with the empirical value of Ala residues in AA, increases with both solvents less than 0.2 Hz in magnitude for all the SSCCs. Our results show that the influence of a dielectric medium is small and constant because the SSCCs are less affected by changes induced by solvation in the molecular geometry.75,76

5. CONCLUSIONS A deep study of Karplus extended equations has been carried out. The number of Fourier coefficients to be included in the Karplus equation and how to calculate them is analyzed. Theoretical WF and DF methods are used to calculate the SSCCs. From a theoretical point of view, a Karplus equation including Fourier coefficients up to C3 and up to S2 seems to predict reasonably the calculated SSCCs. S3 and S1 usually not included in empirical Karplus equations are small. However, the calculated S2 coefficient is not negligible and its incorporation within the Karplus equation will be significant. In order to get appropriate theoretical Fourier coefficients for the side-chain SSCCs, it should be more accurate to use the correct main-chain conformation. This work shows that the relative differences between the SSCCs for α- and βconformers, which are not negligible, amount up to 15%. In this work, the best used WF method is that of SOPPA(CCSD) and the largest basis set is the aug-ccpVTZ-J. Comparing with the remaining methods and basis sets, it seems that the SOPPA(CC2) and the basis set 6-311+ +G**-J (or even the 6-311+G*-J) can be used instead of the SOPPA(CCSD)/aug-cc-pVTZ-J to reduce the computational effort. The SOPPA(CC2)/6-311+G*-J will deteriorate the SSCCs around 5% with respect to those obtained with the SOPPA(CCSD)/aug-cc-pVTZ-J level. The best functional for the three analyzed vicinal SSCCs is the B972. The HF exchange is important for these SSCCs. Concerning the HF exchange, the best results are obtained either adding around 21% to the DF exchange or replacing a larger amount (around 55%) of the DF exchange. Using theoretical Karplus equations, rather than empirical ones for side-chain conformational studies, are promising relationships considering that specific Fourier coefficients can be obtained for specific backbone conformations and for specific substituent models. In summary, the computational protocol presented in this paper provides a rigorous and general approach for evaluating vicinal SSCCs in AAs. In Figure 9, we summarized more significant awrmsd (11 values) between the pair of results. At the top of this Figure 9, two first awrmsd represent the influence of the number of Karplus coefficients and we can conclude that the use of three coefficients (Karplus classical equation) is inadequate showing an awrmsd around 15%. However, 6 coefficients are a good alternative versus 12 ones. The introduction of terminal methyl groups affects scarcely both conformers, but the use of a β conformer produces more differences for the SSCCs. Within the comparison of methods and basis sets, Figure 9 shows that SOPPA (CC2) is the key quantity-quality for computing vicinal coupling constants and that the best functional for this type of calculation produces results at the same level as SOPPA calculations, whether WF or DFT methods are used. We conclude that Pople style basis sets for calculations of SSCCs show a quality of results similar to

Figure 9. Summary of main awrmsd (%) between pairs of results.

larger Dunning style basis sets, reducing the contracted basis functions to the half number ones. Finally, the solvent effect on these SSCCs evaluated using the PCM model are small regardless of whether the solvent directly affects the SSCCs or indirectly through changes in geometry. Future work extending this computational protocol to the rest of the AAs is under way, and we will use the conclusions extracted here for its application.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.9b00131. Additional and numerical data used in the summary tables and figures of this paper (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

J. M. García de la Vega: 0000-0002-9676-7534 Funding

The authors acknowledge the financial support from Project LIQUORGAS-S2013/MAE-2800. Computer time provided by the Centro de Computación Cientifica of Universidad Autónoma de Madrid is gratefully acknowledged. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Berg, J. M.; Tymoczko, J. L.; Stryer, L. Biochemistry; W. H. Freeman: New York, 2002; Chapter Protein Structure and Function. (2) Wright, P. E.; Dyson, H. J. Intrinsically Unstructured Proteins: Re-assessing the Protein Structure-Function Paradigm. J. Mol. Biol. 1999, 293, 321−331. (3) Yang, Y.; Gao, J.; Wang, J.; Heffernan, R.; Hanson, J.; Paliwal, K.; Zhou, Y. Sixty-five years of the long march in protein secondary structure prediction: the final stretch? Briefings Bioinf. 2018, 19, 482− 494. (4) Wuthrich, K. NMR of Proteins and Nucleic Acids; John Wiley & Sons, Inc., 1986. 4261


Article


(26) Tuttle, L. M.; Dyson, H. J.; Wright, P. E. Side-Chain Conformational Heterogeneity of Intermediates in the Escherichia coli Dihydrofolate Reductase Catalytic Cycle. Biochemistry 2013, 52, 3464−3477. (27) Díez, E.; San Fabián, J.; Guilleme, J.; Altona, C.; Donders, L. A. Vicinal proton-proton coupling constants. I. formulation of an equation including interactions between substituents. Mol. Phys. 1989, 68, 49−63. (28) Pachler, K. G. R. Extended Hückel Theory MO Calculations of Proton-Proton Coupling Constants. The Substituent Effect in Fluoroethane. Tetrahedron Lett. 1970, 11, 1955−1958. (29) Donders, A.; De Leeuw, F. A. A. M.; Altona, C. Relationship between Proton−Proton NMR Coupling-Constants and Substituent Electronegativities. IV−An Extended Karplus Equation Accounting for Interactions between Substituents and its Application to CouplingConstant data Calculated by the Extended Hückel-method. Magn. Reson. Chem. 1989, 27, 556−563. (30) Guilleme, J.; San Fabián, J.; Casanueva, J.; Díez, E. Vicinal proton-proton coupling constants: MCSCF ab initio calculations of ethane. Chem. Phys. Lett. 1999, 314, 168−175. (31) Suardíaz, R.; Crespo-Otero, R.; Pérez, C.; San Fabián, J.; García de la Vega, J. M. Communication: Accurate determination of sidechain torsion angle χ1 in proteins: Phenylalanine residues. J. Chem. Phys. 2011, 134, 061101. (32) Wuthrich, K. NMR structures of biological macromolecules. Magn. Reson. Chem. 2003, 41, S89−S98. (33) Hu, J.-S.; Bax, A. Measurement of Three-Bond 13C-13C J Couplings between Carbonyl and Carbonyl/Carboxyl Carbons in Isotopically Enriched Proteins. J. Am. Chem. Soc. 1996, 118, 8170− 8171. (34) Hu, J.-S.; Bax, A. Determination of ϕ and χ(1) angles in proteins from C-13-C-13 three-bond J couplings measured by threedimensional heteronuclear NMR. How planar is the peptide bond? J. Am. Chem. Soc. 1997, 119, 6360−6368. (35) Pérez, C.; Löhr, F.; Rüterjans, H.; Schmidt, J. M. SelfConsistent Karplus Parametrization of 3J Coupling Depending on the Polypeptide Side-Chain Torsion χ1. J. Am. Chem. Soc. 2001, 123, 7081−7093. (36) Geertsen, J.; Oddershede, J.; Raynes, W. T.; Scuseria, G. E. Nuclear Spin-Spin Coupling in the Methane Isotopomers. J. Magn. Reson. 1991, 93, 458−471. (37) Guilleme, J.; San Fabián, J. Basis sets and active space in multiconfigurational self-consistent field calculations of nuclear magnetic resonance spin-spin coupling constants. J. Chem. Phys. 1998, 109, 8168−8181. (38) San Fabián, J.; García de la Vega, J. M.; San Fabián, E. Improvements in DFT Calculations of Spin-Spin Coupling Constants. J. Chem. Theory Comput. 2014, 10, 4938−4949. (39) García de la Vega, J. M.; San Fabián, J. Assessment of DFT functionals with fluorine-fluorine coupling constants. Mol. Phys. 2015, 113, 1924−1936. (40) Hermans, J. The amino acid dipeptide: Small but still influential after 50 years. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 3095−3096. (41) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648−5652. (42) Lee, C.; Yang, W.; Parr, R. G. Development of the ColleSalvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (43) Hariharan, P.; Pople, J. The effect of d-functions on molecular orbital energies for hydrocarbons. Chem. Phys. Lett. 1972, 16, 217− 219. (44) Hariharan, P. C.; Pople, J. A. Influence of polarization functions on molecular-orbital hydrogenation energies. Theor. Chem. Acc. 1973, 28, 213−222. (45) Francl, M. M.; Pietro, W. J.; Hehre, W. J.; Binkley, J. S.; Gordon, M. S.; DeFrees, D. J.; Pople, J. Self-Consistent Molecular Orbital Methods. 23. A polarization-type basis set for 2nd-row elements. J. Chem. Phys. 1982, 77, 3654−3665.

(5) Cavanagh, J.; Fairbrother, W. J.; III, A. G. P; Rance, M.; Skelton, N. J. Protein NMR Spectroscopy. Principles and Practice; Elsevier Academic Press, 2006. (6) Carvalho, A. T. P.; Barrozo, A.; Doron, D.; Kilshtain, A. V.; Major, D. T.; Kamerlin, S. C. L. Challenges in Computational Studies of Enzyme Structure, Function and Dynamics. J. Mol. Graphics Modell. 2014, 54, 62−79. (7) Mantsyzov, A. B.; Shen, Y.; Lee, J. H.; Hummer, G.; Bax, A. MERA: a webserver for evaluating backbone torsion angle distributions in dynamic and disordered proteins from NMR data. J. Biomol. NMR 2015, 63, 85−95. (8) Lee, J. H.; Ying, J.; Bax, A. Quantitative evaluation of positive ϕ angle propensity in flexible regions of proteins from three-bond J couplings. Phys. Chem. Chem. Phys. 2016, 18, 5759−5770. (9) Krivdin, L. B. Carbon-carbon spin-spin coupling constants: Practical applications of theoretical calculations. Prog. Nucl. Magn. Reson. Spectrosc. 2018, 105, 54−99. (10) Krivdin, L. B. Theoretical calculations of carbon-hydrogen spinspin coupling constants. Prog. Nucl. Magn. Reson. Spectrosc. 2018, 108, 17−73. (11) Bystrov, V. F. Spin-spin coupling and the conformational states of peptide systems. Prog. Nucl. Magn. Reson. Spectrosc. 1976, 10, 41− 81. (12) Salvador, P.; Tsai, I.-H. M.; Dannenberg, J. J. J-coupling constants for a trialanine peptide a a function of dihedral angles calculated by density functional theory over the full Ramachandran space. Phys. Chem. Chem. Phys. 2011, 13, 17484−17493. (13) Salvador, P. Dependencies of J-Couplings upon Dihedral Angles on Proteins. Annu. Rep. NMR Spectrosc. 2014, 81, 185−227. (14) Arai, M.; Sugase, K.; Dyson, H. J.; Wright, P. E. Conformational propensities of intrinsically disordered proteins influence the mechanism of binding and folding. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, 9614−9619. (15) Bhowmick, A.; Brookes, D. H.; Yost, S. R.; Dyson, H. J.; Forman-Kay, J. D.; Gunter, D.; Head-Gordon, M.; Hura, G. L.; Pande, V. S.; Wemmer, D. E.; Wright, P. E.; Head-Gordon, T. Finding our way in the dark proteome. J. Am. Chem. Soc. 2016, 138, 9730− 9742. (16) Choi, J.-M.; Pappu, R. V. Experimentally Derived and Computationally Optimized Backbone Conformational Statistics for Blocked Amino Acids. J. Chem. Theory Comput. 2019, 15, 1355−1366. (17) Karplus, M. Contact Electro-Spin Coupling of Nuclear Magnetic Moments. J. Chem. Phys. 1959, 30, 11−15. (18) Karplus, M. Vicinal Proton Coupling in Nuclear Magnetic Resonance. J. Am. Chem. Soc. 1963, 85, 2870−2871. (19) Wang, A. C.; Bax, A. Reparametrization of the Karplus Relation for 3J(Hα − N) and 3J(HN − C′) in Peptides from Uniformly 13C/15N Enriched Human Ubiquitin. J. Am. Chem. Soc. 1995, 117, 1810−1813. (20) Wang, A. C.; Bax, A. Determination of the Backbone Dihedral Angles ϕ in Human Ubiquitin from Reparameterized Empirical Karplus Equations. J. Am. Chem. Soc. 1996, 118, 2483−2494. (21) Case, D. A.; Scheurer, C.; Bruschweiler, R. Static and Dynamic Effects on vicinal Scalar J Couplings in Proteins and Peptides: A MD/ DFT Analysis. J. Am. Chem. Soc. 2000, 122, 10390−10397. (22) Case, D. A. In NMR Parameters in Proteins and Nucleic Acids; Kaupp, M., Bühl, M., Malkin, V. G., Eds.; Wiley-VCH, 2004; pp 341− 351. (23) Schmidt, J. M.; Blümel, M.; Löhr, F.; Ruterjans, H. Selfconsistente 3J coupling analysis for the joint calibration of Karplus coefficients and evaluation of torsion angles. J. Biomol. NMR 1999, 14, 1−12. (24) Chou, J. J.; Case, D. A.; Bax, A. Insights into the Mobility of Methyl-Bearing Side Chains in Proteins from 3JCC and 3JCN Couplings. J. Am. Chem. Soc. 2003, 125, 8959−8966. (25) Vajpai, N.; Gentner, M.; Huang, J.-r; Blackledge, M.; Grzesiek, S. Side-Chain χ1 Conformations in Urea-Denatured Ubiquitin and Protein G from 3J Coupling Constants and Residual Dipolar Couplings. J. Am. Chem. Soc. 2010, 132, 3196−3203. 4262


Article

Journal of Chemical Theory and Computation (46) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Keith, T.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, revision D.01, Gaussian, Inc., Wallingford, CT, 2013. (47) García de la Vega, J.; San Fabián, J.; Crespo-Otero, R.; Suardíaz, R.; Pérez, C. Theoretical DFT Karplus Equations: Amino Acid Side-Chain Torsion Angle χ1. Int. J. Quantum Chem. 2013, 113, 656−660. (48) Anbarasan, R.; Sundar, J. K.; Lakshmi, M. A. Experimental and theoretical insights on the dimeric cation molecule L-alanine Lalaninium chloride for nonlinear optical applications. J. Mol. Struct. 2019, 1179, 154−160. (49) Nielsen, E.; Jørgensen, P.; Oddershede, J. Transition Moments and Dynamic Polarizabilities in a 2nd order Polarization Propagator Approach. J. Chem. Phys. 1980, 73, 6238. (50) Geertsen, J.; Oddershede, J. 2nd-Order polarization propagator calculations of indirect nuclear-spin spin coupling tensors in the water molecule. Chem. Phys. 1984, 90, 301−311. (51) Enevoldsen, T.; Oddershede, J.; Sauer, S. P. A. Correlated calculations of indirect nuclear spin-spin coupling constants using second-order polarization propagator approximations: SOPPA and SOPPA(CCSD). Theor. Chem. Acc. 1998, 100, 275−284. (52) Kjær, H.; Sauer, S. P.; Kongsted, J. Benchmarking NMR indirect nuclear spin-spin coupling constants: SOPPA, SOPPA(CC2), and SOPPA(CCSD) versus CCSD. J. Chem. Phys. 2010, 133, 144106. (53) Sauer, S. P. A. Second-order polarization propagator approximation with coupled-cluster singles and doubles amplitudes SOPPA(CCSD): the polarizability and hyperpolarizability of Li−. J. Phys. B: At., Mol. Opt. Phys. 1997, 30, 3773−3780. (54) Aidas, K.; Angeli, C.; Bak, K. L.; Bakken, V.; Bast, R.; Boman, L.; Christiansen, O.; Cimiraglia, R.; Coriani, S.; Dahle, P.; Dalskov, E. K.; Ekström, U.; Enevoldsen, T.; Eriksen, J. J.; Ettenhuber, P.; Fernández, B.; Ferrighi, L.; Fliegl, H.; Frediani, L.; Hald, K.; Halkier, A.; Hättig, C.; Heiberg, H.; Helgaker, T.; Hennum, A. C.; Hettema, H.; Hjertenæs, E.; Host, S.; Høyvik, I.-M.; Iozzi, M. F.; Jansik, B.; Jensen, H. J. A.; Jonsson, D.; Jørgensen, P.; Kauczor, J.; Kirpekar, S.; Kjærgaard, T.; Klopper, W.; Knecht, S.; Kobayashi, R.; Koch, H.; Kongsted, J.; Krapp, A.; Kristensen, K.; Ligabue, A.; Lutnæs, O. B.; Melo, J. I.; Mikkelsen, K. V.; Myhre, R. H.; Neiss, C.; Nielsen, C. B.; Norman, P.; Olsen, J.; Olsen, J. M. H.; Osted, A.; Packer, M. J.; Pawlowski, F.; Pedersen, T. B.; Provasi, P. F.; Reine, S.; Rinkevicius, Z.; Ruden, T. A.; Ruud, K.; Rybkin, V.; Salek, P.; Samson, C. C. M.; de Merás, A. S.; Saue, T.; Sauer, S. P. A.; Schimmelpfennig, B.; Sneskov, K.; Steindal, A. H.; Sylvester-Hvid, K. O.; Taylor, P. R.; Teale, A. M.; Tellgren, E. I.; Tew, D. P.; Thorvaldsen, A. J.; Thøgersen, L.; Vahtras, O.; Watson, M. A.; Wilson, D. J. D.; Ziolkowski, M.; Ågren, H. The Dalton quantum chemistry program system. WIREs Comput. Mol. Sci. 2014, 4, 269−284. (55) San Fabián, J.; Omar, S.; García de la Vega, J. M. Towards quantifying the role of exact exchange in the prediction hydrogen bond spin-spin coupling constants involving fluorine. J. Chem. Phys. 2016, 145, 084301. (56) García de la Vega, J. M.; Omar, S.; San Fabián, J. Performance of wave function and density functional methods for water hydrogen bond spin−spin coupling constants. J. Mol. Model. 2017, 23, 134.

(57) Dirac, P. A. M. Note on exchange phenomena in the Thomas atom. Math. Proc. Cambridge Philos. Soc. 1930, 26, 376−385. (58) Slater, J. C. A. Simplification of the Hartree-Fock Method. Phys. Rev. 1951, 81, 385−390. (59) Kjær, H.; Sauer, S. P. Pople Style Basis Sets for the Calculation of NMR Spin-Spin Coupling Constants: the 6-31G-J and 6-311G-J Basis Sets. J. Chem. Theory Comput. 2011, 7, 4070−4076. (60) Provasi, P. F.; Aucar, G. A.; Sauer, S. P. A. The effect of lone pairs and electronegativity on the indirect nuclear spin-spin coupling constants in CH2X (X = CH2, NH, O, S): Ab initio calculations using optimized contracted basis sets. J. Chem. Phys. 2001, 115, 1324−1334. (61) Jensen, F. The basis set convergence of spin-spin coupling constants calculated by density functional methods. J. Chem. Theory Comput. 2006, 2, 1360−1369. (62) Benedikt, U.; Auer, A. A.; Jensen, F. Optimization of augmentation functions for correlated calculations of spin-spin coupling constants and related properties. J. Chem. Phys. 2008, 129, 064111. (63) Jensen, F. The optimum contraction of basis sets for calculating spin-spin coupling constants. Theor. Chem. Acc. 2010, 126, 371−382. (64) Miertuš, S.; Scrocco, E.; Tomasi, J. Electrostatic interaction of a solute with a continuum. A direct utilization of ab initio molecular potentials for the prevision of solvent effects. Chem. Phys. 1981, 55, 117−129. (65) Cammi, R.; Tomasi, J. Remarks on the use of the apparent surface charges (ASC) methods in solvation problems: Iterative versus matrix-inversion procedures and the renormalization of the apparent charges. J. Comput. Chem. 1995, 16, 1449−1458. (66) Mennucci, B.; Cancés, E.; Tomasi, J. Evaluation of solvent effects in isotropic and anisotropic dielectrics and in ionic solutions with a unified integral equation method: Theoretical bases, computational implementation, and numerical applications. J. Phys. Chem. B 1997, 101, 10506−10517. (67) Guilleme, J.; San Fabián, J.; Diez, E.; Bermejo, F.; Esteban, A. L. Vicinal protón-protón coupling constants. II. Analysis of the effect of interaction between geminal substituents upon vicinal coupling to methyl groups. Mol. Phys. 1989, 68, 65−85. (68) Coxon, B. Developments in the Karplus equations as they relate to the NMR coupling constants of carbohydrates. Adv. Carbohydr. Chem. Biochem. 2009, 62, 17−82. (69) Yu, C.-H.; Norman, M. A.; Schäfer, L.; Ramek, M.; Peeters, A.; van Alsenoy, C. Ab initio conformational analysis of N-formyl Lalanine amide including electron correlation. J. Mol. Struct. 2001, 567−568, 361−374. (70) Mironov, V.; Alexeev, Y.; Mulligan, V. K.; Fedorov, D. G. A Systematic Study of Minima in Alanine Dipeptide. J. Comput. Chem. 2019, 40, 297−309. (71) Spera, S.; Bax, A. Empirical Correlation between Protein Backbone Conformation and Cα and Cβ 13C Nuclear Magnetic Resonance Chemical Shifts. J. Am. Chem. Soc. 1991, 113, 5490−5492. (72) Suardíaz, R.; Pérez, C.; García de la Vega, J. M.; San Fabián, J.; Contreras, R. H. Theoretical Karplus relationships for vicinal coupling constants around χ1 in Valine. Chem. Phys. Lett. 2007, 442, 119−123. (73) Mastryukov, V. S.; Samdal, S. Asymmetry in Methyl Group of Ethane During Internal Rotation: ab initio Study. J. Comput. Chem. 1998, 19, 1141−1145. (74) Becke, A. D. Density-functional thermochemistry. V. Systematic optimization of exchange-correlation functionals. J. Chem. Phys. 1997, 107, 8554−8560. (75) Ruud, K.; Frediani, L.; Cammi, R.; Mennucci, B. Solvent effects on the indirect spin-spin coupling constants of benzene: The DFTPCM approach. Int. J. Mol. Sci. 2003, 4, 119−134. (76) Astrand, P.-O.; Mikkelsen, K. V.; Jørgensen, P.; Ruud, K.; Helgaker, T. Solvent effects on nuclear shielding and spin-spin couplings of hydrogen selenide. J. Chem. Phys. 1998, 108, 2528−2537.

4263


Computational Protocol to Evaluate Side-Chain Vicinal Spin–Spin

Recommend Documents