Energy Function

Energy calculation is at the core of PELE. It is performed during minimization, ANM, ligand perturbation, side chain prediction and the Metropolis evaluation. Though different approximations are considered in each phase, the function we describe here is the general one. Notice that this section does not yet discuss the binding energy that PELE provides.

The general energy function in PELE has the form:

\[E = E_{\text{MM}} + \Delta G_{\text{solv}} + E_{\text{constraints}}\]

The first term is the molecular mechanics component, and corresponds to the energy of the system in vacuum. Depending on the force field chosen for the simulation (see ForceField), this term will have different functional forms and parameters. The second term is the solvation energy for the system, and depends on the solvation model chosen (see Solvent); of course, this term is null in vacuum studies. The third term depends on the constraints fixed by the user (for example, to keep the ligand at the perturbation target position).

The molecular mechanics term is further divided in bonded and non-bonded energy:

\[E_{\text{MM}} = E_{\text{bnd}} + E_{\text{nb}}\]

The bonded part includes the bond, angle, torsion, improper torsion, and 1-4 interactions terms:

\[E_{\text{bnd}} = E_{\text{bond}} + E_{\text{angle}} + E_{\text{torsion}} + E_{\text{improper torsion}} + E_{1-4}\]

As for the non-bonded energy, when using the multi-scale energy algorithm (see NonBonding in Pele++ control file), and for efficiency reasons, it is calculated only for pairs of neighboring atoms. These neighbors are grouped as short-distance neighbors and long-distance neighbors, and each group is itself divided into two groups, depending on whether all non-bonded interactions are calculated, or only the van der Waals ones. For simplicity, in this section we will consider the non-multiscale algorithm where all interactions are calculated, and therefore:

\[E_{\text{nb}} = E_{\text{vdw}} + E_{\text{ele}}\]

Finally, regarding the solvation energy, it consists of one term for the electrostatic or polar part of that energy, and another for the non-polar part; some solvent models also add some corrections to the electrostatic part as penalty terms:

\[\Delta G_{\text{solv}} = \Delta G_{\text{pol}} + \Delta G_{\text{penalty}} + \Delta G_{\text{npol}}\]

The expanded energy formula is then:

\[E = E_{\text{bond}} + E_{\text{angle}} + E_{\text{torsion}} + E_{\text{improper torsion}} + E_{1-4} + E_{\text{vdw}} + E_{\text{ele}} + \Delta G_{\text{solv,pol}} + \Delta G_{\text{solv,penalty}} + \Delta G_{\text{solv,npol}} + E_{\text{constraints}}\]

and with all terms expanded, for the SGB/NP model:

\[\begin{split}E =& \sum_{\text{bonds}} K_{r}(r - r_{\text{eq}})^{2} + \sum_{\text{angles}} K_{\theta} (\theta - \theta_{\text{eq}})^{2} + \sum_{\text{angle}\, i} \sum_{j} {k_{j,1}^{i}[1 + k_{j,2} \cos (n_{j}\phi_{i})]} \\ &+ \sum_{i < j} \left[ \text{scale}_{\text{coulomb}}\frac{q_i^{s} q_j^{s}}{\varepsilon_{\text{in}(ij)}r_{ij}} + \epsilon_{i}^{s}\epsilon_{j}^{s} \left( \left( \frac{\sigma_{ij}^{s}}{r_{ij}} \right)^{12} - \left( \frac{\sigma_{ij}^{s}}{r_{ij}}\right)^{6} \right) \right] f_{ij} \\ &- \sum_{i<j} {(\frac{1}{\varepsilon_{\text{in}(ij)}} - \frac{e^{-\kappa f_{\text{GB}}}}{\varepsilon_{\text{solv}}}) \frac{q_{i}^{s}q_{j}^{s}}{f_{\text{GB}}}} - \frac{1}{2} \sum_{i} {(\frac{1}{\varepsilon_{\text{in}(ii)}} - \frac{e^{-\kappa \alpha_{i}}}{\varepsilon_{\text{solv}}}) \frac{(q_{i}^{s})^2}{f_{\alpha_{i}}}} + \Delta G_{\text{solv,penalty}} + \sum_{i} {C_{\gamma,i} \text{SASA}_{i} + C_{\alpha,i} S(a/\alpha_{i})} \\ &+ E_{\text{constraints}}\end{split}\]

