The accurate modeling of molecules in solution using molecular mechanics requires realistic models for the interaction of the solvent with the molecule [1-3]. To treat such a medium (usually water) in a molecular calculation, both explicit and implicit models have been developed [4].

Molecular dynamics simulation of a small protein, the N terminal domain of Troponin C.
(A) A view of the inside of the initial box for a simulation which includes water molecules explicitly;
the whole system consits of 1229 atoms in the protein and 6009 atoms for water.
(B) In a implicit solvent simulation, the effect of water is included in a potential of mean force, W(X), vizualized as arrows.
In this case, the simulation only includes the 1229 protein atoms.

1. Explicit Solvent Models

Explicit solvent models rely on using hundreds or thousands of discrete solvent molecules [5]. They are the most widely used methods for carrying out simulations in solvent. Such calculations converge only slowly because of the large number of particles involved. They generally require orders of magnitude more CPU time than corresponding gas phase calculations on the same molecule. Because explicit models are so computationally demanding, there is a significant interest in developing more rapid implicit solvent models.
 

2. Implicit Solvent Models

Implicit models treat the solvent as a continuous medium having the average properties of the real solvent, and surrounding the solute beginning at the van der Waals surface. A variety of continuum models have been described over the years [6-10]. Among these, the generalized Born (GB), Surface Area (SA) model has become very popular [6].
 

2.1 The GB/SA model

In the GB/SA model, the total solvation free energy G_solv is given as the sum of a solvent-solvent cavity term (G_cav), a solute-solvent van der Waals term (G_vdW), and a solute-solvent electrostatics polarization term (G_pol) [6]:

                                                                                                                                             (1)

The GB/SA model computes G_cav + G_vdW together as a linear function of the solvent-accessible surface areas:

                                                                                                                           (2)

where ASA(k) is the total solvent accessible surface area of atom k and ASP(k) is an empirically determined atomic solvation parameter.

Implementation of equation (2) in a molecular mechanics or molecular dynamics simulation requires that the accessible surface area of all atoms as well as their derivatives with respect to the position of all atoms be computed.

The generalized Born equation defines G_pol as [11]:

                                                                                  (3)

where the double sum runs over all pairs of atoms (i and j). q_i and q_j are the partial charges of i and j, respectively, and r_ij is the i,jth atom pair separation.  is the dielectric constant of the medium. b_i and b_j are the so-called Born radius of atom i and j, respectively.
 

2.2 Computing the Born radius of an atom

In the case of a simple ion of radius a and charge q, the electrostatic component of the solvation free energy can be found analytically and the result is the well-known Born formula [12]:

                                                                                                                                                                                 (4)

where  is the dielectric constant of the solvent. The radius a is also known as the Born radius of the ion.

For an atom i in a molecule, the situation is more complicated: it also interacts with the solvent, but part of this solvent has now been replaced by other atoms of the molecule, which are represented by first approximation as sphere. The basic idea of the Generalized Born model is to define an effective radius for i, R_i such that the electrostatics contribution of the solvation free energy of i is given by equation (4). R_i generally depends not only on a_i, the intrinsic atomic radius of i, but also on the radii and relative positions of all other atoms.

Ideally the Born radius should be chosen so that if one were to solve the Poisson equation for a single charge q_i placed at the position of atom i, and a dielectric boundary determined by all of the molecule's atoms and their radii, then the self-energy of charge i would be equal to its Born energy. Obvioulsy, the procedure per se would have no practical advantage over a direct calculation using a numerical solution of the Poisson equation.

To find a more rapidly calculable approximation for the Born radius, approximations based on energy densities have been developped. Using the Coulomb field approximation, an integral formula was proposed for the Born radius:

                                                                                                                                             (5)

where in stands for the interior region of the molecule, excluding a radius a_i around the center of atom i (for more detailed information, see [6]).

Computing the Born radius R of an atom is then equivalent to estimating the integral in equation (5). Three different approximations have been proposed:
 
 

• 2.2.1 The Overlapping Spheres Approach
   

  In the pair-wise versions of the GB theory, the basic idea is to approximate the integral in equation (5) by a sum of contributions for each atom. If the molecule consisted of a set of non-overlapping spheres of radius a_j at position r_ij relative to atom i, then equation 5 could be written as a sum of integrals over spherical volumes,

                                                                                                                       (6)

  The integrals over spheres can then be calculated analytically, leading to:

                                                                               (7)

  A straightforward pair-wise summation as described by equation (7) would overcount the solute region since neighbouring atoms overlap with each other. Hawkins el al [13,14] proposed scaling the neighboring values of R_i as an empirical correction to compensate for this neglect of overlap. Their expression for the GB radiis is:

                                                                                                                               (8)

  where H is a fairly complex expression and the S_j scaling factors, fit either to experiment or to numerical Poisson results. Several groups have adopted this idea, using different training sets to determine how best to scale the neighbouring radii [13-16].
   
   

• 2.2.2 Asymptotic Approach

 

