Finite-difference Poisson-Boltzmann Research Articles

Electrostatic Contributions to Heat Capacity Changes of DNA-Ligand Binding

Significant heat capacity changes ( ΔC p) often accompany protein unfolding, protein binding, and specific DNA-ligand binding reactions. Such changes are widely used to analyze contributions arising from hydrophobic and polar hydration. Current models relate the magnitude of ΔC p to the solvent accessible surface area (ASA) of the molecule. However, for many binding systems—particularly those involving non-peptide ligands—these models predict a ΔC p that is significantly different from the experimentally measured value. Electrostatic interactions provide a potential source of heat capacity changes and do not scale with ASA. Using finite-difference Poisson-Boltzmann methods (FDPB), we have determined the contribution of electrostatics to the ΔC p associated with binding for DNA binding reactions involving the ligands DAPI, netropsin, lexitropsin, and the λ repressor binding domain.

The role of protonation and metal chelation preferences in defining the properties of mercury-binding coiled coils

To define the delicate interplay between metal chelation, protein folding and function in metalloproteins, a family of de novo-designed peptides was synthesized that self-assemble in aqueous solution to form two and three-stranded α-helical coiled coils. Each peptide contains a single Cys residue at an a or d position of the heptad repeat. Peptide association thus produces a Cys-rich coordination environment that has been used to bind Hg(II) ions. These peptides display a pH-dependent association, with trimers observed above the pKa of Glu side-chains and dimers below this value. Finite-difference Poisson-Boltzmann calculations suggest that the dimeric state decreases the unfavorable electrostatic interactions between positively charged Lys side-chains (relative to the trimer). The Cys-containing peptides bind Hg(II) in a position-dependent fashion. Cys at a positions form three-coordinate Hg complexes at high pH where the trimeric aggregation state predominates, and two-coordinate complexes at lower pH. A d position Cys, however, is only able to generate the two-coordinate complex, illustrating the difference in coordination geometry between the two positions in the coiled coil. The binding of Hg(II) was also shown to substantially increase the stability of the helical aggregates.

Investigating the dielectric effects of channel pore water on the electrostatic barriers of the permeation ion by the finite difference Poisson-Boltzmann method.

In this paper, the finite difference Poisson-Boltzmann (FDPB) method with four dielectric constants is developed to study the effect of dielectric saturation on the electrostatic barriers of the permeation ion. In this method, the inner shape of the channel pore is explicitly represented, and the fact that the dielectric constant inside the channel pore is different from that of bulk water is taken into account. A model channel system which is a righthanded twist bundle with four alpha-helical segments is provided for this study. From the FDPB calculations, it is found that the difference of the ionic electrostatic solvation energy for wider domains depends strongly on the pore radius in the vicinity of the ion when the pore dielectric constant is changed from 78 to 5. However, the electrostatic solvation energy of the permeation ion can not be significantly affected by the dielectric constant in regions with small pore radii. Our results indicate that the local electrostatic interactions inside the ion channel are of major importance for ion electrostatic solvation energies, and the effect of dielectric saturation on the electrostatic barriers is coupled to the interior channel dimensions.

Analysis of the stability of looped-out and stacked-in conformations of an adenine bulge in DNA using a continuum model for solvent and ions

A combination of conformational search, energy minimization, and energetic evaluation using a continuum solvent treatment has been employed to study the stability of various conformations of the DNA fragment d(CGCAGAA)/d(TTCGCG) containing a single adenine bulge. The extra-helical (looped-out) bulge conformation derived from a published x-ray structure and intra-helical (stacked bulge base) model structures partially based on nuclear magnetic resonance (NMR) data were used as start structures for the conformational search. Solvent-dependent contributions to the stability of the conformations were calculated from the solvent exposed molecular surface area and by using the finite difference Poisson-Boltzmann approach. Three classes (I-III) of bulge conformations with calculated low energies can be distinguished. The lowest-energy conformations were found in class I, corresponding to structures with the bulge base stacked between flanking helices, and class II, composed of structures forming a triplet of the bulge base and a flanking base pair. All extra-helical bulge structures, forming class III, were found to be less stable compared with the lowest energy structures of class I and II. The results are consistent with NMR data on an adenine bulge in the same sequence context indicating an intra-helical or triplet bulge conformation in solution. Although the total energies and total electrostatic energies of the low-energy conformations show only relatively modest variations, the energetic contributions to the stability were found to vary significantly among the classes of bulge structures. All intra-helical bulge structures are stabilized by a more favorable Coulomb charge-charge interaction but destabilized by a larger electrostatic reaction field contribution compared with all extra-helical and most triplet bulge structures. Van der Waals packing interactions and nonpolar surface-area-dependent contributions appear to favor triplet class II structures and to a lesser degree also the intra-helical stacked bulge conformations. The large conformational variation found for class III conformers might add a favorable entropic contribution to the stability of the extra-helical bulge form.