and for the OBC/ACE model:

\[\begin{split}E =& \sum_{\text{bonds}} K_{r}(r - r_{\text{eq}})^{2} + \sum_{\text{angles}} K_{\theta} (\theta - \theta_{\text{eq}})^{2} + \sum_{\text{angle}\, i} \sum_{j} {k_{j,1}^{i}[1 + k_{j,2} \cos (n_{j}\phi_{i})]} \\ &+ \sum_{i < j} \left[ \text{scale}_{\text{coulomb}}\frac{q_i^{s} q_j^{s}}{\varepsilon_{\text{in}(ij)}r_{ij}} + \epsilon_{i}^{s}\epsilon_{j}^{s} \left( \left( \frac{\sigma_{ij}^{s}}{r_{ij}} \right)^{12} - \left( \frac{\sigma_{ij}^{s}}{r_{ij}}\right)^{6} \right) \right] f_{ij} \\ &- \sum_{i<j} {(\frac{1}{\varepsilon_{\text{in}(ij)}} - \frac{e^{-\kappa f_{\text{GB}}}}{\varepsilon_{\text{solv}}}) \frac{q_{i}^{s}q_{j}^{s}}{f_{\text{GB}}}} - \frac{1}{2} \sum_{i} {(\frac{1}{\varepsilon_{\text{in}(ii)}} - \frac{e^{-\kappa \alpha_{i}}}{\varepsilon_{\text{solv}}}) \frac{(q_{i}^{s})^2}{f_{\alpha_{i}}}} + 4\pi \sum_{i} {\sigma_{i} (R_{i} + R_{s})^{2} \left( \frac{R_{i}}{\alpha_{i}} \right)^{6}} \\ &+ E_{\text{constraints}}\end{split}\]

The Molecular Mechanics energy terms

The different force fields used in PELE (see ForceField) all use the same functional form, though they differ in the parameters, scaling factors, and in some of the expansion terms (such as \(E_{\text{torsion}}\)).

Here come the references for the force fields available in PELE. For further checks, consult Schrödinger’s IMPACT, [hornak:2006], Chimera addcharge documentation and the AMBER manual.

