Simulating electrostatic energies in proteins: perspectives and some recent studies of pKas, redox, and other crucial functional properties

Electrostatic energies provide what is arguably the most effective tool for structure-function correlation of biological molecules. Here, we provide an overview of the current state-of-the-art simulations of electrostatic energies in macromolecules, emphasizing the microscopic perspective but also relating it to macroscopic approaches. We comment on the convergence issue and other problems of the microscopic models and the ways of keeping the microscopic physics while moving to semi-macroscopic directions. We discuss the nature of the protein dielectric "constants" reiterating our long-standing point that the dielectric "constants" in semi-macroscopic models depend on the definition and the specific treatment. The advances and the challenges in the field are illustrated considering different functional properties including pK(a)'s, redox potentials, ion and proton channels, enzyme catalysis, ligand binding, and protein stability. We emphasize the microscopic overcharging approach for studying pK(a) 's of internal groups in proteins and give a demonstration of power of this approach. We also emphasize recent advances in coarse grained models with a physically based electrostatic treatment and provide some examples including further directions in treating voltage activated ion channels.     

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.     

Calculations of antibody-antigen interactions: microscopic and semi-microscopic evaluation of the free energies of binding of phosphorylcholine analogs to McPC603.

The study of antibody-antigen interactions should greatly benefit from the development of quantitative models for the evaluation of binding free energies in proteins. The present work addresses this challenge by considering the test case of the binding free energies of phosphorylcholine analogs to the murine myeloma protein McPC603. This includes the evaluation of the differential binding energy as well as the absolute binding energies and their corresponding electrostatic contributions. Four different approaches are examined: the Protein Dipoles Langevin Dipoles (PDLD) method, the semi-microscopic PDLD (PDLD/S) method, a free energy perturbation (FEP) method based on an adiabatic charging procedure and a linear response approximation that accelerates the FEP calculation. The PDLD electrostatic calculations are augmented by estimates of the relevant hydrophobic and steric contributions. The determination of the hydrophobic energy involves an approach which considers the modification of the effective surface area of the solute by local field effects. The steric contributions are analyzed in terms of the corresponding reorganization energies. This treatment, which considers the protein as a harmonic system, views the steric forces as the restoring forces for the electrostatic interactions. The FEP method is found to give unreliable results with regular cut-off radii and starts to give quantitative results only in very expensive treatment with very large cut-off radii. The PDLD and PDLD/S methods are much faster than the FEP approach and give reasonable results for both the relative and absolute binding energies. The speed and simplicity of the PDLD/S method make it an effective strategy for interactive docking studies and indeed such an option is incorporated in the program MOLARIS. A component analysis of the different energy contributions of the FEP treatment and a similar PDLD analysis indicate that electrostatic effects provide the largest contribution to the differential binding energy, while the hydrophobic and steric contributions are much smaller. This finding lends further support to the idea that electrostatic interactions play a major role in determining the antigen specificity of McPC603.     

Distinguishing native conformations of proteins from decoys with an effective free energy estimator based on the OPLS all-atom force field and the Surface Generalized Born solvent model

Protein decoy data sets provide a benchmark for testing scoring functions designed for fold recognition and protein homology modeling problems. It is commonly believed that statistical potentials based on reduced atomic models are better able to discriminate native-like from misfolded decoys than scoring functions based on more detailed molecular mechanics models. Recent benchmark tests on small data sets, however, suggest otherwise. In this work, we report the results of extensive decoy detection tests using an effective free energy function based on the OPLS all-atom (OPLS-AA) force field and the Surface Generalized Born (SGB) model for the solvent electrostatic effects. The OPLS-AA/SGB effective free energy is used as a scoring function to detect native protein folds among a total of 48,832 decoys for 32 different proteins from Park and Levitt's 4-state-reduced, Levitt's local-minima, Baker's ROSETTA all-atom, and Skolnick's decoy sets. Solvent electrostatic effects are included through the Surface Generalized Born (SGB) model. All structures are locally minimized without restraints. From an analysis of the individual energy components of the OPLS-AA/SGB energy function for the native and the best-ranked decoy, it is determined that a balance of the terms of the potential is responsible for the minimized energies that most successfully distinguish the native from the misfolded conformations. Different combinations of individual energy terms provide less discrimination than the total energy. The results are consistent with observations that all-atom molecular potentials coupled with intermediate level solvent dielectric models are competitive with knowledge-based potentials for decoy detection and protein modeling problems such as fold recognition and homology modeling.     

