Coupling between conformation and proton binding in proteins

13C Chemical Shifts in Proteins: A Rich Source of Encoded Structural Information

Despite the formidable progress in Nuclear Magnetic Resonance (NMR) spectroscopy, quality assessment of NMR-derived structures remains as an important problem. Thus, validation of protein structures is essential for the spectroscopists, since it could enable them to detect structural flaws and potentially guide their efforts in further refinement. Moreover, availability of accurate and efficient validation tools would help molecular biologists and computational chemists to evaluate quality of available experimental structures and to select a protein model which is the most suitable for a given scientific problem. The ¹³C^α nuclei are ubiquitous in proteins, moreover, their shieldings are easily obtainable from NMR experiments and represent a rich source of encoded structural information that makes ¹³C^α chemical shifts an attractive candidate for use in computational methods aimed at determination and validation of protein structures. In this chapter, the basis of a novel methodology of computing, at the quantum chemical level of theory, the ¹³C^α shielding for the amino acid residues in proteins is described. We also identify and examine the main factors affecting the ¹³C^α-shielding computation. Finally, we illustrate how the information encoded in the ¹³C chemical shifts can be used for a number of applications, viz., from protein structure prediction of both α-helical and β-sheet conformations, to determination of the fraction of the tautomeric forms of the imidazole ring of histidine in proteins as a function of pH or to accurate detection of structural flaws, at a residue-level, in NMR-determined protein models.

All-Atom Internal Coordinate Mechanics (ICM) Force Field for Hexopyranoses and Glycoproteins

We present an extension of the all-atom internal-coordinate force field, ICMFF, that allows for simulation of heterogeneous systems including hexopyranose saccharides and glycan chains in addition to proteins. A library of standard glycan geometries containing α- and β-anomers of the most common hexapyranoses, i.e., d-galactose, d-glucose, d-mannose, d-xylose, l-fucose, N-acetylglucosamine, N-acetylgalactosamine, sialic, and glucuronic acids, is created based on the analysis of the saccharide structures reported in the Cambridge Structural Database. The new force field parameters include molecular electrostatic potential-derived partial atomic charges and the torsional parameters derived from quantum mechanical data for a collection of minimal molecular fragments and related molecules. The ϕ/ψ torsional parameters for different types of glycosidic linkages are developed using model compounds containing the key atoms in the full carbohydrates, i.e., glycosidic-linked tetrahydropyran-cyclohexane dimers. Target data for parameter optimization include two-dimensional energy surfaces corresponding to the ϕ/ψ glycosidic dihedral angles in the disaccharide analogues, as determined by quantum mechanical MP2/6-31G** single-point energies on HF/6-31G** optimized structures. To achieve better agreement with the observed geometries of glycosidic linkages, the bond angles at the O-linkage atoms are added to the internal variable set and the corresponding bond bending energy term is parametrized using quantum mechanical data. The resulting force field is validated on glycan chains of 1-12 residues from a set of high-resolution X-ray glycoprotein structures based on heavy atom root-mean-square deviations of the lowest-energy glycan conformations generated by the biased probability Monte Carlo (BPMC) molecular mechanics simulations from the native structures. The appropriate BPMC distributions for monosaccharide-monosaccharide and protein-glycan linkages are derived from the extensive analysis of conformational properties of glycoprotein structures reported in the Protein Data Bank. Use of the BPMC search leads to significant improvements in sampling efficiency for glycan simulations. Moreover, good agreement with the X-ray glycoprotein structures is achieved for all glycan chain lengths. Thus, average/median RMSDs are 0.81/0.68 Å for one-residue glycans and 1.32/1.47 Å for three-residue glycans. RMSD from the native structure for the lowest-energy conformation of the 12-residue glycan chain (PDB ID 3og2) is 1.53 Å. Additionally, results obtained for free short oligosaccharides using the new force field are in line with the available experimental data, i.e., the most populated conformations in solution are predicted to be the lowest energy ones. The newly developed parameters allow for the accurate modeling of linear and branched hexopyranose glycosides in heterogeneous systems.

13C Chemical Shifts in Proteins: A Rich Source of Encoded Structural Information

Despite the formidable progress in Nuclear Magnetic Resonance (NMR) spectroscopy, quality assessment of NMR-derived structures remains as an important problem. Thus, validation of protein structures is essential for the spectroscopists, since it could enable them to detect structural flaws and potentially guide their efforts in further refinement. Moreover, availability of accurate and efficient validation tools would help molecular biologists and computational chemists to evaluate quality of available experimental structures and to select a protein model which is the most suitable for a given scientific problem. The ¹³C^α nuclei are ubiquitous in proteins, moreover, their shieldings are easily obtainable from NMR experiments and represent a rich source of encoded structural information that makes ¹³C^α chemical shifts an attractive candidate for use in computational methods aimed at determination and validation of protein structures. In this chapter, the basis of a novel methodology of computing, at the quantum chemical level of theory, the ¹³C^α shielding for the amino acid residues in proteins is described. We also identify and examine the main factors affecting the ¹³C^α-shielding computation. Finally, we illustrate how the information encoded in the ¹³C chemical shifts can be used for a number of applications, viz., from protein structure prediction of both α-helical and β-sheet conformations, to determination of the fraction of the tautomeric forms of the imidazole ring of histidine in proteins as a function of pH or to accurate detection of structural flaws, at a residue-level, in NMR-determined protein models.

