SZYBKI Theory¶
MMFF94 Force Field¶
The MMFF potential expression is:
where the seven terms respectively describe bond stretching (b), angle bending (a), stretch-bending (s), out-of-plane bending (o), torsion (t), Van der Waals (v) and electrostatic (c) interactions. Their functional forms are given below.
Bond stretching¶
For a bond b between atoms i and j the stretching potential is:
where \(k_{ij}\) is the force constant, \(\Delta r_{ij}\) is the difference between actual and reference bond lengths, and \(c_s\) is so called “cubic-stretch” constant for which the value is -2 Å\(^{-1}\).
Angle bending¶
The bending potential of a bond angle a made by the bonds between atoms i, j and atoms j, k is given by:
where \(k_{ijk}\) is the force constant, \(\Delta \vartheta_{ijk}\) is the difference between actual and reference bond angles, and \(c_b\) is the so called “cubic-bend” constant for which the value is \(-0.4 rad^{-1}\).
Stretch-bend interaction¶
The coupling between the stretching potential of two bonds forming a bond angle and bending that angle is described by:
where \(k_{ijk}\) and \(k_{kji}\) are force constants which couple stretches of i-j and k-j to the i-j-k bend respectively. \(\Delta r\) and \(\Delta \vartheta\) are defined above.
Out-of-plane bending¶
For a trigonal center j, the potential of displacement for an atom l bonded to atom j out of plane i-j-k is:
where \(k_{ijkl}\) is the force constant and \(\chi_{ijkl}\) is an angle formed by the bond j-l and the plane i-j-k.
Torsion interactions¶
For every four bonded atoms i-j-k-l the torsion interaction is described by the term:
where \(V_1\), \(V_2\) and \(V_3\) are constants depending on atoms i, j, k, l, and \(\Phi\) is the dihedral angle formed by bonds i-j and k-l.
Van der Waals interactions¶
For a pair of atoms i and j separated by three or more bonds, where the distance between them is \(r_{ij}\), MMFF adopts the following Van der Waals potential:
where \(R_{ij}\) and \(\epsilon_{ij}\) are defined as follows:
where \(\alpha_i\) is atomic polarizability of atom i, B is 0.2 or 0.0 if one of the atoms is polar hydrogen, \(N_i\) and \(N_j\) are the Slater-Kirkwood effective numbers of valence electrons, \(G_i\), \(G_j\) and \(A_i\) are scale factors.
Electrostatic interactions¶
The electrostatic interaction between two charged atoms i and j separated by at least three bonds is calculated from the standard Coulombic expression:
where D is the dielectric constant for which the default value is 1, \(q_i\) and \(q_j\) are the MMFF partial charges on atoms i and j, \(r_{ij}\) is the interatomic distance, \(\delta\) is the “electrostatic buffering” constant of 0.05 Å. Scaling factor f is 0.75 for 1,4 interactions, and 1.0 otherwise.
Protein-ligand Amber potential¶
Protein-ligand interaction can be described by the Amber force field instead of MMFF potential:
where the summation is over protein-ligand atom pairs. \(r_{ij}\) is the interatomic distance, \(R_{ij}\) is the VdW distance for a pair of atoms, \(q_i\) and \(q_j\) are the Amber partial charges on atoms i and j, and \(\epsilon_0\) is the vacuum permittivity. VdW parameters \(R_{ij}\) and \(\epsilon_{ij}\) and partial charges are taken from work by Wang and coworkers ([Wang-J-2000]).
SMIRNOFF99FROSST and OPENFF¶
SMIRNOFF (SMIRks Natural Open Force Field) is a kind of small-molecule force field
springing from the Open Force Field Initiative (OFF Initiative or OFFI)
(https://openforcefield.org/). The key feature of this force field is its chemical
representation of parameters by using extended SMARTS patterns. “Extended”
here means that the SMARTS patterns also use atom assignment syntax from SMIRKS. For example,
[*:2]~[*:1]~[#1] specifies three atoms, where the first is assigned to be atom 2,
the second atom is assigned to be atom 1, and the third atom is not assigned. This format
frees the force field from using atom types, instead it uses the direct chemical perception
of the extended SMARTS patterns ([Mobley-2018]).
Compared to traditional atom typing approaches, this representation greatly reduces the
number of parameters needed to specify a general small-molecule force field and makes
each parameter more independent.
While it currently has a traditional “Amber” type functional form as specified below, future versions expect to modify the functional form. The OPENFF Initiative will be releasing periodic revisions of the parameters which will include changes to the chemical specifications (ie the extended SMARTS patterns) as well as the numerical parameters. The revision used by the API points will be given with the results.
Currently SMIRNOFF99FROSST and released in October 2019 OPENFF1.0.0 force fields have 6 terms:
describing respectively bond stretching (b), angle bending (a), out-of-plane displacement (o), torsion (t), Van der Waals (v) and electrostatic (c) interactions. Their functional forms are similar to those in AMBER force field and are given below.
Bond stretching¶
For a bond b the stretching potential has a simple harmonic form:
where \(k_b\) is the force constant, \(r\), and \(r_0\) are actual and reference bond length, respectively.
Angle bending¶
For bond angle a, the angle bending potential form is:
where \(k_a\) is the force constant, \(\alpha\), and \(\alpha_0\) are actual and reference bond angle values, respectively.
Torsion interactions¶
Every 4 bonded atoms defined a torsion whose energy is represented with 1 to 6 terms of a Fourier series:
where \(V_n\), \(d_n\), \(p_n\) and \(\gamma_n\) are barrier height, degeneracy, periodicity and phase for the n component of a torsion whose dihedral angle is \(\theta\)
Out-of-plane bending¶
It is defined as potential of a displacement of the trigonal center atom a bonded to 3 other out-of-plane atoms, b, c and d. Its functional form is similar to a regular torsion with a single term. Sometimes it is also called an improper torsion.
Since one can define 3 improper torsions of moving a trigonal atom a out of plane defined by atoms b, c and d, the energy of such a displacement is calculated as 1/3 of the sum of those torsions.
Van der Waals interactions¶
For a pair of nonbonded atoms i and j, SMIRNOFF force field adopts Lennard-Jones potential:
where \(\epsilon_{ij}\) is well depth, \(R_{ij}\) and \(r_{ij}\) are the equilibrium and actual distance between atoms i and j, respectively.
Electrostatic interactions¶
A simple Coulomb term:
represents Coulomb interactions, where \(q_i\) and \(q_j\) are partial atomic charges, \(\epsilon_0\) is the vacuum permittivity, and \(r_{ij}\) is the distance between atoms i and j.
Entropy evaluation¶
Ligand entropy is evaluated as a sum of configurational entropy (\(S_c\)) and solvation entropy (\(\Delta S_s\)):
Configurational entropy¶
Configurational entropy is calculated as:
where q is the conformation dependent partition function:
Here \(q_t\), \(q_ir\) and \(q_iv\) are the translational, rotational and vibrational partition functions respectively, \(n_c\) is the number of unique conformations in the ensemble. All 3 partition functions are calculated from the classical statistical mechanics expressions which could be found in [McQuarrie-1976]. Vibrational frequencies for each conformation, needed for evaluation of \(q_{ir}\) are derived from diagonalization of a Hessian matrix obtained from Quasi-Newton optimization when convergence is achieved. Eigenvalues \(\lambda_i\) of the mass-weighted Hessian:
are converted into wavenumbers \(\tilde{\nu}_i\) according to:
Solvation entropy¶
Solvation entropy is split into electrostatic and hydrophobic parts:
The electrostatic part of solvation entropy is divided in to the bulk component and tight electrostatic polar solute - water interactions (hydrogen bonds). The bulk contribution is estimated from the temperature dependence of the solvent dielectric constant as:
The second term of the electrostatic solvation entropy is estimated as a constant of 28 J/(mol K).
The hydrophobic term, \(\Delta S_{s,hyd}\), is evaluated as:
where \(\Delta G_{s,hyd}\) consists of 3 components:
describing the free energy of cavitation, solute-solvent van der Waals and inductive terms respectively. The cavity formation term is calculated from Scaled Particle Theory [Pierotti-1976]. Analytical expressions for \(\Delta G_{VdW}\) and \(\Delta G_{Ind}\) terms are also taken from the 1976 Pierotti review.
Protein-bound ligand entropy¶
Configurational entropy of a protein bound ligand is calculated totally as vibrational entropy for 3N modes, assuming that 3 rotational and 3 translational degrees of freedom of a solution ligand are transformed into low-vibrational degrees of freedom for the bound ligand.
Solvation entropy for a ligand in the active site is assumed to be a sum of its fractional value in solution determined by the percentage of the ligand surface exposed to the solvent, \(f\Delta S_s\), and a partial desolvation entropy of the protein active site, \(\Delta S_{des}\)
where f is the fraction of ligand surface exposed to the solvent. It is important to notice that \(S_{protein}\) is not an experimentally measurable value, and that only the difference between \(S_{protein}\) and \(S = S_c + \Delta S_s\) might be compared with experimental binding entropy.