INTRODUCTION

Top Abstract Introduction Theory Materials and methods Results and discussion Conclusions Appendix A Appendix B References

Electrostatics plays a key role in biological processes (Honig and Nicholls, 1995; Davis and McCammon, 1990; Davis et al., 1991). The binding of small electrolytes to a biomolecule in solution is kinetically driven by the electrostatic field generated by the molecule and is highly correlated with the electrostatic potential at the surface of the molecule. In many cases the nonobvious dependence of the kinetic constants of association between an enzyme and a substrate on the solution ionic conditions or kinetic pathways could be elucidated by analysis of the electrostatic fields in solution (Gilson et al., 1994; Sharp et al., 1987). Inspection of many molecular complexes has shown a high degree of complementarity in the electrostatic properties of the contacting surfaces (Honig and Nicholls, 1995). The electrostatic properties of biomolecular systems are influenced by pH and ionic conditions. The extent to which a group is ionized depends on the electrostatic potential generated at that site by the molecule (e.g., Antosiewicz et al., 1994). The ionization state of a biomolecule is in turn crucial for its function and stability.

The methods that have been used to simulate electrostatics in biological systems may be broadly classified into those which simulate explicitly all molecules of the system, including salts and solvent, which are by far the more demanding, and those which simulate the solvent and salts through a continuum model. Among the latter, the Poisson-Boltzmann equation (PBE) has been widely and successfully used. In recent years refined theoretical and numerical tools have been developed to apply the PBE to biomolecular systems (Gilson et al., 1987; Sharp and Honig, 1990; Zhou, 1994; Madura et al., 1995) and a large number of results have been achieved (Madura et al., 1994; Honig and Nicholls, 1995).

The reliability of the PBE has been tested for a few models and real systems by means of more sophisticated methods, such as Monte Carlo or hypernetted chain simulations (Fixman, 1979; Murthy et al., 1985, Jayaram and Beveridge, 1996).

The Poisson-Boltzmann equation was first put forward more than 80 years ago by Gouy (1910) and few years later by Chapman (1913). The equation was obtained either by equating to zero the forces acting on a microscopic volume of the ionic solution (Gouy, 1910) or by equating the chemical potential throughout the solution (Chapman, 1913). The same approaches have been followed by other researchers in the field of colloid chemistry (Derjaguin and Landau, 1941; Verwey and Overbeek, 1948) and electrocapillarity (Grahame, 1947).

Except for the simple planar geometry in the presence of symmetrical electrolytes (Gouy, 1910; Chapman, 1913) and the cylindrical geometry in the absence of added salt (Alfrey et al., 1951; Lifson and Katchalsky, 1954; Katchalsky, 1971), no analytical solution is available. Debye and Hückel (1923), who developed the PBE aiming at explicit calculation of the free energy for an ionic system, noticed that under usual experimental conditions the equation can be linearized to a good degree of accuracy for the computation of various thermodynamic quantities.

Although a number of software packages allow for the solution of the nonlinear PBE (Gilson et al., 1987; Madura et al., 1995), it is often mandatory to employ the linear approximation to reduce computation time. Depending on the system, the solution of the nonlinear PBE takes usually more than twice the time needed to solve its linear version. Moreover, since the electrostatic potential in the linear PBE (LPBE) is the superposition of the electrostatic potentials of each partial charge on the molecule, for all those applications where a charge is modified without altering the molecular shape (like in idealized protonation or deprotonation of an ionizable group), additional computing time is saved (Antosiewicz et al., 1994).

There are at least four major applications of the PBE and its linear form:

	(1)  calculation of the electrostatic potential at the surface of a biomolecule, which is expected to give information about the concentration of small charged solutes in the neighborhood of the molecule and whose inspection may suggest docking sites for biomolecules;
	(2)  calculation of the electrostatic potential outside the molecule, which is expected to give information on the free energy of interaction of small molecules at different positions in the surrounding of the molecule. The electrostatic field is therefore used in Brownian dynamics simulations employing the so-called test charge approximation;
	(3)  calculation of the free energy of a biomolecule or of different states of a biomolecule which gives information on the stability of a biomolecule or of its different states (Sharp and Honig, 1990); and
	(4)  calculation of the electrostatic field to derive mean forces to be added in standard molecular dynamics calculations (Gilson et al., 1993).