Potential of mean force of water-proton bath and molecular dynamic simulation of proteins at constant pH

An advanced implicit solvent model of water-proton bath for protein simulations at constant pH is presented. The implicit water-proton bath model approximates the potential of mean force of a protein in water solvent in a presence of hydrogen ions. Accurate and fast computational implementation of the implicit water-proton bath model is developed using the continuum electrostatic Poisson equation model for calculation of ionization equilibrium and the corrected MSR6 generalized Born model for calculation of the electrostatic atom-atom interactions and forces. Molecular dynamics (MD) method for protein simulation in the potential of mean force of water-proton bath is developed and tested on three proteins. The model allows to run MD simulations of proteins at constant pH, to calculate pH-dependent properties and free energies of protein conformations. The obtained results indicate that the developed implicit model of water-proton bath provides an efficient way to study thermodynamics of biomolecular systems as a function of pH, pH-dependent ionization-conformation coupling, and proton transfer events.

Development of a new physics-based internal coordinate mechanics force field and its application to protein loop modeling

We report the development of internal coordinate mechanics force field (ICMFF), new force field parameterized using a combination of experimental data for crystals of small molecules and quantum mechanics calculations. The main features of ICMFF include: (a) parameterization for the dielectric constant relevant to the condensed state (ε = 2) instead of vacuum, (b) an improved description of hydrogen-bond interactions using duplicate sets of van der Waals parameters for heavy atom-hydrogen interactions, and (c) improved backbone covalent geometry and energetics achieved using novel backbone torsional potentials and inclusion of the bond angles at the C(α) atoms into the internal variable set. The performance of ICMFF was evaluated through loop modeling simulations for 4-13 residue loops. ICMFF was combined with a solvent-accessible surface area solvation model optimized using a large set of loop decoys. Conformational sampling was carried out using the biased probability Monte Carlo method. Average/median backbone root-mean-square deviations of the lowest energy conformations from the native structures were 0.25/0.21 Å for four residues loops, 0.84/0.46 Å for eight residue loops, and 1.16/0.73 Å for 12 residue loops. To our knowledge, these results are significantly better than or comparable with those reported to date for any loop modeling method that does not take crystal packing into account. Moreover, the accuracy of our method is on par with the best previously reported results obtained considering the crystal environment. We attribute this success to the high accuracy of the new ICM force field achieved by meticulous parameterization, to the optimized solvent model, and the efficiency of the search method.

Advances in implicit models of water solvent to compute conformational free energy and molecular dynamics of proteins at constant pH

Modern implicit solvent models for macromolecular simulations in water-proton bath are considered. The fundamental quantity that implicit models approximate is the solute potential of mean force, which is obtained by averaging over solvent degrees of freedom. The implicit solvent models suggest practical ways to calculate free energies of macromolecular conformations taking into account equilibrium interactions with water solvent and proton bath, while the explicit solvent approach is unable to do that due to the need to account for a large number of solvent degrees of freedom. The most advanced realizations of the implicit continuum models by different research groups are discussed, their accuracy are examined, and some applications of the implicit solvent models to macromolecular modeling, such as free energy calculations, protein folding, and constant pH molecular dynamics are highlighted.

FAMBE-pH: a fast and accurate method to compute the total solvation free energies of proteins

A fast and accurate method to compute the total solvation free energies of proteins as a function of pH is presented. The method makes use of a combination of approaches, some of which have already appeared in the literature; (i) the Poisson equation is solved with an optimized fast adaptive multigrid boundary element (FAMBE) method; (ii) the electrostatic free energies of the ionizable sites are calculated for their neutral and charged states by using a detailed model of atomic charges; (iii) a set of optimal atomic radii is used to define a precise dielectric surface interface; (iv) a multilevel adaptive tessellation of this dielectric surface interface is achieved by using multisized boundary elements; and (v) 1:1 salt effects are included. The equilibrium proton binding/release is calculated with the Tanford-Schellman integral if the proteins contain more than approximately 20-25 ionizable groups; for a smaller number of ionizable groups, the ionization partition function is calculated directly. The FAMBE method is tested as a function of pH (FAMBE-pH) with three proteins, namely, bovine pancreatic trypsin inhibitor (BPTI), hen egg white lysozyme (HEWL), and bovine pancreatic ribonuclease A (RNaseA). The results are (a) the FAMBE-pH method reproduces the observed pK a's of the ionizable groups of these proteins within an average absolute value of 0.4 p K units and a maximum error of 1.2 p K units and (b) comparison of the calculated total pH-dependent solvation free energy for BPTI, between the exact calculation of the ionization partition function and the Tanford-Schellman integral method, shows agreement within 1.2 kcal/mol. These results indicate that calculation of total solvation free energies with the FAMBE-pH method can provide an accurate prediction of protein conformational stability at a given fixed pH and, if coupled with molecular mechanics or molecular dynamics methods, can also be used for more realistic studies of protein folding, unfolding, and dynamics, as a function of pH.

