Ammm:Gb gk

Why Implicit Solvent Model?

Solvation effect

• Solvent effect

Explicit vs. Implicit Solvent models

Table Explicit vs. Implicit Solvent Model

Models	Explicit Solvent Model	Implicit Solvent Model
Idea	Represent solvent molecules explicitly.	Respresent solvent as a continuous medium.
Representation
Advantage	High level of details is provided; Generally more accurate.	Computationally efficient.

Main approaches to implicit models

Implicit models treat the solvent as a continuous medium having the average properties of the real solvent. A variety of continuum models have been described over the years .

The basic idea is

${G}_{solv}={G}_{nonpolar}+{G}_{pol}$

Among these, the generalized Born (GB)/ Surface Area (SA) model has become very popular .

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}_{disp}$ ), and a solute-solvent electrostatics (polar) term (${G}_{pol}$ ):

${G}_{solv}={G}_{cav}+{G}_{disp}+{G}_{pol}$ (1)

The GB/SA model computes ${G}_{cav}+{G}_{disp}$ together as a linear function of the solvent-accessible surface areas:

${G}_{cav}+{G}_{disp}=\gamma \cdot SA(k)$ (2) where SA(k) is the total solvent accessible surface area of atom k and $\gamma$ is an empirically determined atomic solvation parameter.

When ions are solvated by polar solvents, Gpol dominates; when polar molecules are solvated by polar solvents, Gpol, Gdisp and Gcav are of similar magnitude. Usually Gpol and Gdisp are negative; Gcav is positive.

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 electrostatic (polar)component of solvation free energy can be calculated:

${{G}_{pol}}=-{\frac {1}{2}}\left(1-{\frac {1}{\epsilon }}\right)\sum _{i=1}^{N}\sum _{j=1}^{N}{\frac {{q}_{i}{q}_{j}}{\sqrt {{{r}_{ij}}^{2}+{R}_{i}{R}_{j}exp\left(-{\frac {{{r}_{ij}}^{2}}{4{R}_{i}{R}_{j}}}\right)}}}$

The following expression is the same as the previous one, with the first term called self-energies term and the second is called cross-energies term:

${{G}_{pol}}=-{\frac {1}{2}}\left(1-{\frac {1}{\epsilon }}\right)\sum _{i=1}^{N}{\frac {{{q}_{i}}^{2}}{{R}_{i}}}-{\frac {1}{2}}\left(1-{\frac {1}{\epsilon }}\right)\sum _{i=1}^{N}\sum _{j\neq i}^{N}{\frac {{q}_{i}{q}_{j}}{\sqrt {{{r}_{ij}}^{2}+{R}_{i}{R}_{j}exp\left(-{\frac {{{r}_{ij}}^{2}}{4{R}_{i}{R}_{j}}}\right)}}}$

Models to calculate the electrostatic(polar) comonent of solvation free energy

• Born in 1920 described the electrostatic solvation energy of a charged, sphereical ion in terms of macroscopic continuum theory.
• In 1934, Kirkwood extended the approach to a spherical particle with arbitrary electrostatic multipole moments.
• Around 1990's, Generalized Born approach has been derived and applied to estimate the electrostatic contribution of solvation free energy.
• In 1997, Kong and Ponder revisited Kirkwood's theory to allow analytic treatment of off-center point multipoles.
• In 2007, Schnieders and Ponder applied the Generalized Kirkwood model to the newly developed AMOEBA force field, which includes permanent atomic multipoles through the quadrupole and treats polarization via induced dipole.

Born Model and Kirkwood Model

Born Model

Derivation of the Born Model

Consider the free energy of placing a charge q at the center of a spherical cavity in a solvent.

The electrostatic energy can be calculatd as

$W={\frac {1}{8\pi }}\int _{V}^{}D\cdot EdV={\frac {1}{2}}\int _{R}^{+inf}{\varepsilon }_{0}\varepsilon \left({\frac {q}{{\varepsilon }_{0}\varepsilon {r}^{2}}}\right)^{2}{r}^{2}dr={\frac {1}{2}}{\frac {{q}^{2}}{R\varepsilon }}$

Therefore Free energy can be calculated using

