Calculations of electrostatic energies in proteins. The energetics of ionized groups in bovine pancreatic trypsin inhibitor

The importance of including different energy contributions in calculations of electrostatic energies in proteins is examined by calculating the intrinsic pKa values of the acidic groups of bovine pancreatic trypsin inhibitor. It appears that such calculations provide a powerful and revealing test; the relevant solvation energies of the ionized acids are of the order of -70 kcal/mol (1 cal = 4.184 J), and microscopic calculations that do not attempt to simulate the complete protein dielectric effect (including the surrounding solvent) can underestimate the solvation energy by as much as 50 kcal/mol. Reproducing correctly, by the same set of parameters, the solvation energies of ionized acids in different sites of a protein cannot be accomplished by including only part of the key energy contributions. The problems associated with macroscopic calculations are also considered and illustrated by the specific case of bovine pancreatic trypsin inhibitor. A promising approach is shown to be provided by a refinement of the previously developed Protein Dipoles Langevin Dipoles model. This model seems to represent consistently the microscopic dielectric of the protein and the surrounding water molecules. The model overcomes the problems associated with the macroscopic models (by treating explicitly the solvent molecules) and avoids the convergence problems associated with all-atom solvent models (by treating the average solvent polarization rather than averaging the actual polarization energy). This paper describes in detail the actual implementation of the model and examines its performance in evaluating intrinsic pKa values. Preliminary microscopic considerations of charge-charge interactions are presented.     

Evaluation of the conformational free energies of loops in proteins

In this paper we discuss the problem of including solvation free energies in evaluating the relative stabilities of loops in proteins. A conformational search based on a gas-phase potential function is used to generate a large number of trial conformations. As has been found previously, the energy minimization step in this process tends to pack charged and polar side chains against the protein surface, resulting in conformations which are unstable in the aqueous phase. Various solvation models can easily identify such structures. In order to provide a more severe test of solvation models, gas phase conformations were generated in which side chains were kept extended so as to maximize their interaction with the solvent. The free energies of these conformations were compared to that calculated for the crystal structure in three loops of the protein E. coli RNase H, with lengths of 7, 8, and 9 residues. Free energies were evaluated with a finite difference Poisson-Boltzmann (FDPB) calculation for electrostatics and a surface area-based term for nonpolar contributions. These were added to a gas-phase potential function. A free energy function based on atomic solvation parameters was also tested. Both functions were quite successful in selecting, based on a free energy criterion, conformations quite close to the crystal structure for two of the three loops. For one loop, which is involved in crystal contacts, conformations that are quite different from the crystal structure were also selected. A method to avoid precision problems associated with using the FDPB method to evaluate conformational free energies in proteins is described. © 1994 John Wiley & Sons, Inc.     

Theoretical studies of enzymic reactions: Dielectric,electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme

A general method for detailed study of enzymic reactions is presented. The method considers the complete enzyme-substrate complex together with the surrounding solvent and evaluates all the different quantum mechanical and classical energy factors that can affect the reaction pathway. These factors include the quantum mechanical energies associated with bond cleavage and charge redistribution of the substrate and the classical energies of steric and electrostatic interactions between the substrate and the enzyme. The electrostatic polarization of the enzyme atoms and the orientation of the dipoles of the surrounding water molecules is simulated by a microscopic dielectric model. The solvation energy resulting from this polarization is considerable and must be included in any realistic calculation of chemical reactions involving anything more than an isolated molecule in vacuo. Without it, acidic groups can never become ionized and the charge distribution on the substrate will not be reasonable. The same dielectric model can also be used to study the reaction of the substrate in solution. In this way the reaction in solution can be compared with the enzymic reaction.In this paper we study the stability of the carbonium ion intermediate formed in the cleavage of a glycosidic bond by lysozyme. It is found that electrostatic stabilization is an important factor in increasing the rate of the reaction step that leads to the formation of the carbonium ion intermediate. Steric factors, such as the strain of the substrate on binding to lysozyme, do not seem to contribute significantly.     

Evaluation of catalytic free energies in genetically modified proteins