Does conformational free energy distinguish loop conformations in proteins?

Limitations in protein homology modeling often arise from the inability to adequately model loops. In this paper we focus on the selection of loop conformations. We present a complete computational treatment that allows the screening of loop conformations to identify those that best fit a molecular model. The stability of a loop in a protein is evaluated via computations of conformational free energies in solution, i.e., the free energy difference between the reference structure and the modeled one. A thermodynamic cycle is used for calculation of the conformational free energy, in which the total free energy of the reference state (i.e., gas phase) is the CHARMm potential energy. The electrostatic contribution of the solvation free energy is obtained from solving the finite-difference Poisson-Boltzmann equation. The nonpolar contribution is based on a surface area-based expression. We applied this computational scheme to a simple but well-characterized system, the antibody hypervariable loop (complementarity-determining region, CDR). Instead of creating loop conformations, we generated a database of loops extracted from high-resolution crystal structures of proteins, which display geometrical similarities with antibody CDRs. We inserted loops from our database into a framework of an antibody; then we calculated the conformational free energies of each loop. Results show that we successfully identified loops with a "reference-like" CDR geometry, with the lowest conformational free energy in gas phase only. Surprisingly, the solvation energy term plays a confusing role, sometimes discriminating "reference-like" CDR geometry and many times allowing "non-reference-like" conformations to have the lowest conformational free energies (for short loops). Most "reference-like" loop conformations are separated from others by a gap in the gas phase conformational free energy scale. Naturally, loops from antibody molecules are found to be the best models for long CDRs (> or = 6 residues), mainly because of a better packing of backbone atoms into the framework of the antibody model.

Molecular docking of superantigens with class II major histocompatibility complex proteins.

The molecular recognition of two superantigens with class II major histocompatibility complex molecules was simulated by using protein-protein docking. Superantigens studied were staphylococcal enterotoxin B (SEB) and toxic shock syndrome toxin-1 (TSST-1) in their crystallographic assemblies with HLA-DR1. Rigid-body docking was performed sampling configurational space of the interfacial surfaces by employing a strategy of partitioning the contact regions on HLA-DR1 into separate molecular recognition units. Scoring of docked conformations was based on an electrostatic continuum model evaluated with the finite-difference Poisson-Boltzmann method. Estimates of nonpolar contributions were derived from the buried molecular surface areas. We found for both superantigens that docking the HLA-DR1 surface complementary with the SEB and TSST-1 contact regions containing a homologous hydrophobic surface loop provided sufficient recognition for the reconstitution of native-like conformers exhibiting the highest-scoring free energies. For the SEB complex, the calculations were successful in reproducing the total association free energy. A comparison of the free-energy determinants of the conserved hydrophobic contact residue indicates functional similarity between the two proteins for this interface. Though both superantigens share a common global association mode, differences in binding topology distinguish the conformational specificities underlying recognition.

Monovalent and Divalent Salt Effects on Electrostatic Free Energies Defined by the Nonlinear Poisson−Boltzmann Equation: Application to DNA Binding Reactions

