Implicit solvation

The free energy of solvation of a solute molecule in the simplest ASA-based method is given by:

${\displaystyle \Delta G_{\mathrm {solv} }=\sum _{i}\sigma _{i}\ ASA_{i}}$

where ${\displaystyle ASA_{i}}$ is the accessible surface area of atom i, and ${\displaystyle \sigma _{i}}$ is solvation parameter of atom i, i.e., a contribution to the free energy of solvation of the particular atom i per surface unit area. The needed solvation parameters for different types of atoms (carbon (C), nitrogen (N), oxygen (O), sulfur (S), etc.) are usually determined by a least squares fit of the calculated and experimental transfer free energies for a series of organic compounds. The experimental energies are determined from partition coefficients of these compounds between different solutions or media using standard mole concentrations of the solutes.[5][6]

Notably, solvation energy is the free energy needed to transfer a solute molecule from a solvent to vacuum (gas phase). This energy can supplement the intramolecular energy in vacuum calculated in molecular mechanics. Thus, the needed atomic solvation parameters were initially derived from water-gas partition data.[7] However, the dielectric properties of proteins and lipid bilayers are much more similar to those of nonpolar solvents than to vacuum. Newer parameters have thus been derived from octanol-water partition coefficients[8] or other similar data. Such parameters actually describe transfer energy between two condensed media or the difference of two solvation energies.

Although this equation has solid theoretical justification, it is computationally expensive to calculate without approximations. The Poisson-Boltzmann equation (PB) describes the electrostatic environment of a solute in a solvent containing ions. It can be written in cgs units as:

${\displaystyle {\vec {\nabla }}\cdot \left[\epsilon ({\vec {r}}){\vec {\nabla }}\Psi ({\vec {r}})\right]=-4\pi \rho ^{f}({\vec {r}})-4\pi \sum _{i}c_{i}^{\infty }z_{i}q\lambda ({\vec {r}})e^{\frac {-z_{i}q\Psi ({\vec {r}})}{kT}}}$

or (in mks):

${\displaystyle {\vec {\nabla }}\cdot \left[\epsilon ({\vec {r}}){\vec {\nabla }}\Psi ({\vec {r}})\right]=-\rho ^{f}({\vec {r}})-\sum _{i}c_{i}^{\infty }z_{i}q\lambda ({\vec {r}})e^{\frac {-z_{i}q\Psi ({\vec {r}})}{kT}}}$

where ${\displaystyle \epsilon ({\vec {r}})}$ represents the position-dependent dielectric, ${\displaystyle \Psi ({\vec {r}})}$ represents the electrostatic potential, ${\displaystyle \rho ^{f}({\vec {r}})}$ represents the charge density of the solute, ${\displaystyle c_{i}^{\infty }}$ represents the concentration of the ion i at a distance of infinity from the solute, ${\displaystyle z_{i}}$ is the valence of the ion, q is the charge of a proton, k is the Boltzmann constant, T is the temperature, and ${\displaystyle \lambda ({\vec {r}})}$ is a factor for the position-dependent accessibility of position r to the ions in solution (often set to uniformly 1). If the potential is not large, the equation can be linearized to be solved more efficiently.[9]

A number of numerical Poisson-Boltzmann equation solvers of varying generality and efficiency have been developed,[10][11][12] including one application with a specialized computer hardware platform.[13] However, performance from PB solvers does not yet equal that from the more commonly used generalized Born approximation.[14]

The Generalized Born (GB) model is an approximation to the exact (linearized) Poisson-Boltzmann equation. It is based on modeling the solute as a set of spheres whose internal dielectric constant differs from the external solvent. The model has the following functional form:

${\displaystyle G_{s}=-{\frac {1}{8\pi \epsilon _{0}}}\left(1-{\frac {1}{\epsilon }}\right)\sum _{i,j}^{N}{\frac {q_{i}q_{j}}{f_{GB}}}}$

where

${\displaystyle f_{GB}={\sqrt {r_{ij}^{2}+a_{ij}^{2}e^{-D}}}}$

and ${\displaystyle D=\left({\frac {r_{ij}}{2a_{ij}}}\right)^{2},a_{ij}={\sqrt {a_{i}a_{j}}}}$

where ${\displaystyle \epsilon _{0}}$ is the permittivity of free space, ${\displaystyle \epsilon }$ is the dielectric constant of the solvent being modeled, ${\displaystyle q_{i}}$ is the electrostatic charge on particle i, ${\displaystyle r_{ij}}$ is the distance between particles i and j, and ${\displaystyle a_{i}}$ is a quantity (with the dimension of length) termed the effective Born radius.[15] The effective Born radius of an atom characterizes its degree of burial inside the solute; qualitatively it can be thought of as the distance from the atom to the molecular surface. Accurate estimation of the effective Born radii is critical for the GB model.[16]

Another possibility is to use ad hoc quick strategies to estimate solvation free energy. A first generation of fast implicit solvents is based on the calculation of a per-atom solvent accessible surface area. For each of group of atom types, a different parameter scales its contribution to solvation ("ASA-based model" described above).[19]

Another strategy is implemented for the CHARMM19 force-field and is called EEF1.[20] EEF1 is based on a Gaussian-shaped solvent exclusion. The solvation free energy is

${\displaystyle \Delta G_{i}^{solv}=\Delta G_{i}^{ref}-\sum _{j}\int _{Vj}f_{i}(r)dr}$

The reference solvation free energy of i corresponds to a suitably chosen small molecule in which group i is essentially fully solvent-exposed. The integral is over the volume V_j of group j and the summation is over all groups j around i. EEF1 additionally uses a distance-dependent (non-constant) dielectric, and ionic side-chains of proteins are simply neutralized. It is only 50% slower than a vacuum simulation. This model was later augmented with the hydrophobic effect and called Charmm19/SASA.[21]

It is possible to include a layer or sphere of water molecules around the solute, and model the bulk with an implicit solvent. Such an approach is proposed by M. J. Frisch and coworkers[22] and by other authors.[23][24] For instance in Ref.[23] the bulk solvent is modeled with a Generalized Born approach and the multi-grid method used for Coulombic pairwise particle interactions. It is reported to be faster than a full explicit solvent simulation with the particle mesh Ewald summation (PME) method of electrostatic calculation. There are a range of hybrid methods available capable of accessing and acquiring information on solvation.[25]

All implicit solvation models rest on the simple idea that nonpolar atoms of a solute tend to cluster together or occupy nonpolar media, whereas polar and charged groups of the solute tend to remain in water. However, it is important to properly balance the opposite energy contributions from different types of atoms. Several important points have been discussed and investigated over the years.