It is of interest, therefore, to investigate the limits of applicability of the LPBE for biomolecular systems and for these applications.

In the present study we address some of these issues, in particular:

	(1)  How accurate are the potentials derived via the LPBE for typical biomolecular systems?
	(2)  Is it possible to correct the biomolecular potential maps obtained via the LPBE in order to reproduce more faithfully the PBE results?
	(3)  How accurate is the free energy computed in the linear approximation?
	(4)  Is it possible to employ the LPBE potential to reach a better approximation of the PBE free energy?

We first compare the results obtained from the LPBE and PBE for systems with a simple geometry (i.e. the plane, the cylinder, and the sphere). Because the PBE for these shapes is characterized by a parameter (Gueron and Weisbuch, 1980) (m = 0, 1, 2 for the plane, the cylinder, and the sphere, respectively) we can heuristically set this parameter to intermediate values which could represent behaviors in intermediate cases.

Then we examine some biological systems and see how well the considerations for the simple shapes translate to these highly asymmetrical systems.

	THEORY

Top Abstract Introduction Theory Materials and methods Results and discussion Conclusions Appendix A Appendix B References

The PBE

In the Poisson-Boltzmann approach the macromolecule is treated as a low dielectric cavity with embedded atomic partial charges. The dielectric constant of the cavity is typically set between 2 and 4 to take into account electronic polarization and the limited flexibility of the macromolecule (Sharp et al., 1992; Gilson and Honig, 1986). The effects of the solvent molecules, whose motions are much faster than those of the molecule and the ions, are taken into account on average through a continuum of high dielectric constant (McCammon and Harvey, 1987).

The average electrostatic potential () is determined by the charge density embedded in the molecule (^f) and by the average charge density due to the mobile ions ^m, via the Poisson equation:

(1)

where  is the position-dependent dielectric constant and all terms are expressed in centimeter gram second-electrostatic units. The charge density ^m can be expressed in terms of the bulk concentrations and a potential of mean force:

(2)

where c_i is the concentration of ion i at an infinite distance from the molecule (or at any reference position where the potential of mean force w_i is set to zero), z_i is its charge number, q is the proton charge, k is the Boltzmann constant and T is the temperature.

The key assumptions to obtain the PBE are that the potentials of mean force are given by w_i = z_iqU and that U is equal to the average electrostatic potential :

(3)

When the term (z_iqU/kT) 1 the exponential can be expanded in a Taylor series, retaining only the first two terms. Due to electroneutrality, _ic_iz_iq = 0, the LPBE is obtained:

(4)

The most serious inconsistency of the PBE (Eq. 3) stems from the lack of reciprocity, i.e., different distributions are obtained for an ion pair by switching the definition of the central ion (Onsager, 1933; Fowler and Guggenheim, 1939). For some time this was regarded as an issue in favor of linearization.

Electrostatic free energy from the PBE

The electrostatic free energy for the hypothetical process of charging a sphere, organizing and charging the ionic atmosphere was earlier calculated according to the adiabate principle (Onsager, 1933; Verwey and Overbeek, 1948) where the free energy is obtained from the charging integral:

(5)

where q is the final charge on the sphere.

Another expression for the free energy of the process of charging the system, put forward by Marcus (1955), employs standard expressions for the chemical potential of solute molecules and is closely related to the expression we give below.

Sharp and Honig (1990) and, independently, Reiner and Radke (1990) derived the electrostatic free energy from a variational principle. They considered the PBE and built the Euler-Lagrange functional, which is extremized by the solution of the PBE. With an appropriate choice of multiplicative and additive constants, this functional could easily be interpreted as the free energy of the system.

The expression for the free energy is

(6)

though other forms, not involving derivatives of the potential, may be derived by exploiting the basic relationships _V((U)²/8)dV = _V(U/2)dV and c_i = c_i exp(z_iqU/kT) (Sharp and Honig, 1990).