We have extended the finite-difference Poisson−Boltzmann (FDPB) equation method to incorporate the treatment of mixed salts (i.e. NaCl/MgCl2). In this context, we have derived an expression for the total electrostatic free energy for mixed salt systems. We use the theory to study nonspecific mixed salt effects on the binding free energies of the minor groove binding antibiotic DAPI, and λ repressor, with DNA. We find that in a pure salt solution the electrostatic contribution to binding varies linearly with log[Mn+], (where Mn+ represents an n-valent cation) and that the effect is uncorrelated with either the valence of the binding ligand or the number of counterions in the binding site. In mixed salt solution, the monovalent and divalent counterions “compete” for the immediate vicinity of the DNA. As experimentally observed in mixed salt solutions, a pronounced curvature appears in the plot of the electrostatic binding free energy vs log[Mn+]. The curvature for DAPI−DNA binding in mixed salts reflects th...

Simulation of electrostatic and hydrodynamic properties of Serratia endonuclease.

We analyze the electrostatic and hydrodynamic properties of a nuclease from the pathogenic gram-negative bacterium Serratia marcescens using finite-difference Poisson-Boltzmann methods for electrostatic calculations and a bead-model approach for diffusion coefficient calculations. Electrostatic properties are analyzed for the enzyme in monomeric and dimeric forms and also in the context of DNA binding by the nuclease. Our preliminary results show that binding of a double-stranded DNA dodecamer by nuclease causes an overall shift in the charge of the protein by approximately three units of elementary charge per monomer, resulting in a positively charged protein at physiologic pH. In these calculations, the free enzyme was found to have a negative (-1 e) charge per monomer at pH 7. The most dramatic shift in pKa involves His 89 whose pKa increases by three pH units upon DNA binding. This shift leads to a protonated residue at pH 7, in contrast to the unprotonated form in the free enzyme. DNA binding also leads to a decrease in the energetic distances between the most stable protonation states of the enzyme. Dimerization has no significant effect on the electrostatic properties of each of the monomers for both free enzyme and that bound to DNA. Results of hydrodynamic calculations are consistent with the dimeric form of the enzyme in solution. The computed translational diffusion coefficient for the dimer model of the enzyme is in very good agreement with measurements from light scattering experiments. Preliminary electrooptical calculations indicate that the dimer should possess a large dipole moment (approximately 600 Debye units) as well as substantial optical anisotropy (limiting reduced linear electric dichroism of about 0.3). Therefore, this system may serve as a good model for investigation of electric and hydrodynamic properties by relaxation electrooptical experiments.

Electrostatic interactions in hirudin-thrombin binding

Hirudin is a good anticoagulant owing to potent inhibition of the serine protease thrombin. An aspartate- and glutamate-rich portion of hirudin plays an important part in its tight binding to thrombin through a ladder of salt bridges, and these residues have previously been mutated to asparagine or glutamine. Detailed calculations of the electrostatic contribution to changes in binding from these mutations have been performed using the finite-difference Poisson-Boltzmann method which include charge-charge interactions, solvation interactions, the residual electrostatic interaction of mutant residues, p K a shifts, and ionic strength. Single mutant effects on binding energy were close to experimental values, except for the D55N mutant whose effect is overestimated, perhaps because of displacement of a bound chloride ion from the site where it binds. Multiple mutation values were generally overestimated. The effect of p K a shifts upon the binding is significant for one hirudin residue E58, but this appears to be due to a poor salt bridge with thrombin caused by crystal contacts. Electrostatic interaction between the acidic residues is unfavorable. However, analysis of experimental multiple mutation/single mutation data shows apparently negative interactions between these residues, from which it is concluded that structural changes can occur in the complex to relieve an unfavorable interaction when more than one acidic residue is mutated. In all cases, there is a loss in stability of the complex from mutations due to loss of favorable charge-charge interactions with thrombin, but this is largely compensated for by reduced unfavorable desolvation interactions, and by residual polar interactions in the Asn/Gln mutants.

Free Energy Determinants of Secondary Structure Formation: III. β-Turns and their Role in Protein Folding