The long distance limit of the integrals over spheres in equation (6) is just V_jr_ij^-4, where V_j is the volume of the jth atom. The direct use of this functional form for pair-wise atom contributions would lead to significant errors in the short distance range. Qiu et al [17] have modified this idea, scaling each atomic contribution by a factor that depends on the number of covalent bonds between atom j and atom i, so that the sum over neighbouring atoms becomes:

                                                                                         (9)

the P_k are adjustable parameters and C is a close-contact function that adjusts radii for non-bonded atoms that are very close to the central atom i. This approach has alos been adopted by several groups [17-19].
 
 

• 2.2.3 The Analytical Continuum Electrostatics Model (ACE)

 

Schaefer and Karplus [20] presented a general formalism for decomposing energy functions based on integration of the Couilomb field energy density into pair-wise atomic terms. The integration over the solute interior in equation (5) can be rewritten as an integral over space with the integrand multiplied by a Step function P(r) whose value is 1 in the molecular interior and zero elsewhere. This function can then be written as a sum of atomic terms,

                                                                                                                                                         (10)

where, for example, the P_j might be step function corresponding to the Voronoi volumes of the atoms. Schaefer and Karplus [20] proposed a Gaussian form for the P_j, normalized according to the effective volume of each atom. Based on this idea, they developed an analytical, continuous, and differentiable pair-wise atomic expression for the electrostatics energy of the solute, called analytical treatment of continuum electrostatics (ACE).

It is not yet clear if there are intrinsic differences among the models that would systematically favor one over the others.

References

1. Cheatham, TE and Kollman, PA. Molecular dynamics simulation of nucleic acids. Ann. Rev. Phys. Chem., 51, 435-471 (2000).

2. Makarov, V, Pettit, BM and Feig, M. Solvation and hydration of proteins and nucleic acids: A theoretical view of simulation and experiment. Accounts of Chem. Res., 35, 376-384 (2002).

3. Karplus, M. Molecular dynamics simulations of biomolecules. Accounts of Chem. Res., 35, 321-323 (2002).

4. Xia, B, Tsui, V, Case, DA, Dyson, HJ and Wright, PE. Comparison of protein solution structures refined b molecular dynamics simulation in vacuum, with a generalized Born model, and with explicit water. J. Biol. NMR, 22, 317-331 (2002).

5. Bizzarri, AR and Cannistraro, S. Molecular dynamics of water at the protein-solvent interface. J. Phys. Chem. B, 106, 6617-6633 (2002).

6. Tsui, V and Case DA. Theory and applications of the generalized Born solvation model in macromolecular simulations. Biopolymers, 56, 275-291 (200).

7. Simonson, T. Macromolecular electrostatics: continuum models and their growing pains. Curr. Opin. Struct. Biol. , 11, 243-252 (2001).

8. Hassan, SA and Mehler, EL. A critical analysis of continuum electrostatics: the screened Coulomb potential-implicit solvent model and the study of the alanine dipeptide and discrimination of misfolded structures of proteins. Proteins: Struct. Funct. Genet., 47, 45-61 (2002).

9. Lee, MS, Salsbury, FR and Brooks, CL. Novel generalized Born methods. J. Chem. Phys., 116, 10606-10614 (2002).

10. Orzco, M and Luque, FJ. Theoretical methods for the description of the solvent effect in biomolecular systems. Chem. Rev., 100, 4187-4225 (2000).

11. Still, WC, Tempczyk, A, Hawley, RC, Hendrickson, T. Semianalytical treatment of solvation for molecular mechanics and dynamics. J. Am. Chem. Soc., 112, 6127-6129 (1991).

12. Born, M. Z. Phys., 1, 45-48 (1920).

13. Hawkins, GD, Cramer, CJ and Truhlar, DG. Pairwise solute descreening of solute charges from a dielectric medium. Chem. Phys. Lett., 246, 122-129 (1995).

14. Hawkins, GD, Cramer, CJ and Truhlar, DG. Parametrized models of aqueous free energies of solvation based on pairwise descreening of solute atomic charges from a dielectric medium. J. Phys. Chem., 100, 19824-19839 (1995).

15. Dudek, MJ and Ponder, JW. Accurate modeling of the intramolecular electrostatic energy of proteins. J. Comput. Chem., 16, 791-816 (1995).

16. Tsui, V and Case, DA. Molecular dynamics simulations of nucleic acids with a generalized born solvation model. J. Am. Chem. Soc., 122, 2489-2498 (2000).

17. Qiu, D, Shenkin, PS, Hollinger, FP and Still, WC. The GB/SA continuum model for solvation. A fast analytical method for the calculation of approximate Born radii. J. Phys. Chem. A, 101, 3005-3014 (1997).

18. Dominy, BN and Brooks, CL. Development of a generalized Born model parametrization for proteins and nucleic acids. J. Phys. Chem. B, 103, 3765-3773 (1999).

19. Zhu, J, Shi, Y and Liu, H. Parametrization of a generalized Born/solvent accessible surface area model and applications to the simulation of protein dynamics. J. Phys. Chem. B, 106, 4844-4853 (2002).

20. Schaefer, M and Karplus, M. A comprehensive analytical treatment of continuum electrostatics. J. Phys. Chem., 100, 1578-1599 (1996).