The derivation faces several problems, however, including the paradoxical observation that the functional is not minimized but maximized. Nevertheless, it is possible to show that a proper free energy functional, defined by combining standard thermodynamics and the usual Poisson-Boltzmann approximations, is minimized by the ionic distribution obtained via the PBE (Fogolari and Briggs, 1997). Zhou (1994) showed that the free energy given by Eq. 6 may be alternatively obtained by a standard charging process (Eq. 5), and that the free energy is independent of the charging pathway.

For practical reasons we may rewrite the electrostatic free energy in terms of different contributions due to the electrostatic energy obtained by integrating U/2 over two regions entailing the fixed (G^ef) and mobile charges (G^em), and the entropic (for a discussion of the entropy in electrostatic systems see Sharp, 1995) free energy of mixing of mobile species (G^mob) and solvent (G^solv),

(7)

where the different contributions read:

(8)

(9)

(10)

(11)

The latter three terms may be further grouped into a single term to indicate the outer space contribution to the free energy density integral:

(12)

This decomposition of the free energy does not correspond to any thermodynamic pathway but, in fact, it is closely related to the way software packages compute the electrostatic free energy. Misra et al. (1994) considered a thermodynamic pathway for charging the molecule and organizing and charging the ionic atmosphere that allows identification of the non-salt-dependent contribution to the free energy of the system (G^ns), the contribution arising from the ionic atmosphere interaction with the molecule (G^im), the contribution from the ion-ion interactions (Gⁱⁱ), and the contribution from the entropy cost of organizing the ionic atmosphere around the solute (G^org).

The relationship of such a decomposition with the one given above (Eq. 7) is straightforward and is reported in Fogolari et al. (1997).

In the LPBE approach the only term contributing electrostatic free energy is G^ef (Sharp and Honig, 1990) up to the order of the linear approximation, though some simple corrections may be devised, as we discuss below.

Applications of the LPBE to biomolecular systems

It is generally recognized that when (qU/kT) 1 the PBE can be approximated by the LPBE which results from the approximation sinh(qU/kT)  qU/kT. But it is common experience, at least in biomolecular simulations, that the solution of the LPBE is close to the solution of the PBE even when qU/kT at the molecular surface is in the range of 1 to 2, although in such cases the hyperbolic sine is 20% to 80% larger than the corresponding linear approximation. For higher potentials, even when the potential is several kT/q, the solutions of the LPBE and the PBE are not as dramatically distant as sinh(qU/kT) and qU/kT are.

The LPBE solution is usually larger than the PBE one. For centrosymmetrical ions in symmetrical solutions Gronwall, La Mer, and Sandved (1928) have given a series correction to the solution of the LPBE, but such rather involved expansion is of little use when dealing with irregularly shaped molecules possessing uneven charge distributions.

Before approaching complex biomolecular systems we consider systems with a simple geometry, which can be highly idealized models for proteins, nucleic acids, and membranes. For these systems we find a general correction rule that brings the LPBE potential close to the PBE potential at the surface. We also define some simple rules to derive free energies from the solution of the LPBE, which include contributions to the free energy integral from the outer volume of the molecule.

Systems with a simple geometry

The PBE and LPBE have been numerically solved and compared for systems with a simple geometry (SSG) where the corresponding equations are dependent on either one or two variables, depending on the symmetry of the system. Much attention has been given to planar, cylindrical, and spherical shapes (Gueron and Weisbuch, 1980; Stigter 1978) and, more recently, to spheroidal geometries (Yoon and Kim, 1989; Hsu and Liu, 1996a,b).

Usually the equations are solved for simple boundary conditions, like constant surface charge or potential, or for a mix of these. This is an excellent approximation in the fields of colloid chemistry, where the surface charge is often controlled via ionizable groups sensitive to changes in pH, or electrocapillarity, where the electrode potential is externally controlled. However, it is bound to give only a very rough picture of biomolecules.

Moreover, SSG are very rough representations of real biomolecules. For instance the cylindrical model does an excellent job for regular biopolymers like DNA, but it is very difficult to model proteins with spheres or ellipsoids of constant charge. A more sophisticated approach was proposed by Kirkwood (1934), but still it appears too simplistic to represent real biomolecules. Nevertheless, SSG may be easily and extensively studied and conclusions reached about these systems may apply to complex systems. For these reasons SSG have received much attention in the past as model systems.