The stability of β-turns is calculated as a function of sequence and turn type with a Monte Carlo sampling technique. The conformational energy of four internal hydrogen-bonded turn types, I, I′, II and II′, is obtained by evaluating their gas phase energy with the CHARMM force field and accounting for solvation effects with the Finite Difference Poisson- Boltzmann (FDPB) method. All four turn types are found to be less stable than the coil state, independent of the sequence in the turn. The free-energy penalties associated with turn formation vary between 1.6 kcal/mol and 7.7 kcal/mol, depending on the sequence and turn type. Differences in turn stability arise mainly from intraresidue interactions within the two central residues of the turn. For each combination of the two central residues, except for -Gly-Gly-, the most stable β-turn type is always found to occur most commonly in native proteins. The fact that a model based on local interactions accounts for the observed preference of specific sequences suggests that long-range tertiary interactions tend to play a secondary role in determining turn conformation. In contrast, for β-hairpins, long-range interactions appear to dominate. Specifically, due to the right-handed twist of β-strands, type I′ turns for -Gly-Gly- are found to occur with high frequency, even when local energetics would dictate otherwise. The fact that any combination of two residues is found able to adopt a relatively low-energy turn structure explains why the amino acid sequence in turns is highly variable. The calculated free-energy cost of turn formation, when combined with related numbers obtained for α-helices and β-sheets, suggests a model for the initiation of protein folding based on metastable fragments of secondary structure.

Calculations of the electrostatic free energy contributions to the binding free energy of sulfonamides to carbonic anhydrase

The interactions between biologically important enzymes and drugs are of great interest. In order to address some aspects of these interactions we have initiated a program to investigate enzymedrug interactions. Specifically, the interactions between one of the isozymes of carbonic anhydrase and a family of drugs known as sulfonamides have been studied using computational methods. In particular the electrostatic free energy of binding of carbonic anhydrase II with acetazolamide, methazolamide,p-chlorobenzenesulfonamide,p-aminobenzenesulfonamide and three new compounds (MK1, MK2, and MK3) has been computed using finite-difference Poisson-Boltzmann (FDPB) [1] method and the semimacroscopic version [2, 3] of the protein dipole Langevin dipole (PDLD) method [4]. Both methods, FDPB and PDLD, give similar results for the electrostatic free energy of binding even though different charges and different treatments were used for the protein. The calculated electrostatic binding free energies are in reasonable agreement with the experimental data. The potential and the limitation of electrostatic models for studies of binding energies are discussed.

Computation of the dipole moments of proteins

A simple and computationally feasible procedure for the calculation of net charges and dipole moments of proteins at arbitrary pH and salt conditions is described. The method is intended to provide data that may be compared to the results of transient electric dichroism experiments on protein solutions. The procedure consists of three major steps: (i) calculation of self energies and interaction energies for ionizable groups in the protein by using the finite-difference Poisson-Boltzmann method, (ii) determination of the position of the center of diffusion (to which the calculated dipole moment refers) and the extinction coefficient tensor for the protein, and (iii) generation of the equilibrium distribution of protonation states of the protein by a Monte Carlo procedure, from which mean and root-mean-square dipole moments and optical anisotropies are calculated. The procedure is applied to 12 proteins. It is shown that it gives hydrodynamic and electrical parameters for proteins in good agreement with experimental data.

Free Energy Determinants of Secondary Structure Formation: I. α-Helices

The Zimm-Bragg parameters sand σ are calculated for the helix-coil transition of poly- l-alanine. The theoretical approach involves evaluating gas phase conformational energies for both coil and helical states using the CHARMM potential function and accounting for solvation effects with various continuum solvation models. Conformational free energies are then incorporated into a formalism developed by Go et al.for the calculation of s and σ. Calculated values for both s and σ as well as the enthalpy change associated with helix formation are in good agreement with experimental data when the Finite Difference Poisson-Boltzmann (FDPB) method is used to treat solvent effects. The driving force for the helix-coil transition is analyzed in terms of individual free energy components. Hydrogen bond formation is found to contribute little to helix stability because the internal hydrogen bonding energy is largely canceled by the large free energy cost associated with removing polar groups from water. The entropic cost associated with fixing backbone dihedral angles in the helical conformation is found to be ∼7 e.u./residue (about 2 kcal/mol at room temperature). The major driving force favoring helix formation can be associated with interactions including enhanced van der Waals interactions in the close-packed helix conformation and the hydrophobic effect. These contribute about 2 kcal/mol favoring the helical state. The differences in helical propensities between alanine and glycine are attributed primarily to hydrophobic and packing interactions involving the C βwith a smaller contribution arising from increased conformational freedom for glycine in the coil state. The description of helix formation presented here is consistent with previous conclusions regarding tertiary structure formation which suggest that hydrophobic and close-packed interactions provide stability while hydrogen bond formation constitutes a structural constraint imposed by the high free energy cost associated with burying unsatisfied hydrogen bonding groups. α-Helix formation may thus be viewed as a form of hydrophobic collapse constrained by the requirement that polar groups be either exposed to solvent or form hydrogen bonds. More generally it appears from this study that for a folding model to be a realistic, it must properly account for the chemical nature of the polypeptide chain, particularly the solvation energetics of amide groups.