Use of decoys to optimize an all-atom force field including hydration

A novel method of parameter optimization is proposed. It makes use of large sets of decoys generated for six nonhomologous proteins with different architecture. Parameter optimization is achieved by creating a free energy gap between sets of nativelike and nonnative conformations. The method is applied to optimize the parameters of a physics-based scoring function consisting of the all-atom ECEPP05 force field coupled with an implicit solvent model (a solvent-accessible surface area model). The optimized force field is able to discriminate near-native from nonnative conformations of the six training proteins when used either for local energy minimization or for short Monte Carlo simulated annealing runs after local energy minimization. The resulting force field is validated with an independent set of six nonhomologous proteins, and appears to be transferable to proteins not included in the optimization; i.e., for five out of the six test proteins, decoys with 1.7- to 4.0-A all-heavy-atom root mean-square deviations emerge as those with the lowest energy. In addition, we examined the set of misfolded structures created by Park and Levitt using a four-state reduced model. The results from these additional calculations confirm the good discriminative ability of the optimized force field obtained with our decoy sets.

Factors affecting the use of 13C(alpha) chemical shifts to determine, refine, and validate protein structures

Interest centers here on the analysis of two different, but related, phenomena that affect side-chain conformations and consequently 13C(alpha) chemical shifts and their applications to determine, refine, and validate protein structures. The first is whether 13C(alpha) chemical shifts, computed at the DFT level of approximation with charged residues is a better approximation of observed 13C(alpha) chemical shifts than those computed with neutral residues for proteins in solution. Accurate computation of 13C(alpha) chemical shifts requires a proper representation of the charges, which might not take on integral values. For this analysis, the charges for 139 conformations of the protein ubiquitin were determined by explicit consideration of protein binding equilibria, at a given pH, that is, by exploring the 2(xi) possible ionization states of the whole molecule, with xi being the number of ionizable groups. The results of this analysis, as revealed by the shielding/deshielding of the 13C(alpha) nucleus, indicated that: (i) there is a significant difference in the computed 13C(alpha) chemical shifts, between basic and acidic groups, as a function of the degree of charge of the side chain; (ii) this difference is attributed to the distance between the ionizable groups and the 13C(alpha) nucleus, which is shorter for the acidic Asp and Glu groups as compared with that for the basic Lys and Arg groups; and (iii) the use of neutral, rather than charged, basic and acidic groups is a better approximation of the observed 13C(alpha) chemical shifts of a protein in solution. The second is how side-chain flexibility influences computed 13C(alpha) chemical shifts in an additional set of ubiquitin conformations, in which the side chains are generated from an NMR-derived structure with the backbone conformation assumed to be fixed. The 13C(alpha) chemical shift of a given amino acid residue in a protein is determined, mainly, by its own backbone and side-chain torsional angles, independent of the neighboring residues; the conformation of a given residue itself, however, depends on the environment of this residue and, hence, on the whole protein structure. As a consequence, this analysis reveals the role and impact of an accurate side-chain computation in the determination and refinement of protein conformation. The results of this analysis are: (i) a lower error between computed and observed 13C(alpha) chemical shifts (by up to 3.7 ppm), was found for approximately 68% and approximately 63% of all ionizable residues and all non-Ala/Pro/Gly residues, respectively, in the additional set of conformations, compared with results for the model from which the set was derived; and (ii) all the additional conformations exhibit a lower root-mean-square-deviation (1.97 ppm < or = rmsd < or = 2.13 ppm), between computed and observed 13C(alpha) chemical shifts, than the rmsd (2.32 ppm) computed for the starting conformation from which this additional set was derived. As a validation test, an analysis of the additional set of ubiquitin conformations, comparing computed and observed values of both 13C(alpha) chemical shifts and chi(1) torsional angles (given by the vicinal coupling constants, 3J(N-Cgamma) and 3J(C'-Cgamma), is discussed.

Assessment of two theoretical methods to estimate potentiometric titration curves of peptides: comparison with experiment

We compared the ability of two theoretical methods of pH-dependent conformational calculations to reproduce experimental potentiometric titration curves of two models of peptides: Ac-K5-NHMe in 95% methanol (MeOH)/5% water mixture and Ac-XX(A)7OO-NH2 (XAO) (where X is diaminobutyric acid, A is alanine, and O is ornithine) in water, methanol (MeOH), and dimethyl sulfoxide (DMSO), respectively. The titration curve of the former was taken from the literature, and the curve of the latter was determined in this work. The first theoretical method involves a conformational search using the electrostatically driven Monte Carlo (EDMC) method with a low-cost energy function (ECEPP/3 plus the SRFOPT surface-solvation model, assumming that all titratable groups are uncharged) and subsequent reevaluation of the free energy at a given pH with the Poisson-Boltzmann equation, considering variable protonation states. In the second procedure, molecular dynamics (MD) simulations are run with the AMBER force field and the generalized Born model of electrostatic solvation, and the protonation states are sampled during constant-pH MD runs. In all three solvents, the first pKa of XAO is strongly downshifted compared to the value for the reference compounds (ethylamine and propylamine, respectively); the water and methanol curves have one, and the DMSO curve has two jumps characteristic of remarkable differences in the dissociation constants of acidic groups. The predicted titration curves of Ac-K5-NHMe are in good agreement with the experimental ones; better agreement is achieved with the MD-based method. The titration curves of XAO in methanol and DMSO, calculated using the MD-based approach, trace the shape of the experimental curves, reproducing the pH jump, while those calculated with the EDMC-based approach and the titration curve in water calculated using the MD-based approach have smooth shapes characteristic of the titration of weak multifunctional acids with small differences between the dissociation constants. Nevertheless, quantitative agreement between theoretically predicted and experimental titration curves is not achieved in all three solvents even with the MD-based approach, which is manifested by a smaller pH range of the calculated titration curves with respect to the experimental curves. The poorer agreement obtained for water than for the nonaqueous solvents suggests a significant role of specific solvation in water, which cannot be accounted for by the mean-field solvation models.

A New Force Field (ECEPP-05) for Peptides, Proteins, and Organic Molecules

Parametrization and testing of a new all-atom force field for organic molecules and peptides with fixed bond lengths and bond angles are described. The van der Waals parameters for both the organic molecules and the peptides were taken from J. Phys. Chem. B 2003, 107, 7143 and J. Phys. Chem. B 2004, 108, 12181. First, the values of the 1-4 nonbonded and electrostatic scale factors appropriate to the new force field were determined by computing the conformational energies of six model molecules, namely, ethanol, ethylamine, propanol, propylamine, 1,2-ethanediol, and 1,3-propanediol with different values of these factors. The partial atomic charges of these molecules were obtained by fitting to the electrostatic potentials calculated with the HF/6-31G quantum-mechanical method. Two different charge models (single- and multiple-conformation-derived) were also considered. We demonstrated that the charge model has a stronger effect on the conformational energies than the 1-4 scaling. The choice of a charge model affected the conformational energies of even the smallest molecules considered, whereas the effect of the 1-4 electrostatic or nonbonded scaling was apparent only for 1,3-propanediol. The best agreement with high-level ab initio data was obtained with the multiple-conformation-derived charges and with no scaling of the 1-4 nonbonded or electrostatic interactions (scale factors of 1.0). Next, the torsional parameters of a large number of neutral and charged organic molecules, assumed to be models of the side chains of the 20 naturally occurring amino acids, were computed by fitting to rotational energy profiles obtained from ab initio MP2/6-31G calculations. The quality of the fits was high with average errors for torsional profiles of less than 0.2 kcal/mol. To derive the torsional parameters for the peptide backbone, the partial atomic charges of the 20 neutral and charged amino acids were obtained by fitting to the electrostatic potentials of terminally blocked amino acids using the HF/6-31G quantum-mechanical method. Then, the phi-psi energy maps of Ac-Ala-NMe and Ac-Gly-NMe were computed using MP2/6-31G//HF/6-31G quantum-mechanical methods. The phi-psi energy map of Ac-Ala-NMe was used for refinement of the nonbonded parameters for the backbone nitrogen and hydrogen bonded to it. Subsequently, the main-chain torsional parameters were obtained by fitting the molecular mechanics energies to the phi-psi energy maps of Ac-Ala-NMe and Ac-Gly-NMe. The transferability of the entire force field was demonstrated by reproducing the main energy minima of terminally blocked Ala3 from the literature. The performance of the force field was also evaluated by simulating crystal structures of small peptides. By comparison of simulated and experimental data, examination of the torsional-angle and atom-positional root-mean-square deviations of the energy-minimized crystal structures from the corresponding X-ray model structures demonstrated high accuracy of the force field.