The relevant equations and definitions for SSG are reported in Appendix A. It is apparent that the solution of the PBE and all the derived thermodynamic quantities depend on the boundary conditions which may be imposed through the reduced electric field '(x₀) at the surface position expressed in Debye lengths. These are in turn determined by the interplay of three relevant length scales: the radius of curvature, the Debye length, and the electric field scale length, defined in Appendix B. Previous results obtained on SSG, summarized in Appendix B, showed similarities between the behavior of the PBE solution for systems with different geometry and showed that for all systems the ratio _D/ appears critical for the applicability of the linearization. Rather than studying the solution of the PBE, which depends on the shape and on the radius of curvature, we reasoned that the relationship between the solutions of the LPBE and the PBE should depend, in addition to the absolute values they can take, on the parameter _D/ itself. In particular they should be coincident when (_D/) 1, whereas for (_D/) 1 we have _PBE  2 ln(|_PBE/x|_x0|). Therefore we searched for a correction to be applied to the solution of the LPBE which depends only on the ratio _D/, to recover the PBE solution. For its simple connection with the boundary conditions we rewrite _D/ in reduced units: (_D/) = (2/'(x₀)), where the derivative is taken with respect to r/_D.

	MATERIALS AND METHODS

Top Abstract Introduction Theory Materials and methods Results and discussion Conclusions Appendix A Appendix B References

Calculation protocols

For SSG the one-variable PBE and LPBE were solved numerically using an adaptive Runge-Kutta fourth-order algorithm (Press et al., 1990). Tentative values were put forward for the value of the potential at the surface and the behavior of the potential or its derivative at 5 Debye lengths from the surface was checked. The guess value for the surface potential was reset until the reduced potential and its derivative were <0.005 at 5 Debye lengths. All thermodynamic quantities were then obtained using the discretized analogs of the equations reported in the theory section.

All biomolecular simulations were performed with the software package UHBD (Madura et al., 1995) using standard procedures. The calculations employed a grid of 110 × 110 × 110 points with a grid mesh of 1.37 Å and one focusing step for a final grid mesh of 0.51 Å. In all calculations the dielectric constants of the solvent and solute molecules were 78 and 4, respectively. The radius of the ions was 2.0 Å and the solvent probe radius was 1.4 Å.

For the test of the electrostatic potential inside the molecule we used a grid of 110 × 110 × 110 points with a grid mesh of 1.0 Å in order to have all surface points inside the grid.

We have run a few tests on different conformers of amino acids in model dipeptide and tripeptide compounds studied by Fogolari et al. (1998) and on some anthracycline drugs studied by Baginski et al. (1997). In all these cases, studied at 150 mM ionic strength, the LPBE and the PBE gave virtually identical results.

We have chosen the following systems as test cases: a complex between the Antennapedia homeodomain with Cys 39 substituted by a serine (Antp C39S HD) (Billeter et al., 1993) and a stretch of 31 base pairs of DNA as a highly charged system with positive and negative regions of irregular shape (for details on the construction of the molecular model see Fogolari et al. (1997)), the isolated homeodomain which possesses an extended arm with positively charged residues as a highly positively charged mainly globular but irregularly shaped system, the isolated DNA as a highly charged cylindrical system and monomeric bovine -lactoglobulin at pH 2, as a highly charged overall globular system. For the last system the most probable protonation state was obtained following the protocol of Antosiewicz et al. (1994) applied on the structure of the monomeric unit A, recently obtained by Sawyer and coworkers (Brownlow et al., 1996), but using the partial charges taken from the forcefield CHARMM (version 22) (Brooks et al., 1983). For this protein the presence of a stable core in the monomer with most native connectivities at pH 2 was established by Ragona et al. (1997) via NMR spectroscopy. Because in the most probable protonation state only few carboxylic groups are still deprotonated making the overall net charge positive and very high, we have decided to keep all the ionizable sites protonated, since in the present context this theoretical model is chosen only for the purpose of comparing the LPBE and PBE solutions.