What are the dielectric “constants” of proteins and how to validate electrostatic models?

Implicit models for evaluation of electrostatic energies in proteins include dielectric constants that represent effect of the protein environment. Unfortunately, the results obtained by such models are very sensitive to the value used for the dielectric constant. Furthermore, the factors that determine the optimal value of these constants are far from being obvious. This review considers the meaning of the protein dielectric constants and the ways to determine their optimal values. It is pointed out that typical benchmarks for validation of electrostatic models cannot discriminate between consistent and inconsistent models. In particular, the observed pK(a) values of surface groups can be reproduced correctly by models with entirely incorrect physical features. Thus, we introduce a discriminative benchmark that only includes residues whose pK(a) values are shifted significantly from their values in water. We also use the semimacroscopic version of the protein dipole Langevin dipole (PDLD/S) formulation to generate a series of models that move gradually from microscopic to fully macroscopic models. These include the linear response version of the PDLD/S models, Poisson Boltzmann (PB)-type models, and Tanford Kirkwwod (TK)-type models. Using our different models and the discriminative benchmark, we show that the protein dielectric constant, epsilon(p), is not a universal constant but simply a parameter that depends on the model used. It is also shown in agreement with our previous works that epsilon(p) represents the factors that are not considered explicitly. The use of a discriminative benchmark appears to help not only in identifying nonphysical models but also in analyzing effects that are not reproduced in an accurate way by consistent models. These include the effect of water penetration and the effect of the protein reorganization. Finally, we show that the optimal dielectric constant for self-energies is not the optimal constant for charge-charge interactions.     

The barrier for proton transport in aquaporins as a challenge for electrostatic models: the role of protein relaxation in mutational calculations

The origin of the barrier for proton transport through the aquaporin channel is a problem of general interest. It is becoming increasingly clear that this barrier is not attributable to the orientation of the water molecules across the channel but rather to the electrostatic penalty for moving the proton charge to the center of the channel. However, the reason for the high electrostatic barrier is still rather controversial. It has been argued by some workers that the barrier is due to the so-called NPA motif and/or to the helix macrodipole or to other specific elements. However, our works indicated that the main reason for the high barrier is the loss of the generalized solvation upon moving the proton charge from the bulk to the center of the channel and that this does not reflect a specific repulsive electrostatic interaction but the absence of sufficient electrostatic stabilization. At this stage it seems that the elucidation and clarification of the origin of the electrostatic barrier can serve as an instructive test case for electrostatic models. Thus, we reexamine the free-energy surface for proton transport in aquaporins using the microscopic free-energy perturbation/umbrella sampling (FEP/US) and the empirical valence bond/umbrella sampling (EVB/US) methods as well as the semimacroscopic protein dipole Langevin dipole model in its linear response approximation version (the PDLD/S-LRA). These extensive studies help to clarify the nature of the barrier and to establish the "reduced solvation effect" as the primary source of this barrier. That is, it is found that the barrier is associated with the loss of the generalized solvation energy (which includes of course all electrostatic effects) upon moving the proton charge from the bulk solvent to the center of the channel. It is also demonstrated that the residues in the NPA region and the helix dipole cannot be considered as the main reasons for the electrostatic barrier. Furthermore, our microscopic and semimacroscopic studies clarify the problems with incomplete alternative calculations, illustrating that the effects of various electrostatic elements are drastically overestimated by macroscopic calculations that use a low dielectric constant and do not consider the protein reorganization. Similarly, it is pointed out that microscopic potential of mean force calculations that do not evaluate the electrostatic barrier relative to the bulk water cannot be used to establish the origin of the electrostatic barrier. The relationship between the present study and calculations of pK(a)s in protein interiors is clarified, pointing out that approaches that are applied to study the aquaporin barrier should be validated by pK(a)s calculations. Such calculations also help to clarify the crucial role of solvation energies in establishing the barrier in aquaporins.     