${G}_{elec}=W\left(\varepsilon =\epsilon \right)-W\left(\varepsilon =1\right)={\frac {1}{2}}{\frac {{q}^{2}}{R\epsilon }}-{\frac {1}{2}}{\frac {{q}^{2}}{R}}=-{\frac {{q}^{2}}{2R}}(1-{\frac {1}{\varepsilon }})$

Calculation of the Born Radii

Computing the Born radius of an atom

In the case of a simple ion of radius R and charge q, the electrostatic component of the solvation free energy is:

${{G}_{Born}}=-{\frac {{q}^{2}}{2R}}(1-{\frac {1}{\varepsilon }})$ (4)

where $\varepsilon$ is the dielectric constant of the solvent. The radius R 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, Ri such that the electrostatics contribution of the solvation free energy of i is given by equation (4). Ri generally depends not only on ai, 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 qi 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.

Considering the work of placing a charge qi at the orgin within a molecule whose interior dielectric constant is ${\varepsilon }_{in}$, surrounded by a medium of dielectric constant ${\varepsilon }_{ex}$and in which no other charges have yet been placed,

$W={\frac {1}{8\pi }}\int _{V}^{}D\cdot EdV\approx {\frac {1}{8\pi }}\int _{in}^{}{\frac {{{q}_{i}}^{2}}{{r}^{4}{\varepsilon }_{in}}}dV+{\frac {1}{8\pi }}\int _{ex}^{}{\frac {{{q}_{i}}^{2}}{{r}^{4}{\varepsilon }_{ex}}}dV$

Taking the difference in W when ${\varepsilon }_{ex}$ is changed from 1 to ${\varepsilon }_{w}$,

${\Delta G}_{{pol}_{i}}=-{\frac {1}{8\pi }}\left(1-{\frac {1}{{\varepsilon }_{w}}}\right)\int _{ex}^{}{\frac {{{q}_{i}}^{2}}{{r}^{4}{\varepsilon }_{ex}}}dV$

${{R}_{i}}^{-1}={\frac {1}{4\pi }}\int _{ex}^{}{\frac {1}{{r}^{4}}}dV$ (5)

Rewrite the eq(5) in terms of an integration over the interior region, excluding a radius ai around the origin,

${{R}_{i}}^{-1}={{a}_{i}}^{-1}-{\frac {1}{4\pi }}\int _{in,r>{a}_{i}}^{}{\frac {dV}{{r}^{4}}}$ (6)

where in stands for the interior region of the molecule, excluding a radius ai around the center of atom i.

Note that in the case of a monoatomic ion, where the molecular boundary is simply the sphere of radius ai, the equation becomes Ri=ai and the Born formula is recovered exactly. For a roughly spherical molecule with the atom i in its center, the Born radius Ri should be approxiately equal to Li.

Computing the Born radius R of an atom is then equivalent to estimating the integral in equation (6). Many different approximations have been proposed with the following implemented in TINKER: STILL,ONION, HCT, OBC, ACE, etc.

The Overlapping Spheres Approach

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

${{R}_{i}}^{-1}={{a}_{i}}^{-1}-{\frac {1}{4\pi }}\sum _{j}^{}\int _{Spherej}^{}{\frac {dV}{{r}^{4}}}$ (7)

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

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle {{R}_{i}}^{-1}={{a}_{i}}^{-1}-\sum_{j}^{}\left(\frac{{a}_{j}}{2\left({{r}_{ij}}^{2}-{{a}_{j}}^{2} \right)} -\frac{1}{4{r}_{ij}}\log\frac{{r}_{ij}-{a}_{j}}{{r}_{ij}+{a}_{j}}\right)} (8)

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

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle {{R}_{i}}^{-1}={{a}_{i}}^{-1}-\sum_{j}^{}H({r}_{ij},{S}_{j}{a}_{j})} (9)

where H is a fairly complex expression and the Sj 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 radus .

Kirkwood Model

Derivation of the Kirkwood Model

Kirkwood Expansion

When n=0,

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle {G}_{ion}=\frac{1}{2}q{\Phi }_{RF}(r, n=0)=\frac{{q}^{2}\left({\varepsilon }_{in}-{\varepsilon }_{out} \right)}{8\pi {\varepsilon }_{in}a}\frac{0+1}{0+{\varepsilon }_{out}}=\frac{{q}^{2}\left(1-\varepsilon \right)}{8\pi \varepsilon a}=-\frac{{q}^{2}}{8\pi a}\left(1-\frac{1}{\varepsilon } \right)}