Electrostatics of hemoglobins from measurements of the electric dichroism and computer simulations

Hemoglobins from normal human cells, from sickle cells, and from horse were investigated by electrooptical methods in their oxy and deoxy forms. The reduced linear dichroism measured as a function of the electric field strength demonstrates the existence of permanent dipole moments in the range of 250-400 Debye units. The reduced limiting dichroism is relatively small (< or = 0.1); it is negative for hemoglobin from sickle cells and positive for the hemoglobins from normal human cells and from horse. The dichroism decay time constants are in the range from about 55 to 90 ns. Calculations of the electrooptical data from available crystal structures are given according to models of various complexity, including Monte Carlo simulations of proton fluctuations with energies evaluated by a finite difference Poisson-Boltzmann procedure. The experimental dipole moments are shown to be consistent with the results of the calculations. In the case of human deoxyhemoglobin, the root mean square dipole is higher than the mean dipole by a factor of about 4.5, indicating a particularly large relative contribution due to proton fluctuations. The ratio of the root mean square dipole to the mean dipole is much smaller (approximately 1.1 to approximately 1.5) for the other hemoglobin molecules. The calculations demonstrate that the dichroism decay time constants are not simply determined by the size/shape of the proteins, but are strongly influenced by the orientation of the dipole vector with respect to the axis of maximal absorbance. The comparison of experimental and calculated electrooptical data provides a useful test for the accuracy of electrostatic calculations and/or for the equivalence of structures in crystals and in solutions.

Intrinsic pKas of ionizable residues in proteins: an explicit solvent calculation for lysozyme.

Molecular dynamics simulations of triclinic hen egg white lysozyme in aqueous solution were performed to calculate the intrinsic pKas of 14 ionizable residues. An all-atom model was used for both solvent and solute, and a single 180 ps simulation in conjunction with a Gaussian fluctuation analysis method was used. An advantage of the Gaussian fluctuation method is that it only requires a single simulation of the system in a reference state to calculate all the pKas in the protein, in contrast to multiple simulations for the free energy perturbation method. pKint shifts with respect to reference titratable residues were evaluated and compared to results obtained using a finite difference Poisson-Boltzmann (FDPB) method with a continuum solvent model; overall agreement with the direction of the shifts was generally observed, though the magnitude of the shifts was typically larger with the explicit solvent model. The contribution of the first solvation shell to the total charging free energies of the titratable groups was explicitly evaluated and found to be significant. Dielectric shielding between pairs of titratable groups was examined and found to be smaller than expected. The effect of the approximations used to treat the long-range interactions on the pKint shifts is discussed.

Electrostatic potentials and electrostatic interaction energies of rat cytochrome b5 and a simulated anion-exchange adsorbent surface

Electrostatic potentials were determined for the soluble tryptic core of rat cytochrome b5 (using a structure derived from homology modeling) and a simulated anion-exchange surface through application of the linearized finite-difference Poisson-Boltzmann equation with the simulation code UHBD. Objectives of this work included determination of the contributions of the various charged groups on the protein surface to electrostatic interactions with a simulated anion-exchange surface as a function of orientation, separation distance, and ionic strength, as well as examining the potential existence of a preferred contact orientation. Electrostatic interaction free energies for the complex of the model protein and the simulated surface were computed using the electrostatics section of UHBD employing a 110(3) grid. An initial coarse grid spacing of 2.0 A was required to obtain correct boundary conditions. The boundary conditions of the coarse grid were used in subsequent focusing steps until the electrostatic interaction free energies were relatively independent of grid spacing (at approximately 0.5 A). Explicit error analyses were performed to determine the effects of grid spacing and other model assumptions on the electrostatic interaction free energies. The computational results reveal the presence of a preferred interaction orientation; the interaction energy between these two entities, of opposite net charge, is repulsive over a range of orientations. The electrostatic interaction free energies appear to be the summation of multiple fractional interactions between the protein and the anion-exchange surface. The simulation results are compared with those of ion-exchange adsorption experiments with site-directed mutants of the recombinant protein. Comparisons of the results from the computational and experimental studies should lead to a better understanding of electrostatic interactions of proteins and charged surfaces.