An evaluation of Poisson-Boltzmann electrostatic free energy calculations through comparison with experimental mutagenesis data

For systems involving highly and oppositely charged proteins, electrostatic forces dominate association and contribute to biomolecular complex stability. Using experimental or theoretical alanine-scanning mutagenesis, it is possible to elucidate the contribution of individual ionizable amino acids to protein association. We evaluated our electrostatic free energy calculations by comparing calculated and experimental data for alanine mutants of five protein complexes. We calculated Poisson-Boltzmann electrostatic free energies based on a thermodynamic cycle, which incorporates association in a reference (Coulombic) and solvated (solution) state, as well as solvation effects. We observe that Coulombic and solvation free energy values correlate with experimental data in highly and oppositely charged systems, but not in systems comprised of similarly charged proteins. We also observe that correlation between solution and experimental free energies is dependent on dielectric coefficient selection for the protein interior. Free energy correlations improve as protein dielectric coefficient increases, suggesting that the protein interior experiences moderate dielectric screening, despite being shielded from solvent. We propose that higher dielectric coefficients may be necessary to more accurately predict protein-protein association. Additionally, our data suggest that Coulombic potential calculations alone may be sufficient to predict relative binding of protein mutants.     

Mechanisms of tryptophan fluorescence shifts in proteins

Tryptophan fluorescence wavelength is widely used as a tool to monitor changes in proteins and to make inferences regarding local structure and dynamics. We have predicted the fluorescence wavelengths of 19 tryptophans in 16 proteins, starting with crystal structures and using a hybrid quantum mechanical-classical molecular dynamics method with the assumption that only electrostatic interactions of the tryptophan ring electron density with the surrounding protein and solvent affect the transition energy. With only one adjustable parameter, the scaling of the quantum mechanical atomic charges as seen by the protein/solvent environment, the mean absolute deviation between predicted and observed fluorescence maximum wavelength is 6 nm. The modeling of electrostatic interactions, including hydration, in proteins is vital to understanding function and structure, and this study helps to assess the effectiveness of current electrostatic models.     

Predicting the rate enhancement of protein complex formation from the electrostatic energy of interaction

The rate of association of proteins is dictated by diffusion, but can be enhanced by favorable electrostatic forces. Here the relationship between the electrostatic energy of interaction, and the kinetics of protein-complex formation was analyzed for the protein pairs of: hirudin-thrombin, acetylcholinesterase-fasciculin and barnase-barstar, and for a panel of point mutants of these proteins. Electrostatic energies of interaction were calculated as the difference between the electrostatic energy of the complex and the sum of the energies of the two individual proteins, using the computer simulation package DelPhi. Calculated electrostatic energies of interaction were compared to experimentally determined rates of association. One kcal/mol of Coulombic interaction energy increased the rate of association by a factor of 2.8, independent of the protein-complex or mutant analyzed. Electrostatic energies of interaction were also determined from the salt dependence of the association rate constant, using the same basic equation as for the theoretical calculation. A Br?nsted analysis of the electrostatic energies of interactions plotted versus experimentally determined ln(rate)s of association shows a linear relation between the two, with a beta value close to 1. This is interpreted as the energy of the transition state varies according to the electrostatic interaction energy, fitting a two state model for the association reaction. Calculating electrostatic rate enhancement from the electrostatic interaction energy can be used as a powerful tool to design protein complexes with altered rates of association and affinities.     

Electrostatic effects in unfolded staphylococcal nuclease