A combination of the empirical valence bond method and a free energy perturbation approach is used to simulate the activity of genetically modified enzymes. The simulations reproduce in a semiquantitative way the observed effects of mutations on the activity and binding free energies of trypsin and subtilisin. This suggests that we are approaching a stage of quantitative structure-function correlation of enzymes. The analysis of the calculations points towards the electrostatic energy of the reacting system as the key factor in enzyme catalysis. The changes in the charges of the reacting system and the corresponding changes in "solvation" free energy (generalized here as the interaction between the charges and the given microenvironment) are emphasized. It is argued that a reliable evaluation of these changes might be sufficient for correlating structure and catalysis. The use of free energy perturbation methods and thermodynamic cycles for evaluation of solvation energies and reactivity is discussed, pointing out our early contributions. The apparent elaborated nature of our treatment is clarified, explaining that such a treatment is essential for consistent calculations of chemical reactions in polar environments. The problems associated with seemingly more rigorous quantum mechanical methods are discussed, emphasizing the inconsistency associated with using gas phase charge distributions. The importance of dynamic aspects is examined by evaluating the autocorrelation of the protein "reaction field" on the reacting substrate. It is found that, at least in the present case, dynamic effects are not important. The nature of the catalytic free energy is considered, arguing that the protein provides preoriented dipoles (polarized to stabilize the transition state charge distribution) and small reorganization energy, thus reducing the activation free energy. The corresponding catalytic free energy is related to the folding free energy, which is being invested in aligning the active site dipoles.     

A CFF 91-based Continuum Solvation Model: Solvation Free Energies of Small Organic Molecules and Conformations of the Alanine Dipeptide in Solution

Abstract

For molecular mechanics simulations of solvated molecules, it is important to use a consistent approach for calculating both the force field energy and the solvation free energy. A continuum solvation model based upon the atomic charges provided with the CFF91 force field is derived. The electrostatic component of the solvation free energy is described by the Poisson-Bolzmann equation while the nonpolar comonent of the solvation energy is assumed to be proportional to the solvent accessible surface area of the solute. Solute atomic radii used to describe the interface between the solute and solvent are fitted to reproduce the energies of small organic molecules. Data for 140 compounds are presented and compared to experiment and to the results from the well-characterized quantum mechanical solvation model AM1-SM2. In particular, accurate results are obtained for amino acid neutral analogues (mean unsigned error of 0.3 kcal/mol). The conformational energetics of the solvated alanine dipeptide is discussed.     

On the calculation of electrostatic interactions in proteins

In this paper we present a classical treatment of electrostatic interactions in proteins. The protein is treated as a region of low dielectric constant with spherical charges embedded within it, surrounded by an aqueous solvent of high dielectric constant, which may contain a simple electrolyte. The complete analysis includes the effects of solvent screening, polarization forces, and self energies, which are related to solvation energies. Formulae, and sample calculations of forces and energies, are given for the special case of a spherical protein. Our analysis and model calculations point out that any consistent treatment of electrostatic interactions in proteins should account for the following. Solvent polarization is an important factor in the calculation of pairwise electrostatic interactions. Solvent polarization substantially affects both electrostatic energies and forces acting upon charges. No simple expression for the effective dielectric constant, Deff, can generally be valid, since Deff is a sensitive function of position. Solvent screening of pairwise interactions involving dipolar groups is less effective than the screening of charges. In fact for many interactions involving dipoles, solvent screening can be essentially ignored. The self energy of charges makes a large contribution to the total electrostatic energy of a protein. This must be compensated by specific interactions with other groups in the protein. Strategies for applying our analysis to proteins whose structures are known are discussed.     

Modeling loop reorganization free energies of acetylcholinesterase: a comparison of explicit and implicit solvent models

The treatment of hydration effects in protein dynamics simulations varies in model complexity and spans the range from the computationally intensive microscopic evaluation to simple dielectric screening of charge-charge interactions. This paper compares different solvent models applied to the problem of estimating the free-energy difference between two loop conformations in acetylcholinesterase. Molecular dynamics (MD) simulations were used to sample potential energy surfaces of the two basins with solvent treated by means of explicit and implicit methods. Implicit solvent methods studied include the generalized Born (GB) model, atomic solvation potential (ASP), and the distance-dependent dieletric constant. By using the linear response approximation (LRA), the explicit solvent calculations determined a free-energy difference that is in excellent agreement with the experimental estimate, while rescoring the protein conformations with GB or the Poisson equation showed inconsistent and inferior results. While the approach of rescoring conformations from explicit water simulations with implicit solvent models is popular among many applications, it perturbs the energy landscape by changing the solvent contribution to microstates without conformational relaxation, thus leading to non-optimal solvation free energies. Calculations applying MD with a GB solvent model produced results of comparable accuracy as observed with LRA, yet the electrostatic free-energy terms were significantly different due to optimization on a potential energy surface favored by an implicit solvent reaction field. The simpler methods of ASP and the distance-dependent scaling of the dielectric constant both produced considerable distortions in the protein internal free-energy terms and are consequently unreliable.     

Conformational statistics of pectin substances in solution by a Metropolis Monte Carlo study

