Theoretical Background

Water has many biologically necessary properties, one of which is as a dielectric. As a dielectric, water screens (lessens) electrostatic interactions between charged particles. Water can therefore be crudely modeled as a dielectric continuum. In this manner, the electrostatic forces of a biological system can be expressed as a system of differential equations which can be solved for the electric field caused by a collection of charges.

Poisson Boltzmann Equation

The Poisson Boltzmann equation (PBE),

$$\vec{\nabla} \cdot \left[ \epsilon (\vec{r}) \vec{\nabla} \Psi(\v... ..._i \lambda(\vec{r}) \cdot \exp \left[ \frac{-q_i \Psi(\vec{r})}{k_B T} \right]$$

is a nonlinear equation which solves for the electrostatic field, $\Psi(\vec{r})$, based on the position dependent dielectric, $\epsilon(\vec{r})$, the position-dependent accessibility of position $\vec{r}$ to the ions in solution, $\lambda(\vec{r})$, the solute charge distribution, $\rho^f(\vec{r})$, and the bulk charge density, $c_i^\infty$, of ion $q_i$. While this equation does exactly solve for the electrostic field of a charge distribution in a dielectric, it is very expensive to solve, and therefore not suitable for molecular dynamics.

Generalized Born

The Generalized Born (GB) equation is an approximation of the PBE. It models atoms as charged spheres whose internal dielectric is lower than that of the environment. The screening which each atom, $i$, experiences is determined by the local environment; the more atom $i$ is surrounded by other atoms, the less it's electrostatics will be screened since it is more surrounded by low dielectric; this property is called one atom descreening another. Different GB models calculate atomic descreening differently. Descreening is used to calculate the Born radius, $\alpha_i$, of each atom. The Born radius of an atom measures the degree of descreening. A large Born radius represents small screening (strong electric field) as if the atom were in vacuum. A small Born radius represents large screening (weak electric field) as if the atom were in bulk water. The next section describes how the Born radius is calculated and how it is used to calculate electrostatics.

Generalized Born Equations

In a GB simulation, the total electrostatic force on an atom, $i$, is the net Coulomb force on atom $i$ (from nearby atoms) minus the GB force on atom $i$ (also caused by nearby atoms):

$$\vec{F}_i = \vec{F}_i^{Coulomb} - \vec{F}_i^{GB}.$$

Forces are contributed by other nearby atoms within a cutoff. The GB force on atom $i$ is the derivative of the total GB energy with respect to relative atom distances $r_{ij}$,

$$\vec{F}_i^{GB} = -\sum_j \left[ \frac{d E_T^{GB}}{d r_{ij}} \right] \hat{r}_{ji}$$ (8)

$$= -\sum_j \left[ \sum_k \frac{\partial E_T^{GB}}{\partial \alpha_k}... ...}{d r_{ij}} + \frac{\partial E_{ij}^{GB}}{\partial r_{ij}} \right] \hat{r}_{ji}$$ (9)

$$= -\sum_j \left[ \frac{\partial E_T^{GB}}{\partial \alpha_i}\frac{d... ...r_{ij}}+ \frac{\partial E_{ij}^{GB}}{\partial r_{ij}} \right] \hat{r}_{ji} \; .$$ (10)

where the partial derivatives are included since the Born radius, $\alpha$, is a function of all relative atom distances. The total GB energy of the system is

$$E_T^{GB} = \sum_i \sum_{j>i} E_{ij}^{GB} + \sum_{i} E_{ii}^{GB} \; ,$$ (11)

where $E_{ii}^{GB}$ is the Born radius dependent self energy of atom $i$, and the GB energy between atoms $i$ and $j$ is given by

$$E^{GB}_{ij} = - k_e D_{ij} \frac{q_i q_j}{f_{ij}} \; .$$ (12)

The dielectric term [64] is

$$D_{ij} = \left( \frac{1}{\epsilon_p} - \frac{\exp{\left(-\kappa f_{ij}\right)}}{\epsilon_s} \right) \; ,$$ (13)

and the GB function [65] is

$$f_{ij} = \sqrt{r_{ij}^2 + \alpha_i \alpha_j \exp{\left(\frac{-r_{ij}^2}{4 \alpha_i \alpha_j}\right)}} \; .$$ (14)

As the Born radii of atoms $i$ and $j$ decrease (increasing screening), the effective distance between the atoms ($f_{ij}$) increases. The implicit solvent implemented in NAMD is the model of Onufriev, Bashford and Case [53,54] which calculates the Born radius as

$$\alpha_k = \left[ \frac{1}{\rho_{k0}} - \frac{1}{\rho_k}\textrm{tanh}\left(\delta\psi_k - \beta\psi_k^2 + \gamma\psi_k^3\right)\right]^{-1}$$ (15)

where

$$\psi_k = \rho_{k0} \sum_l H_{kl} \; .$$ (16)

$H_{ij}$ is the piecewise descreening function [54,31,59]; the seven piecewise regimes are

$$\textrm{Regimes} = \left\{ \begin{array}{r l l} \textrm{0} & r_{i... ...m{otherwise} &(\textrm{sphere}~j~\textrm{inside~sphere}~j)\\ \end{array}\right.$$ (17)

and the values of $H_{ij}$ are

$$H_{ij} = \left\{ \begin{array}{r l} \textrm{0} & 0 \\ \textrm{I} ... ...ho_{js}-r_{ij}}{r_{ij}+\rho_{js}} \right]\\ \textrm{VI} & 0 \end{array} \right.$$ (18)

Below are defined the derivatives of the above functions which are required for force calculations.

$$\frac{\partial E_{ij}}{\partial r_{ij}} = - k_e \left[ \frac{q_i ... ...\frac{q_i q_j D_{ij}}{f_{ij}^2} \frac{\partial f_{ij}}{\partial r_{ij}} \right]$$ (19)

$$\frac{\partial D_{ij}}{\partial r_{ij}} = \frac{\kappa}{\epsilon_s} \exp{\left(-\kappa f_{ij}\right)\frac{\partial f_{ij}}{\partial r_{ij}}}$$ (20)

$$\frac{\partial f_{ij}}{\partial r_{ij}} = \frac{r_{ij}}{f_{ij}} \... ... - \frac{1}{4} \exp{\left(\frac{-r_{ij}^2}{4 \alpha_i \alpha_j}\right)} \right]$$ (21)

$$\frac{d \alpha_k}{d r_{ij}} = \frac{\alpha_k^2}{\rho_k}\left(1-\t... ...t) \left( \delta - 2\beta\psi_k+3\beta\psi_k^2\right) \frac{d \psi_k}{d r_{ij}}$$ (22)

$$\frac{d \psi_k}{d r_{ij}} = \rho_{k0} \sum_l \frac{d H_{kl}}{d r_{ij}}$$ (23)

$$= \rho_{k0} \sum_l \frac{\partial H_{kl}}{\partial r_{kl}}\frac{d r_{kl}}{d r_{ij}}$$ (24)

$$= \rho_{k0} \left[ \frac{\partial H_{kj}}{\partial r_{kj}}\delta_{ki} + \frac{\partial H_{ki}}{\partial r_{ki}}\delta_{kj} \right]$$ (25)

$$\frac{d \alpha_k}{d r_{ij}} = \frac{\alpha_i^2\rho_{i0}}{\rho_i}\left(1-\textrm{tanh}^2\left(\d... ...psi_i+3\beta\psi_i^2\right) \frac{\partial H_{ij}}{\partial r_{ij}} \delta_{ki}$$

$$+ \frac{\alpha_j^2\rho_{j0}}{\rho_j}\left(1-\textrm{tanh}^2\left(\d... ...psi_j+3\beta\psi_j^2\right) \frac{\partial H_{ji}}{\partial r_{ij}} \delta_{kj}$$ (26)

$$\frac{\partial E_{ij}}{\partial \alpha_i} = -\frac{1}{\alpha_i}\f... ...rac{r_{ij}^2}{4}\right)\exp{\left(\frac{-r_{ij}^2}{4 \alpha_i \alpha_j}\right)}$$ (27)

$$\frac{\partial E_{ij}}{\partial \alpha_j} = -\frac{1}{\alpha_j}\f... ...rac{r_{ij}^2}{4}\right)\exp{\left(\frac{-r_{ij}^2}{4 \alpha_i \alpha_j}\right)}$$ (28)

$$\frac{\partial H_{ij}}{\partial r_{ij}} = \left\{ \begin{array}{r... ...o_{js}-r_{ij}}{r_{ij}+\rho_{js}} \right] \\ \textrm{VI} & 0 \end{array} \right.$$ (29)

Other variables referenced in the above GB equations are $r_{ij}$ - distance between atoms i and j; calculated from atom coordinates.

• $\kappa$ - debye screening length; calculated from ion concentration, $\kappa^{-1} = \sqrt{\frac{\epsilon_0 \epsilon_p k T}{2 N_A e^2 I}}$; $\kappa^{-1} = 10$ Å for 0.1 M monovalent salt.
• $\epsilon_s$ - dielectric constant of solvent.
• $\epsilon_p$ - dielectric constant of protein.
• $\alpha_i$ - Born radius of atom $i$.
• $\rho_i$ - intrinsic radius of atom $i$ taken from Bondi [9].
• $\rho_0$ - intrinsic radius offset; $\rho_0 = 0.09$ Å by default [54].
• $\rho_{i0} = \rho_i - \rho_0$
• $\rho_{is} = \rho_{i0} S_{ij}$
• $S_{ij}$ - atom radius scaling factor [31,64].
• $k_e$ - Coulomb's constant, $\frac{1}{4 \pi \epsilon_0}$, 332.063711 kcal Å / e$^2$.
• $\{\delta, \beta, \gamma\} = \{0.8, 0, 2.91\}~\textrm{or}~\{1.0, 0.8, 4.85\}$ [54]