Optimized parameters for liquid simulation charges and atomic radii (Jorgensen and Tirado-Rives, 1988, Pranata et al., 1991) were employed in the calculations on the homeodomain-DNA complex, and isolated DNA and homeodomain, while for -lactoglobulin the set of CHARMM charges and radii was used (Brooks et al., 1983). The temperature was set to 300 K. The net charge of the molecules is 47 for the homeodomain-DNA complex, 62 for the DNA, 15 for the homeodomain and 21 for -lactoglobulin. -lactoglobulin (2580 atoms) and homeodomain (790 atoms) have a radius of ~25 Å, while DNA (2754 atoms) is approximately a cylinder with radius 10 Å and length 100 Å. Thermodynamic quantities were computed from the output accessibility and potential maps. Surface points were obtained as the interfacial points in the solvent.

Computation times

Typically, 3800 to 6200 s were required by UHBD on a Silicon Graphics, Inc. (Mountain View, CA) O2, R5000, 180 MHz computer with 128 Mbytes RAM to solve the larger and focused grid. Corresponding times for the LPBE ranged from 1800-2600 s. Generating the corrected potential grid map and extracting thermodynamic quantities from the map takes <120 s, so that the correction procedure is negligible on the overall computation time. The generation of the potential inside the molecule, tested only for -lactoglobulin (2580 atoms and 6328 interfacial points) takes ~200 s, but this time could be greatly optimized by properly selecting the interfacial points and possibly by choosing faster ways to solve Laplace's equation with Dirichlet boundary conditions.

	RESULTS AND DISCUSSION

Top Abstract Introduction Theory Materials and methods Results and discussion Conclusions Appendix A Appendix B References

SSG

We have solved numerically the PBE and LPBE for a large number of boundary conditions and for different values of m (0, 0.4, 0.5, 0.8, 1.0, 1.2, 1.5, 1.6, and 2.0). Although noninteger values of m do not have a general physical counterpart, we expect these to represent intermediate cases between the three limiting simple shapes.

Surface electrostatic potentials

(13)

(14)

View larger version (21K): [in this window] [in a new window]   FIGURE 1   o = |x0 obtained from the PBE versus o = |x0 obtained from LPBE for various values of the shape parameters m and x0. The data computed for each value of (/x)|x0 and different values of m and x0 group along curves which have been described by Eq. 14: where PBE indicates the estimate for the correct (PBE) potential and the function A((/x)|x0) has the following form: A(u) = 3.037 + 0.1940u + 0.00227u2. The coefficients in the latter equation have been determined by best fit of all points in the plot.

Electrostatic potentials outside the model molecules

View larger version (15K): [in this window] [in a new window]   FIGURE 2   and /x plotted as a function of the reduced distance from the center of a charged cylinder (m = 1.0, x0 = 0.5, and (/x)|x0 = 20.0) for the PBE (solid lines) and LPBE (dotted lines). The LPBE potential, corrected according to the local LPBE electric field, is shown (dashed line). The plot of the "effective" LPBE potential and of its derivative is also shown (long dashed lines).

Free energies

View larger version (14K): [in this window] [in a new window]   FIGURE 3   Gel/kT from the PBE versus Gef/kT from the LPBE.

View larger version (18K): [in this window] [in a new window]   FIGURE 4   Gef/kT (circles) and Gout/kT (squares) contribution to the electrostatic free energy versus the electrostatic free energy Gel/kT from the PBE.

(15)

View larger version (13K): [in this window] [in a new window]   FIGURE 5   Gout/kT from the LPBE using the correcting formula (Eq. 14) (squares) or estimated through an effective potential (Eq. 15) from the LPBE (circles) versus Gout/kT obtained from the PBE. Only values in the range 0.0 to 3.0 for both variables are shown.

Biomolecular systems

Next we discuss the limits of applicability of the LPBE for a few systems of biological interest that may be representative of a diversity of shapes and that correspond to high charge densities. For low-charge systems the LPBE is usually in striking agreement with the PBE.

As for the SSG we discuss separately three areas where we can extend and put to use the results obtained for the SSG models: first, calculation of the electrostatic potential at the surface of the molecule; second, calculation of the electrostatic potential outside the molecule; and third, calculation of the free energies of the system. In addition we discuss the correction of the electrostatic potential inside the molecule.