Using a Metropolis Monte Carlo algorithm, various average properties of several pectic polysaccharide models were calculated based on the conformational energies for parent disaccharides. The relaxed potential energy surfaces of disaccharides were calculated using the CHARMm force field. Solvent effects were evaluated by calculating a solvation energy for each conformational state by estimating contributions from a cavity formation, and from the electrostatic and dispersion interactions between solvent and solute molecules. The behavior of the mean characteristic ratio, the squared radius of gyration, and the persistence length versus chain length is discussed for various structural models, temperature, and solvent. It is found that the unrefined model of the alternating co-polymer polygalacto-galacturonic acid in vacuum is consistent with the experimentally measured dimension of pectin in salt excess. This model is used to generate computer images of the characteristic conformation of pectin chain.     

Comparison of end-point continuum-solvation methods for the calculation of protein-ligand binding free energies

We have compared the predictions of ligand‐binding affinities from several methods based on end‐point molecular dynamics simulations and continuum solvation, that is, methods related to MM/PBSA (molecular mechanics combined with Poisson–Boltzmann and surface area solvation). Two continuum‐solvation models were considered, viz., the Poisson–Boltzmann (PB) and generalised Born (GB) approaches. The nonelectrostatic energies were also obtained in two different ways, viz., either from the sum of the bonded, van der Waals, nonpolar solvation energies, and entropy terms (as in MM/PBSA), or from the scaled protein–ligand van der Waals interaction energy (as in the linear interaction energy approach, LIE). Three different approaches to calculate electrostatic energies were tested, viz., the sum of electrostatic interaction energies and polar solvation energies, obtained either from a single simulation of the complex or from three independent simulations of the complex, the free protein, and the free ligand, or the linear‐response approximation (LRA). Moreover, we investigated the effect of scaling the electrostatic interactions by an effective internal dielectric constant of the protein (?_int). All these methods were tested on the binding of seven biotin analogues to avidin and nine 3‐amidinobenzyl‐1H‐indole‐2‐carboxamide inhibitors to factor Xa. For avidin, the best results were obtained with a combination of the LIE nonelectrostatic energies with the MM+GB electrostatic energies from a single simulation, using ?_int = 4. For fXa, standard MM/GBSA, based on one simulation and using ?_int = 4–10 gave the best result. The optimum internal dielectric constant seems to be slightly higher with PB than with GB solvation. © Proteins 2012; © 2012 Wiley Periodicals, Inc.     

Ligand solvation in molecular docking

Solvation plays an important role in ligand‐protein association and has a strong impact on comparisons of binding energies for dissimilar molecules. When databases of such molecules are screened for complementarity to receptors of known structure, as often occurs in structure‐based inhibitor discovery, failure to consider ligand solvation often leads to putative ligands that are too highly charged or too large. To correct for the different charge states and sizes of the ligands, we calculated electrostatic and non‐polar solvation free energies for molecules in a widely used molecular database, the Available Chemicals Directory (ACD). A modified Born equation treatment was used to calculate the electrostatic component of ligand solvation. The non‐polar component of ligand solvation was calculated based on the surface area of the ligand and parameters derived from the hydration energies of apolar ligands. These solvation energies were subtracted from the ligand‐receptor interaction energies. We tested the usefulness of these corrections by screening the ACD for molecules that complemented three proteins of known structure, using a molecular docking program. Correcting for ligand solvation improved the rankings of known ligands and discriminated against molecules with inappropriate charge states and sizes. Proteins 1999;34:4–16. © 1999 Wiley‐Liss, Inc.     

A precise analytical method for calculating the electrostatic energy of macromolecules in aqueous solution

A new method for calculating the total electrostatic free energy of a macromolecule in solution is presented. It is applicable to molecules of arbitrary shape and size, including membranes or macromolecular assemblies with substrate molecules and ions. The method is derived from integrating the energy density of the electrostatic field and is termed the field energy method. It is based on the dielectric model, in which the solute and the surrounding water are regarded as different continuous dielectrics. The field energy method yields both the interaction energy between all charge pairs and the self energy of single charges, effectively accounting for the interaction with water. First, the dielectric boundary and mirror charges are determined for all charges of the solute. The energy is then given as a simple function of the interatomic distances, and the standard atomic partial charges and volumes. The interaction and self energy are shown to result from three-body and pairwise interactions. Both energy terms explicitly involve apolar atoms, revealing that apolar groups are also subject to electrostatic forces. We applied the field energy method to a spherical model protein. Comparison with the Kirkwood solution shows that errors are within a small percentage. As a further test, the field energy method was used to calculate the electrostatic potential of the protein superoxide dismutase. We obtained good agreement with the result from a program that implements the numerical finite difference algorithm. The field energy method provides a basis for energy minimization and dynamics programs that account for the solvent and screening effect of water at little computational expense.     