Combined Conformational Search and Finite-Difference Poisson?Boltzmann Approach for Flexible Docking: Application to an Operator Mutation in the λ Repressor-Operator Complex

The N-terminal domain of the phage λ repressor binds as a dimer to its palindromic DNA operator sequence. In addition to a helix-turn-helix DNA recognition motif, the first six amino acids of the phage λ repressor form a flexible peptide segment which wraps around DNA. Site-directed mutagenesis studies have shown that amino acid replacements or partial removal of the arm structure, or changes in the DNA sequence contacting the N-terminal arm, can lower the repressor-operator binding affinity by several orders of magnitude. The finite-difference Poisson-Boltzmann approach in combination with a conformational search procedure was used to study energetic contributions of the λ arm to repressor-operator recognition based on the high resolution X-ray structure. It allows for the local relaxation of the structure upon changing the DNA sequence in the λ arm binding region. A simplified potential energy function including torsional, truncated Lennard-Jones and approximate electrostatic terms is used in the initial step to screen out energetically unfavorable structures. The electrostatic energy of selected conformations is subsequently calculated more accurately using the finite-difference Poisson-Boltzmann approach. The method was applied to study the effect of a C→T mutation at position 6 of the consensus half-site of the operator. This base-pairs contacts Lys4 which is part of the arm segment. Keeping only the Lys4 side-chain mobile and with the wild-type DNA operator sequence, several conformations close to the X-ray structure were identified as those with lowest energy. In the case of the DNA mutation, lowest energy conformations differed significantly from those selected for the wild-type sequence. These initial calculations indicate that the approach might be a useful tool to estimate conformational and energetic effects upon mutagenesis of protein-DNA complexes.

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.

The role of electrostatic charge in the membrane insertion of colicin A. Calculation and mutation.

The bacterial toxin colicin A binds spontaneously to the surfaces of negatively charged membranes. The surface-bound toxin must subsequently, however, become an acidic 'molten globule' before it can fully insert into the lipid bilayer. Clearly, electrostatic interactions must play a significant role in both events. The electrostatic field around the toxin in solution was calculated using the finite-difference Poisson-Boltzmann method of the Delphi programme and the known X-ray structure. A large positively charged surface was identified which could be involved in the binding of colicin to negatively charged membranes. The applicability of the result was tested by also calculating the fields around modelled structures of the closely related colicins B and N. Surprisingly, colicin N showed a similar charge distribution in spite of its isoelectric point of pI 10.20 (colicin A has pI 5.44). One reason for this is the strong conservation of certain negative charges in all colicins. There is a single highly conserved aspartate residue (Asp78) on the positively charged face which provides a small but discrete region of negative charge. This residue, Asp78, was replaced by asparagine in the mutant D78N. D78N binds faster to negatively charged vesicles but inserts only half as fast as the wild-type protein into the membrane core. This indicates that, first, the initial membrane binding has a significant electrostatic component and, second, that the isolated charge on Asp78 plays a role in the formation of the insertion intermediate.

The inclusion of electrostatic hydration energies in molecular mechanics calculations

The problem of including solvent effects in molecular mechanics calculations is discussed. It is argued that the neglect of charge-solvent (solvation) interactions can introduce significant errors. The finite difference Poisson-Boltzmann (FDPB) method for calculating electrostatic interactions is summarized and is used as a basis for introducing a new pairwise energy term which accounts for charge-solvent interactions. This term acts between all pairs of atoms usually considered in molecular mechanics calculations and can be easily incorporated into existing force fields. As an example, a parameterization is developed for the CHARMm force field and the results compared to the predictions of the FDPB method. An approach to the realistic incorporation of solvent screening into force fields is also outlined.