• OPLS2005: [banks:2005] for ligands (it is the citing paper for OPLS2005 according to Schrödinger, see http://www.schrodinger.com/kb/546 ; in Maestro, the OPLS2005 parameters are at $SCHRODINGER/mmshare-vversion/data/f14/) and [kaminski:2001] for proteins (affects torsional parameters and, for cysteine and methionine, some nonbonded parameters). Original atom types for OPLS are those of AMBER ([cornell:1995], table 1). The original OPLS family of force fields is based on [jorgensen:1996] (for carbohydrates there was another publication: [damm:1997]) and [kaminski:2001]. The first publication on OPLS, for united atom, was [jorgensen:1988].

• AMBER99: [cornell:1995] (the base of the force field) + [wang:2000] (changes in torsional parameters). Not thoroughly tested in production in PELE, but should work in principle.

• AMBER99SB: AMBER99 + [hornak:2006] (corrections for protein backbone dihedral terms). Not thoroughly tested in production in PELE, but should work in principle.

• AMBER99SBBSC0: AMBER99SB + [perez:2007] (corrections for nucleic acids). This is the AMBER force field used in production in PELE, though most simulations have been with DNA until now.

Following are the different terms, with comments depending on the actual force field used.

Bond and angle terms are the same in all force fields, only differing in the parameters used.

\[ \begin{align}\begin{aligned}E_{\text{bond}} = \sum_{\text{bonds}} K_{r}(r - r_{\text{eq}})^{2}\\E_{\text{angle}} = \sum_{\text{angles}} K_{\theta} (\theta - \theta_{\text{eq}})^{2}\end{aligned}\end{align} \]

The AMBER force fields use this formula for the torsional term:

\[E_{\text{torsion,AMBER}} = \sum_{dihedrals} {\frac{V_{n}}{2}[1 + \cos (n \phi_{i} - \gamma_{n})]}\]

while OPLS uses this other one:

\[E_{\text{torsion,OPLS}} = \sum_{i} {(\frac{V_{1}^{i}}{2}[1 + \cos (\phi_{i} + f_1^{i})] + \frac{V_{2}^{i}}{2}[1 - \cos (2\phi_{i} + f_2^{i})] + \frac{V_{3}^{i}}{2}[1 + \cos (3\phi_{i} + f_3^{i})])}\]

As stated in p. 11227 of [jorgensen:1996], the phase angles \(f_1^{i}\), \(f_2^{i}\) and \(f_3^{i}\) are zero.

To accomodate both types, the actual PELE function is like this:

\[E_{\text{torsion}} = \sum_{\text{angle}\, i} \sum_{j} {k_{j,1}^{i}[1 + k_{j,2} \cos (n_{j}\phi_{i})]}\]

Apparently, the values of \(\gamma_{n}\) or \(f_{n}^{i}\) can only be \(0^{\circ}\) or \(180^{\circ}\), so they are actually modelled by the prefactor value \(\pm 1\) (\(k_{j,2}\)) in the parameter files (since \(\cos (\pi - \phi) = - \cos (\phi)\)). Notice that \(n_{j}\) can only take values 1,2,3,4 and 6. There seem to be force fields or special dihedral atoms for which \(E_{\text{torsion}}\) has more than 3 terms, while this was the number of terms in the original OPLS force field.

Also notice that \(E_{\text{improper torsion}}\) uses the same formula.

For the non bonded term (including the 1-4 interactions), AMBER uses:

\[ \begin{align}\begin{aligned}E_{\text{nb,AMBER}} = \sum_{i < j}[ \text{scale}_{\text{coulomb}}\frac{q_i q_j}{\varepsilon r_{ij}} + \text{scale}_{\text{VDW}}(\frac{A_{ij}}{r_{ij}^{12}} - \frac{B_{ij}}{r_{ij}^{6}})]f_{ij}\\A_{ij} = 4\epsilon_{ij} \sigma_{ij}^{12}\\B_{ij} = 4\epsilon_{ij} \sigma_{ij}^{6}\\\sigma_{ij} = \sqrt {\sigma_{ii} \sigma_{jj}}\\\epsilon_{ij} = \sqrt {\epsilon_{ii} \epsilon_{jj}}\\\text{scale}_{\text{coulomb}} = \frac{1}{1.2} \quad \text{1-4 pairs}\\\text{scale}_{\text{coulomb}} = 1.0 \quad \text{all other pairs}\\\text{scale}_{\text{VDW}} = \frac{1}{2.0} \quad \text{1-4 pairs}\\\text{scale}_{\text{VDW}} = 1.0 \quad \text{all other pairs}\\f_{ij} = 0 \quad \text{1-2 pairs, 1-3 pairs}\\f_{ij} = 1.0 \quad \text{all other pairs}\end{aligned}\end{align} \]

where the permittivity constant, \(\varepsilon\), is 1 since the Molecular Mechanics term is in vacuum, and the scaling factors for the electrostatic and VDW terms in AMBER are documented in [cornell:1995]. However, OPLS uses:

\[ \begin{align}\begin{aligned}E_{\text{nb,OPLS}} = \sum_{i < j}[ \frac{q_i q_j}{r_{ij}} + 4\epsilon_{ij}(\frac{\sigma_{ij}^{12}}{r_{ij}^{12}} - \frac{\sigma_{ij}^{6}}{r_{ij}^{6}})]f_{ij}\\\sigma_{ij} = \sqrt {\sigma_{ii} \sigma_{jj}} \quad (\text{collision diameter})\\\epsilon_{ij} = \sqrt {\epsilon_{ii} \epsilon_{jj}} \quad (\text{well depth energy})\\f_{ij} = 0 \quad \text{1-2 pairs, 1-3 pairs}\\f_{ij} = 0.5 \quad \text{(OPLS)} \text{1-4 pairs}\\f_{ij} = 1.0 \quad \text{all other pairs}\end{aligned}\end{align} \]

In PELE, the following formula is used:

\[ \begin{align}\begin{aligned}E_{\text{nb}} = \sum_{i < j} \left[ \text{scale}_{\text{coulomb}}\frac{q_i^{s} q_j^{s}}{\varepsilon_{\text{in}(ij)}r_{ij}} + \epsilon_{i}^{s}\epsilon_{j}^{s} \left( \left( \frac{\sigma_{ij}^{s}}{r_{ij}} \right)^{12} - \left( \frac{\sigma_{ij}^{s}}{r_{ij}}\right)^{6} \right) \right] f_{ij}\\f_{ij} = 0 \quad \text{1-2 pairs, 1-3 pairs}\\f_{ij} = 0.5 \quad \text{1-4 pairs}\\f_{ij} = 1.0 \quad \text{all other pairs}\\q_{i}^{s} = q_{i} \sqrt{\frac{1}{4 \pi \varepsilon_{0}} N_{A} e^{2} \frac{10^{10}}{4.184 \times 1000}}\end{aligned}\end{align} \]

Notice that special charges \(q_{i}^{s}\) are obtained from the actual charges \(q_{i}\) by multiplying by a factor (a squared root, so that when multiplying one charge by the other, we obtain the factor that converts the charge product into the appropriate units) that includes the common Coulomb constant (\(\frac{1}{4 \pi \varepsilon_{0}}\)) that appears in the Coulomnb formulation of charge interaction, as well as the Avogradro number (\(N_{A}\)), unit charge in Coulombs (\(e\)) and other factors that convert from a radius in \(\AA{}\) and a charge in atomic units to kcal/mol.

The parameters in the previous formula resolve differently for the different force fields; to use the same formula for different forcefields, special symbols \(\epsilon_{i}^{s}\) and \(\sigma_{ij}^{s}\) are used. Also, \(\text{scale}_{\text{coulomb}}\) is \(\frac{2}{1.2}\) for 1-4 pairs in AMBER, and 1 in all other cases, including all cases in OPLS. As for the VDW parameters, for OPLS:

\[ \begin{align}\begin{aligned}\epsilon_{i}^{s} = 2 \sqrt {\epsilon_{ii}}\\(\sigma_{ij}^{s})^{2} = 4 \frac{\sigma_{ii}}{2} \frac{\sigma_{jj}}{2}\end{aligned}\end{align} \]

Notice that, in the OPLS template files, the values provided are \(\epsilon_{ii}\) and \(\sigma_{ii}\).

For AMBER:

\[ \begin{align}\begin{aligned}\epsilon_{i}^{s} = 2 \sqrt {\epsilon_{ii}}\\(\sigma_{ij}^{s})^{2} = \left[ \frac{r_{i} + r_{j}}{\sqrt[6]{2}} \right]^{2}\end{aligned}\end{align} \]

Notice that, in the AMBER template files, the values provided are \(\epsilon_{ii}\) and \(r_{i}\) (also named \(R^{*}\)). Also note that \(r_{0} = r_{ij} = r_{i} + r_{j} = \sqrt[6]{2} \sigma_{ij}\) (see http://ambermd.org/Questions/vdwequation.pdf). Remember that \(\sigma_{ii}\) is the intermolecular distance where the potential of two particles of type \(i\) is zero, while \(r_{0}\) is the intermolecular distance where the potential is at a minimum (of value \(\epsilon_{ij}\)).

Note

The VDGBNP solvent model affects the \(E_{\text{nb}}\) term through \(\varepsilon_{\text{in}(ij)}\) (see [zhu:2007]). For all other models (including vacuum), this internal permittivity constant is 1.

The solvation energy

PELE allows working with implicit solvent (otherwise, it works in vacuum), in particular the Generalized Born approximation (see [bashford:2000]). The configuration is described at Solvent Block in Pele++ control file.

When in vacuum, the \(\Delta G_{\text{solv}}\) term is zero. With implicit solvent, the user can select among the following variants:

• Surface Generalized Born + Non-Polar (SGBNP): [ghosh:1998]) gives the detail on how to calculate the Born radii (though the empirical corrections are not followed), and [gallicchio:2002] describes the non-polar term. It uses a penalty term for ion pairs.

• Variable-Dielectric Generalized Born + Non-Polar (VDGBNP): [zhu:2007]; uses the same non polar term as SGBNP.

• OBC: [onufriev:2004] (based on the HCT model of [hawkins:1996]. The actual formulas for calculating the original HCT alphas in the model appears in [hawkins:1995]). It uses the non polar term of [schaefer:1998] as implemented in Tinker.

The polar energy term

All implicit solvent models in PELE use the same formula for the \(\Delta G_{\text{solv,pol}}\) term (see, for instance, eq. 2 in [onufriev:2004], eq. 2 in [zhu:2007], and eq. 19 for [ghosh:1998]), varying in how the dielectic permittivity constant is calculated. Notice that the formula is sometimes called the Still energy from [still:1990].

\[\Delta G_{\text{solv,pol}} = -\sum_{i<j} {(\frac{1}{\varepsilon_{\text{in}(ij)}} - \frac{e^{-\kappa f_{\text{GB}}}}{\varepsilon_{\text{solv}}}) \frac{q_{i}^{s}q_{j}^{s}}{f_{\text{GB}}}} - \frac{1}{2} \sum_{i} {(\frac{1}{\varepsilon_{\text{in}(ii)}} - \frac{e^{-\kappa \alpha_{i}}}{\varepsilon_{\text{solv}}}) \frac{(q_{i}^{s})^2}{\alpha_{i}}}\]

where \(\varepsilon_{\text{in}(ij)}\) depends on the pair of atoms for VDGBNP and is 1 in the other solvent models (as explained in [zhu:2007], a value of 1 is used on the theory that optical polarization of the protein is incorporated into the force field); \(\varepsilon_{\text{solv}}\) is \(80.0\); the second sumation in the formula is what is called self energy in PELE. Notice that the first summation does not have the factor \(\frac{1}{2}\) that usually appears in articles; this is because in this formula each pair of atoms appears only once. Regarding the treatment of ions in solution, if not using the Debye approximation then \(\kappa\) is zero, otherwise, it depends on the ionic strength as follows:

\[\kappa = 0.73 \times 10^{-10} \sqrt{\frac{2N_{A}e^{2}I}{\varepsilon_{r}\varepsilon_{0}RT}}\]

with \(N_{A}\) the avogadro number, \(e\) the elementary charge of a proton, \(I\) the ionic strength of the solution, \(\varepsilon_{r}\) is the relative permittivity of the solvent (\(80.0\)), and \(\varepsilon_{0}\) is the vacuum permittivity. The factor \(10^{-10}\) is to pass the units to \(\AA{}\). The factor \(0.73\) is a correction introduced in [srinivasan:1999a].

The distance function \(f_{\text{GB}}\) depends on the Born radii of atoms \(i\) and \(j\) (also called alphas) in the following way (all models in PELE follow the Still formulation):

\[f_{\text{GB}} = \sqrt{r_{ij}^{2} + \alpha_{ij}^{2}e^{-D}}\]

with \(\alpha_{ij} = \sqrt{\alpha_{i} \alpha_{j}}\) and \(D=\frac{r_{ij}^2}{(2\alpha_{ij})^2}\).

To calculate Born radii (alphas), the atomic radii used to define the solute-solvent dielectric interface are those described in eq. 6 of [gallicchio:2002], with some modifications in OPLS2005, and stated in column Rad. Non Polar SGB of the NBON section in the IMPACT template files (IMPACT template file format).

Warning

Despite eq. 6 of [gallicchio:2002], it can be noticed in file f14_sgbnp.param that this equation is not always applied, at least in OPLS2005. For example, both atom types CA and CA2 (aromatic C) have a standard OPLS radius of \(1.775\AA{}\), but the radii used to calculate Born radii differ for both types (\(2.002\AA{}\) and \(2.500\AA{}\) respectively, in Schrödinger 2017-1 release), following eq. 6 only in the first case.

The penalty term

This term corrects the tendency to form salt bridges of implicit solvent models (see [yu:2004]). In PELE, this term is only present for the SGB solvent, since it is already included in the VDGBNP model through the internal dielectric constants. The OBC model does not use this penalty term either. The penalties are only calculated for 1-4 or non-bonded interactions.

The formula for the penalty term is:

\begin{cases} w_{i} w_{j} V_{\text{max},ij} & \text{if $\text{scale} < -1$} \\ w_{i} w_{j} V_{\text{max},ij} (0.25 \text{scale}^3 - 0.75 \text{scale} + 0.5) & \text{if $ -1 \le \text{scale} \le 1$} \\ 0 & \text{if $ \text{scale} > 1$} \end{cases}

where

\[ \begin{align}\begin{aligned}w_{i} = \max (0.0, a_{0,ij} + a_{1,ij} \alpha_{i})\\w_{j} = \max (0.0, a_{0,ij} + a_{1,ij} \alpha_{j})\\\text{scale} = \frac{2.0 r_{ij} - (\text{cutoff}^{H}_{ij} + \text{cutoff}^{L}_{ij})}{\text{cutoff}^{H}_{ij} - \text{cutoff}^{L}_{ij}}\end{aligned}\end{align} \]

The variable dielectric

The variable dielectric solvent model modifies the internal dielectric of the protein for specific atoms (see [zhu:2007]). This affects the polar part of the solvation energy, as well as the molecular mechanics interactions through the \(\varepsilon_{\text{in}(ij)}\) parameter, which follows this formula:

\[\varepsilon_{\text{in}(ij)} = \max (\varepsilon_{i}, \varepsilon_{j})\]

Notice that this rule is the same as in [zhu:2007], but in the code a \(\min\) function is used, since the actual formula used is:

\[\frac{1}{\varepsilon_{\text{in}(ij)}} = \min (\frac{1}{\varepsilon_{i}}, \frac{1}{\varepsilon_{j}})\]

Also, notice the specific dielectric values for the atoms match for some atoms, but not for all (not for the whole side chain, and also notice the different values depending on the protonation state of the histidine) with those in the article, since the article provides dielectric constants, while PELE works with inverse dielectric values. The values in PELE are (actually, inverse dielectric values):

• N-terminal N and H atoms: inverse dielectric of 0.25

• C-terminal C, O and OXT atoms: inverse dielectric of 1.0

• Lysine atoms NZ, HZ1, HZ2, HZ3: inverse dielectric of 0.25

• Arginine atoms N2, H3, CZ: inverse dielectric of 0.5

• Histidine (HIS, HIE, HID) atoms CG, ND1, HD1, CD2, HD2, CE1, HE1, NE2, HE2: inverse dielectric of 0.5

• HIP atoms CG, ND1, HD1, CD2, HD2, CE1, HE1, NE2, HE2: inverse dielectric of \(\frac{1}{3}\)

• ASP atoms O2, CG: inverse dielectric of 0.25

• GLU atoms O2, CD: inverse dielectric of \(\frac{1}{3}\)

• All other atoms: inverse dielectric of 1.0

As said above, and following [zhu:2007], a value of 1 is used for the dielectric constant of all other atoms (in this case, it corresponds with an inverse value of 1) on the theory that optical polarization of the protein is incorporated into the force field

The non polar energy term

SGB and VDGBNP models

These solvent models follow the SGB/NP model (eq. 13 of [gallicchio:2002]):

\[\Delta G_{solv,npol} = \sum_{i} {C_{\gamma,i} \text{SASA}_{i} + C_{\alpha,i} S(a/\alpha_{i})}\]

where \(C_{\gamma, i}\) and \(C_{\alpha,i}\) are atomic parameters, \(\text{SASA}_{i}\) is the solvent accessible surface area (in \(\AA{}^{2}\)) of atom \(i\) (calculated using radii stated for eq. 7 in [gallicchio:2002], and listed in column Rad. Non Polar Type of the NBON section in the template files, see IMPACT template file format), \(S()\) is a switching function, \(\alpha_{i}\) is the Born radius of atom \(i\), and the adjustable parameters have values \(a=1.5 \AA{}\), (and \(b = 0.4\) and \(c = 12\), not shown).

OBC model

This solvent model follows the ACE model of [schaefer:1998] (eq. 2) as implemented in Tinker (for example, in Tinker 7.1.2), which uses an exponent of 6 for the term \(\frac{R_{i}}{\alpha_{i}}\), as well as uses per atom \(\sigma\) parameters (see https://simtk.org/forum/message.php?msg_id=6151 commenting how OpenMM implements the Tinker version of ACE non-polar solvation).

The actual formula used in PELE:

\[\Delta G_{solv,npol} = 4\pi \sum_{i} {\sigma_{i} (R_{i} + R_{s})^{2} \left( \frac{R_{i}}{\alpha_{i}} \right)^{6}}\]

where \(R_{i}\) is the atomic radii as parameterized for ACE, \(\sigma_{i}\) is a coefficient, \(R_{s}\) is the solvent probe molecule radius (\(1.4 \AA{}\)), and \(\alpha_{i}\) is the atom Born radius.