Molecular dynamics simulation of highly charged proteins: comparison of the particle-particle particle-mesh and reaction field methods for the calculation of electrostatic interactions

Molecular dynamics (MD) simulations of the activation domain of porcine procarboxypeptidase B (ADBp) were performed to examine the effect of using the particle-particle particle-mesh (P3M) or the reaction field (RF) method for calculating electrostatic interactions in simulations of highly charged proteins. Several structural, thermodynamic, and dynamic observables were derived from the MD trajectories, including estimated entropies and solvation free energies and essential dynamics (ED). The P3M method leads to slightly higher atomic positional fluctuations and deviations from the crystallographic structure, along with somewhat lower values of the total energy and solvation free energy. However, the ED analysis of the system leads to nearly identical results for both simulations. Because of the strong similarity between the results, both methods appear well suited for the simulation of highly charged globular proteins in explicit solvent. However, the lower computational demand of the RF method in the present implementation represents a clear advantage over the P3M method.     

Solvation effects in calculated electrostatic association free energies for the C3d‐CR2 complex and comparison with experimental data

The complement system is an integral part of the innate immune system that participates in the clearance of pathogens from the body. The association between complement protein fragment C3d and B or T cell‐receptor complement receptor (CR) 2 represents a crucial link between innate and adaptive immunities. The goal of this study is to predict association abilities of C3d and CR2 mutants by theoretically calculating electrostatic free energies of association and to assess the importance of solvation effects in the calculations. We demonstrate that calculated solvation free energy differences and Coulombic free energies of association are more sensitive than electrostatic free energies of association in solution and, thus, more accurate in predicting previously published experimental data for the association abilities (relative to the parent proteins) of specific C3d and CR2 mutants. We show that a proportional relationship exists between the predicted solvation free energy differences and the experimental data, while an inversely proportional relationship exists between the predicted Coulombic free energies of association and the experimental data. Our results yield new insights into the physicochemical properties underlying C3d‐CR2 association. We discuss the predictive validity of Coulombic, solvation, and solution electrostatic free energies of association and the generalization of our method for theoretical mutagenesis studies of other systems. This is a basic study, aimed toward improving our understanding of the theoretical basis of immune system regulation at the molecular level. Such insight can serve as the groundwork for the design of regulators with tailored properties, vaccines, and other biotechnology products. © 2010 Wiley Periodicals, Inc. Biopolymers 93:509–519, 2010. This article was originally published online as an accepted preprint. The “Published Online” date corresponds to the preprint version. You can request a copy of the preprint by emailing the Biopolymers editorial office at biopolymers@wiley.com     

Investigating protein-protein interaction surfaces using a reduced stereochemical and electrostatic model

A method of calculating the electrostatic potential energy between two molecules, using finite difference potential, is presented. A reduced charge set is used so that the interaction energy can be calculated as the two static molecules explore their full six-dimensional configurational space. The energies are contoured over surfaces fixed to each molecule with an interactive computer graphics program. For two crystal structures (trypsin-trypsin inhibitor and anti-lysozyme Fab-lysozyme), it is found that the complex corresponds to highly favourable interacting regions in the contour plots. These matches arise from a small number of protruding basic residues interacting with enhanced negative potential in each case. The redox pair cytochrome c peroxidase-cytochrome c exhibits an extensive favourably interacting surface within which a possible electron transfer complex may be defined by an increased electrostatic complementarity, but a decreased electrostatic energy. A possible substrate transfer configuration for the glycolytic enzyme pair glyceraldehyde phosphate dehydrogenase-phosphoglycerate kinase is presented.     

Free energy determinants of tertiary structure and the evaluation of protein models

We develop a protocol for estimating the free energy difference between different conformations of the same polypeptide chain. The conformational free energy evaluation combines the CHARMM force field with a continuum treatment of the solvent. In almost all cases studied, experimentally determined structures are predicted to be more stable than misfolded "decoys." This is due in part to the fact that the Coulomb energy of the native protein is consistently lower than that of the decoys. The solvation free energy generally favors the decoys, although the total electrostatic free energy (sum of Coulomb and solvation terms) favors the native structure. The behavior of the solvation free energy is somewhat counterintuitive and, surprisingly, is not correlated with differences in the burial of polar area between native structures and decoys. Rather. the effect is due to a more favorable charge distribution in the native protein, which, as is discussed, will tend to decrease its interaction with the solvent. Our results thus suggest, in keeping with a number of recent studies, that electrostatic interactions may play an important role in determining the native topology of a folded protein. On this basis, a simplified scoring function is derived that combines a Coulomb term with a hydrophobic contact term. This function performs as well as the more complete free energy evaluation in distinguishing the native structure from misfolded decoys. Its computational efficiency suggests that it can be used in protein structure prediction applications, and that it provides a physically well-defined alternative to statistically derived scoring functions.     

