advanced protocols in order to refine the model of the organic
molecule if the available parameters are suspected to be inadequate
[51, 52].4.3 Defining the
Intermediate States
The core step in the setup of the calculations is possibly the defini-
tion of the alchemical pathway. As depicted in Fig.1, several inter-
mediate states need to be used to link the bound and unbound
states and recover a precise estimate of the binding free energy.
Since we will be sampling discrete alchemical states along the path,
we need to choose how many intermediate states to have, how they
should be distributed along the alchemical path, and for how long
to simulate them. From Fig.1 it is also possible to gather how there
will be two sets of simulations during the calculations: one set in
which the ligand in solution is simulated (left leg of the cycle in
Fig.1), and one in which the protein–ligand complex is simulated
(right leg of the cycle in Fig.1).
The coupling parameterλdefines the thermodynamic state of
the system along the alchemical pathway. This parameter can take
any value between zero and one and is used to scale coulombic
charges, Lennard–Jones parameters and restraint force constants.
Since these are the three sets of parameters that need to be scaled, it
is convenient to have distinct coupling parameters:λcoul,λvdw,λrestr.
In fact, it is often easier to carry out these three transformation
separately, i.e., using a set of windows where only the charges are
changed, then a set where only the LJ are changed, and then a third
set where only the restraints are modified, as exemplified in Fig.3.
However, as mentioned in Subheading2.3.1,ΔGsolvrestrcan be calcu-
lated analytically so thatλrestrapply only to the protein–ligand
complex simulations on the right-hand side of the cycle in Fig.1.
For convenience, and to avoid confusion, in Fig.3 and throughout
the textλ¼0 always indicates the state where the ligand is unre-
strained and fully coupled (both in solution and the complex), and
λ¼1 always indicate the state where the ligand is restrained and
fully decoupled. We can see in Fig.1 that for the protein–ligand
complex simulations we need to calculate the free energy of cou-
pling the ligand to the environment, i.e., going from the decoupled
and restrained state to the coupled and unrestrained state. None-
theless, it does not matter what is defined asλ¼0 andλ¼1, since
the equilibrium free energy of coupling (state D!F in Fig.1) the
ligand is simply the opposite of the free energy of decoupling it
(state F!D in Fig.1); we just need to make sure the correct signs
are used. When using the EXP estimator, there are also considera-
tions about the forward and reverse calculations, as mentioned in
Subheading2.2.1. However, this is not a concern when using the
more robust BAR and MBAR estimators that consider information
from multiple states at once.
As mentioned, there are two main sets of simulations to be run:
the ligand and complex simulations. Figure 3 shows a simple212 Matteo Aldeghi et al.