4. Scoring Functions

In order to rank different docking ‘poses’ for a given ligand (or to compare the expected binding strengths of a set of ligands) it is necessary to have the ability to make an estimate of the binding strength that reflects the features of the ligand and protein active site (the types of atoms and their positions) as well as the strength of the interactions between them. Thus, each pose can be assigned a docking score and the equation used to achieve this is known as a scoring function (SF).

Scoring functions have traditionally been split into three main groups: physics-based, empirical and knowledge-based. Recently a new, and potentially very powerful, class of SF has been added to this list that is based on the use of machine learning models. Each of these types of SF have their own strengths and weaknesses and the quality and reliability of the results of any protein-ligand docking study will depend on understanding at least the general principles behind the construction of your SF of choice.

For the purposes of this course it is important to introduce the ideas behind physics-based and empirical SF. Knowledge-based and machine learning-based approaches are beyond the scope of the course but information on these can be easily found in the literature if so desired.

4.1. Physics-based scoring functions

Perhaps the most obvious approach to estimating the strength of the interactions between a ligand and a protein binding site is to make use of the well understood physical processes involved. However, due to the size of the molecular systems involved it is always necessary to make approximations.

For example, at the simplest level it is possible to decompose the binding energy, \(E_{bind}\), by taking only the key contributions to the binding between the two molecules. A sensible place to start would be to look at the noncovalent interactions between the protein and ligand which should be amongst the most important things determining the strength of binding.

This corresponds to taking only the classical force field terms for van der Waals (\(E_{vdW}\)) and electrostatic (\(E_{el}\)) interactions (see Force Fields):

\[E_{bind} = E_{vdW} + E_{el}\]

This approach includes the important enthalpy contributions to \(E_{bind}\) but it fails to account for changes in entropy on binding of protein and ligand. This entropy contribution is difficult to include in an efficient manner and really requires that extended molecular dynamics simulations be carried out which is far too computationally demanding for routine applications in protein-ligand docking studies.

One further correction that can be used to improve the classical force-field prediction of \(E_{bind}\) in a relatively simple manner is to include a term for the change in free energy of solvation when the ligand moves from bulk solvent to the protein-bound state:

\[E_{bind} = E_{vdW} + E_{el} + \Delta G_{solv}\]

The term \(\Delta G_{solv}\) is most commonly obtained from an implicit solvent model where the effect of the solvent environment is treated as a constant background but can in principle also be more accurately modelled using explicit solvent molecules which gives a much more realistic picture of solvent-solute interactions (especially important for solvents such as water where a correct description of hydrogen-bonding is vital).

Unfortunately, even if these extra terms are included scoring functions based on classical force field methods are limited in the accuracy that they can achieve. This comes from the nature of the classical approximation itself where, e.g. polarisation effects, are not usually incorporated in the force field design. In this case, where a contribution to \(E_{bind}\) would be expected from polarisation of both the ligand and the active site upon binding none will be observed.

Perhaps more fundamental is the problem that different force field parameterisations cannot be expected to yield the same energies for a given interaction. Thus, scoring functions based on different force fields also cannot be expected to provide consistent values of \(E_{bind}\).

A work-around exists to take care of these problems: replace the force field energies with those obtained from ab initio electronic structure calculations (note that the \(\Delta G_{solv}\) contribution will still generally come from an implicit solvent model)

\[E_{bind} = E_{QM} + \Delta G_{solv}\]

Whilst the quantum mechanical approach can, in principle, provide much more accurate energies for e.g. polarisation and charge-transfer effects which are important in binding the computational cost is so high for these methods that they are effectively unusable except on large high-performance computing systems. One approach to ease this problem is to combine a quantum mechanics treatment of the important part of the system (ligand and active site residues) with a classical force field molecular mechanics treatment of the rest of the protein.

Despite this, any scheme that employs quantum mechanics will be very much slower than force field based approaches and for this reason it is still uncommon to find examples of the use of electronic structure methods in docking.

The classical energy expression is still by far the most common of the physics-based approaches to scoring function construction and a good example is the one used [1] in the UCSF DOCK software which takes the form:

\[E_{bind} = \sum_{i=1}^{lig} \sum_{j=1}^{rec} \left( \frac{A_{ij}}{r_{ij}^{12}} - \frac{B_{ij}}{r_{ij}^6} + \frac{q_iq_j}{\epsilon r_{ij}} \right)\]

where for atoms i and j

\[ \begin{align}\begin{aligned}\frac{A_{ij}}{r_{ij}^{12}} - \frac{B_{ij}}{r_{ij}^6} = E_{vdW}\\\frac{q_iq_j}{\epsilon r_{ij}} = E_{el}\end{aligned}\end{align} \]


The last term in this equation \(E_{el}\) depends on atomic partial charges, \(q\). This is a potential problem because these charges can be calculated using various different methods and so there is no way that these can all produce correct results when fed into the scoring function.

For this reason care must be taken when using physics-based scoring functions to ensure that the correct procedure is followed in order to obtain meaningful values.

4.2. Empirical scoring functions

As can seen with the physics-based models, there must be a trade-off between accuracy and speed. In order to obtain an accurate description of the interactions that determine the strength of binding between a protein and ligand more (and often very significant) computing time and power is required.

An alternative to simply adding more physically motivated terms to the equation for \(E_{bind}\) is to generalise and scale the existing ones to improve the agreement to a set of known experimental binding data. The resulting empirical scoring functions retain the speed of the simple physics-based models but, if well constructed, can provide much more accurate results.

An example of this approach is the X-Score SF [2]:

\[E_{bind} = w_0 + w_1\Delta G_{vdW} + w_2\Delta G_{Hbond} + w_3\Delta G_{rot} + w_4\Delta G_{hydro}\]

Here \(w_0, w_1, \ldots\) are empirical weights used to scale the various contributions. In the context of this course it should be noted that Autodock Vina uses a functional form that was based on X-Score and is broadly similar with the main difference being that the X-Score SF only contains intermolecular contributions whilst Vina SF also has terms relating to intramolecular interactions [3] and can therefore rank different ligand conformations more accurately.

It is important to remember when choosing/using empirical SF is that their performance reflects the training data that was used in their creation. Fortunately this is usually not a concern for ‘standard’ drug-like ligands but a general empirical SF cannot be expected to perform well for highly unusual ligands (this is similar to the problem of parameterisation in force fields).