Most of the results are illustrated for the homeodomain-DNA complex (Fig. 6), the most irregular of the systems studied in terms of shape and charge density. Results for the other systems are summarized in few figures and in one table. Physiological (145 mM) and low (14.5 mM) ionic strengths have been considered. The agreement between the LPBE and PBE increases with increasing ionic strength (i.e., with better screening of the potential).

View larger version (190K): [in this window] [in a new window]   FIGURE 6   The electrostatic potential (in kcal/q · mol units) at the surface of the Antp C39S HD-DNA complex at 14.5 mM ionic strength as obtained from the PBE, visualized with the software GRASP (Nicholls, 1993). A similar, but reduced in magnitude, potential pattern is obtained at 145 mM ionic strength. Only the potential from the last focused region is shown.

Surface electrostatic potentials

View larger version (59K): [in this window] [in a new window]   FIGURE 7   The electrostatic field (in reduced (kT/q)kD units) from the PBE (upper panel) and from the LPBE (lower panel) in a plane orthogonal to the DNA axis going through the center of geometry of the Antp C39S HD-DNA complex.

View larger version (18K): [in this window] [in a new window]   FIGURE 8   The average electrostatic potential (in kT/q units) from the LPBE (dashed lines) and from the corrected (Eq. 14) LPBE (solid lines) at the surface of the Antp C39S HD-DNA complex, the isolated DNA and HD and fully protonated -lactoglobulin at 14.5 mM ionic strength versus the electrostatic potential obtained from the PBE. CPLX, HD, DNA and BLG stand for the homeodomain-DNA complex, the isolated homeodomain, the isolated DNA stretch and fully protonated -lactoglobulin. The RMSD are also given as vertical bars only for the Antp C39S HD-DNA complex to avoid excessive plot crowding. The corrected and uncorrected LPBE can be easily paired because they obviously span the same range of the x-axis.

Electrostatic potentials inside the molecule

(16)

Electrostatic potentials outside the molecule

View larger version (39K): [in this window] [in a new window]   FIGURE 9   Maps of the electrostatic potential (in kT/q units) from the PBE (upper panel) in a plane orthogonal to the DNA axis going through the center of geometry of the Antp C39S HD-DNA complex at 14.5 mM ionic strength. The difference in absolute value between the potential obtained from the LPBE (middle panel) and the corrected (Eq. 14) LPBE (lower panel) are shown.

View larger version (28K): [in this window] [in a new window]   FIGURE 10   The distribution of the difference (in kT/q units) in absolute value between the PBE potential and the potential obtained from the LPBE and the corrected (Eq. 14) LPBE for the Antp C39S HD-DNA complex and isolated DNA (upper panel), and isolated homeodomain and fully protonated -lactoglobulin (lower panel) at 14.5 mM ionic strength, shown as percentage of points in small intervals of the x-axis. The integral of the distribution is also shown. The uncorrected LPBE curves are recognizable because they span a larger error range.

Electrostatic free energies

	CONCLUSIONS

Top Abstract Introduction Theory Materials and methods Results and discussion Conclusions Appendix A Appendix B References

The LPBE is a widely used approximation of the PBE for biomolecular simulations. Results on SSG show that the potential at the surface obtained using the LPBE may be easily corrected to reproduce fairly accurately that obtained from the PBE. The correction depends on the reduced electric field at the surface, but does not involve any geometrical parameter. Although it might be surprising that systems as different as spheres, cylinders, and planes behave very similarly in this respect, the results may be rationalized taking into account that the deviations from the linearization condition depend on the ratio between the electric field and the Debye scale length and not on geometrical parameters.

Results on high charge density systems show that there are limitations to the applicability of the LPBE to biomolecular systems at low ionic strength (14.5 mM) and to a lesser extent at higher ionic strength (145 mM). For systems in physiological ionic strength with lower charge density the LPBE gives virtually the same results as the PBE. The range of the error in the potential for the LPBE (compared to the PBE) spans few kT/q for the systems studied here. The LPBE can be corrected with a simple formula that does not involve any geometrical parameter, as inferred from the study of SSG. The correction allows for more accurate calculation of the electrostatic free energy of the systems.