This G is known as the Born energy.

When n=1, the energy of a dipole in its reaction field (the negative gradient of the reaction potential) is minus one half of the dot product of the dipole and the reaction field.

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle {G}_{ion}=-\frac{1}{2}\mathit{\mu}\cdot {E }_{RF}(r, n=1)=\frac{q\left({\varepsilon }_{in}-{\varepsilon }_{out} \right)}{8\pi {\varepsilon }_{in}a}\frac{1+1}{{\varepsilon }_{in}+{2\varepsilon }_{out}}\frac{{r}_{s}cos\theta\mu }{{a}^{2}} =-\frac{1}{8\pi }\left[\frac{2\left(\varepsilon -1 \right)}{\left(2\epsilon +1 \right)} \right]\frac{{\mu }^{2}}{{a}^{3}}}

Several names (Bell, Onsager, Kirkwood) have been associated with this energy( when dipole Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mu} occupies the sphere a).

Generalized Born and Generalized Kirkwood

Generalized Born

The generalized Born equation defines Gele as :

where qi and qj are the partial charges of i and j, respectively, and rij is the i,jth atom pair separation. is the dielectric constant of the medium. Ri and Rj are the so-called Born radius of atom i and j, respectively.

Generalized Kirkwood

According to the Reaction Potential Approximation

The reaction potential at r due to an off-center charge at r inside a spherical dielectric cavity of permittivity  Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle {\varepsilon }_{in}} surrounded by solvent with permittivity  Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle {\varepsilon }_{out}} is given by Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Phi(r)=\frac{q}{a{\varepsilon }_{in}}\sum_{l=0}^{inf}\frac{\left(l+1 \right)\left({\varepsilon }_{in}-{\varepsilon }_{out} \right)}{{\left(l+1 \right)\varepsilon }_{out}+l{\varepsilon }_{in}}{\left(\frac{r{r}_{0}}{{a}^{2}} \right)}^{l} {P}_{l}\left(cos\theta \right)} (10)

where a is the cavity radius, q is the magnitude of the charge, and P is the Legendre polynomial The self-energy of a charge based on eq (10) is Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle W(d)=\frac{{q}^{2}}{2a{\varepsilon }_{in}}\sum_{l=0}^{inf}\frac{\left(l+1 \right)\left({\varepsilon }_{in}-{\varepsilon }_{out} \right)}{{\left(l+1 \right)\varepsilon }_{out}+l{\varepsilon }_{in}}{\left(\frac{{d}^{2}}{{a}^{2}} \right)}^{l}{P}_{l}(1) }

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \approx \frac{1}{2} \left(\frac{1}{{\varepsilon }_{out}}-\frac{1}{{\varepsilon }_{in}} \right)\frac{{q}^{2}}{a}\sum_{l=0}^{inf} {\left(\frac{{d}^{2}}{{a}^{2}} \right)}^{l}=\frac{1}{2} \left(\frac{1}{{\varepsilon }_{out}}-\frac{1}{{\varepsilon }_{in}} \right){q}^{2}\frac{a}{\left({a}^{2}-{d}^{2} \right)}}

As shown by Grycuk, it is possible to calculate the factor Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle {a}_{r}=\frac{a}{\left({a}^{2}-{d}^{2} \right)}} , which is equivalent to the inverse of an effective radius.

After determining effective radii, the self-energy for each permanent atomic multipole can be evaluated using

For the cross-energies term,

Limitations of the GB Model

• Coulomb field approximation Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle {D}_{i}\approx \frac{{q}_{i}\mathbf{r}}{{r}^{3}}} is the main source of its deviation from solvation energies calculated using solutions of the Poisson-Boltzmann equation. If the dielectric boundary between solvent and solute is sharp, then D can be thought of as arising from an induced surface charge density on the dielectric boundary.

• In some cases, instead of Coulomb field approximation, an image-charge approximation for D can recover better energetic behavior.

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

GK With AMOEBA force field

A little demo using tinker and AMOEBA

TINKER 5.0 provides a selection of continuum solvation treatments including several variants of the generalized Born (GB/SA) model and the generalized Kirkwood implicit solvation model for AMOEBA.

References