Continuum electrostatic analysis of preferred solvation sites around proteins in solution

To understand water-protein interactions in solution, the electrostatic field is calculated by solving the Poisson-Boltzmann equation, and the free energy surface of water is mapped by translating and rotating an explicit water molecule around the protein. The calculation is applied to T4 lysozyme with data available on the conservation of solvent binding sites in 18 crystallographically independent molecules. The free energy maps around the ordered water sites provide information on the relationship between water positions in crystal structure and in solution. Results show that almost all conserved sites and the majority of nonconserved sites are within 1.3 A of local free energy minima. This finding is in sharp contrast to the behavior of randomly placed water molecules in the boundary layer, which, on the average, must travel more than 3 A to the nearest free energy minimum. Thus, the solvation sites are at least partially determined by protein-water interactions rather than by crystal packing alone. The characteristic water residence times, obtained from the free energies at the local minima, are in good agreement with nuclear magnetic resonance experiments. Only about half of the potential sites show up as ordered water in the 1.7 A resolution X-ray structure. Crystal packing interactions can stabilize weak or mobile potential sites (in fact, some ordered water positions are not close to free energy minima) or can prevent water from occupying certain sites. Apart from a few buried water molecules that are strong binders, the free energies are not very different for conserved and nonconserved sites. We show that conservation of a water site between two crystals occurs if the positions of protein atoms, primarily contributing to the free energy at the local minimum, do not substantially change from one structure to the other. This requirement can be correlated with the nature of the side chain contacting the water molecule in the site.     

Computational method for relative binding energies of enzyme-substrate complexes.

A computational method for estimating the relative binding free energies of enzyme-substrate complexes is described that combines electrostatic and solvation models and X-ray crystallographic data. The polar contribution is evaluated by the Poisson-Boltzman equation. The nonpolar contribution is evaluated by solvent transfer data and surface area calculations. This algorithm was used to calculate the relative binding energies of 63 pairs of nine different mutant proteins with seven different substituted R-malate substrates of Escherichia coli isocitrate dehydrogenase. Comparison of calculated values with the experimentally observed values shows a high degree of correlation.     

Concepts and protocols for electrostatic free energies

Electrostatic free energies play an essential role in numerous biomolecular processes occurring in solution. Difficulties arise when the long-range Coulomb interaction is computed for idealised infinite simulation models with periodic boundary conditions. To maintain a neutral simulation box and a finite per-box energy, a neutralising charge density or ‘gellium’ is commonly used, leading to a mean box potential that is constrained to be rigorously equal to zero at all times. Thus, in considering quantities such as ion solvation free energy, the potential drop to move from solvent into the usual, gas phase reference state is missing. In fact, for an infinite molecular system, the electrostatic potential itself is not uniquely defined, but takes the form of an infinite series that is only conditionally convergent. This leads to several possible computational conventions that give different values for the potential and field, all mathematically valid. For experimentally measurable quantities, however, unique results are obtained when sufficiently large simulation boxes are utilised. These concepts are detailed, as well as a fundamental, linear response theoretical framework that provides qualitative understanding of the physical processes involved, especially dielectric relaxation of the environment in response to a new solute charge. Illustrative applications to ligand binding and biomolecular electron transfer are described.     

Excluded volume approximation to protein-solvent interaction. The solvent contact model.

Important properties of globular proteins, such as the stability of its folded state, depend sensitively on interactions with solvent molecules. Existing methods for estimating these interactions, such as the geometrical surface model, are either physically misleading or too time consuming to be applied routinely in energy calculations. As an alternative, we derive here a simple model for the interactions between protein atoms and solvent atoms in the first hydration layer, the solvent contact model, based on the conservation of the total number of atomic contacts, a consequence of the excluded-volume effect. The model has the conceptual advantage that protein-protein contacts and protein-solvent contacts are treated in the same language and the technical advantage that the solvent term becomes a particularly simple function of interatomic distances. The model allows rapid calculation of any physical property that depends only on the number and type of protein-solvent nearest-neighbor contacts. We propose use of the method in the calculation of protein solvation energies, conformational energy calculations, and molecular dynamics simulations.