Structure-based calculations of pKa values and electrostatic free energies of proteins assume that electrostatic effects in the unfolded state are negligible. In light of experimental evidence showing that this assumption is invalid for many proteins, and with increasing awareness that the unfolded state is more structured and compact than previously thought, a detailed examination of electrostatic effects in unfolded proteins is warranted. Here we address this issue with structure-based calculations of electrostatic interactions in unfolded staphylococcal nuclease. The approach involves the generation of ensembles of structures representing the unfolded state, and calculation of Coulomb energies to Boltzmann weight the unfolded state ensembles. Four different structural models of the unfolded state were tested. Experimental proton binding data measured with a variant of nuclease that is unfolded under native conditions were used to establish the validity of the calculations. These calculations suggest that weak Coulomb interactions are an unavoidable property of unfolded proteins. At neutral pH, the interactions are too weak to organize the unfolded state; however, at extreme pH values, where the protein has a significant net charge, the combined action of a large number of weak repulsive interactions can lead to the expansion of the unfolded state. The calculated pKa values of ionizable groups in the unfolded state are similar but not identical to the values in small peptides in water. These studies suggest that the accuracy of structure-based calculations of electrostatic contributions to stability cannot be improved unless electrostatic effects in the unfolded state are calculated explicitly.     

Molecular Dynamics Simulations of a DMPC Bilayer Using Nonadditive Interaction Models

We present a polarizable force field based on the charge-equilibration formalism for molecular dynamics simulations of phospholipid bilayers. We discuss refinement of headgroup dihedral potential parameters to reproduce ab initio conformational energies of dimethylphosphate calculated at the MP2/cc-pVTZ level of theory. We also address the refinement of electrostatic and Lennard-Jones (van der Waals) parameters to reproduce ab initio polarizabilities and water interaction energies of dimethylphosphate and tetramethylammonium. We present results of molecular dynamics simulations of a solvated dimyristoylphosphatidylcholine bilayer using this polarizable force field as well as the nonpolarizable, fixed-charge CHARMM27 and CHARMM27r force fields for comparison. Calculated atomic and electron-density profiles, deuterium order parameters, and headgroup orientations are found to be consistent with previous simulations and with experiment. Polarizable interaction models for solvent and lipid exhibit greater water penetration into the lipid interior; this is due to the variation of water molecular dipole moment from a bulk value of 2.6 Debye to a value of 1.9 Debye in the membrane interior. The reduction in the electrostatic component of the desolvation free-energy penalty allows for greater water density. The surface dipole potential predicted by the polarizable model is 0.95 V compared to the value of 0.8 V based on nonpolarizable force-field calculations. Effects of inclusion of explicit polarization are discussed in relation to water dipole moment and varying charge distributions. Dielectric permittivity profiles for polarizable and nonpolarizable interactions exhibit subtle differences arising from the nature of the individual component parameterizations; for the polarizable force field, we obtain a bulk dielectric permittivity of 79, whereas the nonpolarizable force field plateaus at 97 (the value for pure TIP3P water). In the membrane interior, both models predict unit permittivities, with the polarizable models contributing from one to two more units due to the optical dielectric (high-frequency dipole fluctuations). This contribution is a step toward the continuing development of a CHARMM (Chemistry at Harvard Molecular Mechanics) polarizable force field for simulations of biomacromolecular systems.     

When cells lose water: Lessons from biophysics and molecular biology

Organisms living in deserts and anhydrobiotic species are useful models for unraveling mechanisms used to overcome water loss. In this context, late embryogenesis abundant (LEA) proteins and sugars have been extensively studied for protection against desiccation stress and desiccation tolerance. This article aims to reappraise the current understanding of these molecules by focusing on converging contributions from biochemistry, molecular biology, and the use of biophysical tools. Such tools have greatly advanced the field by uncovering intriguing aspects of protein 3-D structure, such as folding upon stress. We summarize the current research on cellular responses against water deficit at the molecular level, considering both plausible water loss-sensing mechanisms and genes governing signal transduction pathways. Finally, we propose models that could guide future experimentation, for example, by concentrating on the behavior of selected proteins in living cells.     

Multiple-site titration and molecular modeling: Two rapid methods for computing energies and forces for ionizable groups in proteins

Computer models of proteins frequently treat the energies and forces associated with ionizable groups as if they were purely electrostatic. This paper examines the validity of the purely electrostatic approach, and concludes that significant errors in energies can result from the neglect of ionization changes. However, a complete treatment of ionizable groups presents substantial computational obstacles, because of the large number of ionization states which must be examined in systems having multiple interacting titratable groups. In order to address this problem, two novel methods for treating the energetics and forces associated with ionizable groups with a minimum of computer time have been developed. The most rapid method yields approximate energies by computing the free energy of a single highly occupied ionization state. The second method separates ionizable groups into clusters, and treats intracluster interactions exactly, but intercluster interactions approximately. This method yields both accurate energies and fractional charges. Good results are obtained in tests of both methods on proteins having has many as 123 ionizable groups. The more rapid method requires computer times of 0.01 to 0.34 sec, while the more accurate method requires 0.7 to 15 sec. These methods may be fast enough to permit the incorporation of ionization effects in iterative computations, such as energy minimizations and conformational searches. © 1993 Wiley-Liss, Inc.     

Electrostatic potential of nucleotide-free protein is sufficient for discrimination between adenine and guanine-specific binding sites

Despite sharing many common features, adenine-binding and guanine-binding sites in proteins often show a clear preference for the cognate over the non-cognate ligand. We have analyzed electrostatic potential (ESP) patterns at adenine and guanine-binding sites of a large number of non-redundant proteins where each binding site was first annotated as adenine/guanine-specific or non-specific from a survey of primary literature. We show that more than 90% of ESP variance at the binding sites is accounted for by only two principal component ESP vectors, each aligned to molecular dipoles of adenine and guanine. Projected on these principal component vectors, the adenine/guanine-specific and non-specific binding sites, including adenine-containing dinucleotides, show non-overlapping distributions. Adenine or guanine specificities of the binding sites also show high correlation with the corresponding electrostatic replacement (cognate by non-cognate ligand) energies. High correlation coefficients (0.94 for 35 adenine-binding sites and 1.0 for 20 guanine-binding sites) were obtained when adenine/guanine specificities were predicted using the replacement energies. Our results demonstrate that ligand-free protein ESP is an excellent indicator for discrimination between adenine and guanine-specific binding sites and that ESP of ligand-free protein can be used as a tool to annotate known and putative purine-binding sites in proteins as adenine or guanine-specific.     

An overview of data-driven HADDOCK strategies in CAPRI rounds 38-45

Our information-driven docking approach HADDOCK has demonstrated a sustained performance since the start of its participation to CAPRI. This is due, in part, to its ability to integrate data into the modeling process, and to the robustness of its scoring function. We participated in CAPRI both as server and manual predictors. In CAPRI rounds 38-45, we have used various strategies depending on the available information. These ranged from imposing restraints to a few residues identified from literature as being important for the interaction, to binding pockets identified from homologous complexes or template-based refinement/CA-CA restraint-guided docking from identified templates. When relevant, symmetry restraints were used to limit the conformational sampling. We also tested for a large decamer target a new implementation of the MARTINI coarse-grained force field in HADDOCK. Overall, we obtained acceptable or better predictions for 13 and 11 server and manual submissions, respectively, out of the 22 interfaces. Our server performance (acceptable or higher-quality models when considering the top 10) was better (59%) than the manual (50%) one, in which we typically experiment with various combinations of protocols and data sources. Again, our simple scoring function based on a linear combination of intermolecular van der Waals and electrostatic energies and an empirical desolvation term demonstrated a good performance in the scoring experiment with a 63% success rate across all 22 interfaces. An analysis of model quality indicates that, while we are consistently performing well in generating acceptable models, there is room for improvement for generating/identifying higher quality models.