180

Simonson

Because the preceding protocols use a ligand with two side chains, they are often referred to as ‘‘dual topology’’ methods. A slightly different approach is shown in Figure 3c. Here, the ligand has a single side chain; the transformation changes the OD2 atom of Asp into the ND2 of Asn, and the Asn hydrogens HD21, HD22 are represented by dummy atoms in Asp. This is known as a ‘‘single-topology’’ approach; i.e., a single side chain is used and the number of dummy atoms is kept to a minimum. This approach imposes the same position on the two side chains, which can be a problem if the two ligands occupy different rotameric states in the complex, for example. The internal parameters of the side chain (bond, angle, and dihedral parameters; atomic van der Waals parameters; and charges) must all be altered during the mutation. A change in the length of a covalent bond, for example, contributes the following term to the free energy derivative [32]:

	∂F	2kh(b λ b0	(λ))	db0		(15)
	∂λ			dλ

where kh is the bond force constant, b the instantaneous bond length, and b0 the target bond length, which is a function of the coupling parameter λ. If the bond is stretched on average [b λ b0(λ)], then the free energy to increase b0 is negative. Changes in covalent bond lengths require careful sampling, as a bond’s vibrational frequency is usually distant from the frequencies of most of the degrees of freedom around it, making equipartition of energy with the bond’s surroundings inefficient. The use of Langevin dynamics can significantly improve the sampling. We return to free energy changes as a function of a conformational coordinate in Section IV; for a detailed discussion of singlevs. dual topology approaches, see Refs. 28 and 33.

D. Other Thermodynamic Functions

The free energy is the most important equilibrium thermodynamic function, but other quantities such as the enthalpy and entropy are also of great interest. Thermodynamic integration and perturbation formulas can be derived for them as well. For example, the derivative of the entropy can be written [24]

	∂S		1 ∂U					∂U
T						U λ			λ U λ	(16)
	∂λ		kT		∂λ			∂λ

In the (N, p, T) ensemble, U is replaced by U pV in the right-hand terms. For the energy E E F TS, one obviously has

	∂					∂F			∂S
			E					T		(17)
	∂λ					∂λ			∂λ
More directly,
		B					A UB B UA A			(18)
	E				E

(A and B have the same number of degrees of freedom, so the mean kinetic energies cancel). In the (N, p, T) ensemble, the enthalpy change is

HB HA UB B UA A p(VB B VA A)

(19)

Free Energy Calculations

181

The second pV term is normally negligible [atmospheric pressure corresponds to an en-

thalpy of 1.5 10	5	˚ 3	)]. Thus, the enthalpy can safely be calculated in practice
		kcal/(mol. A

by integrating Eq. (17). Calculations of S, H, and E are intrinsically less precise than free energy calculations; for example, the direct enthalpy calculation [Eq. (19)] involves subtracting two large energies averaged over two separate simulations. However, with increasing computer power, such calculations are rapidly becoming routine [12].

E. Free Energy Component Analysis

In current force fields, the potential energy U is usually a sum of pairwise non-bonded terms (electrostatic and van der Waals) and bonded terms involving groups of two to four atoms. Thus, the energy can be directly decomposed into contributions from groups of atoms and from different energy terms (electrostatic, van der Waals, bonded). From Eq. (11), the free energy derivative ∂F/∂λ can be decomposed in the same way. This leads to a decomposition of the free energy based on groups of atoms and/or different energy terms [27]. For the Asp → Asn example of Figure 3, we have

U(λ) [Unbi2 (λ) Unbi3 (λ)] Ub12 Ub13 U11 U22 U33	(20)
i

where the indices 1, 2, 3 represent the blocks of the system (defined above) and the sum is over atoms of block 1. Superscripts b and nb refer to bonded and non-bonded energy terms, respectively. Atom i of block 1 contributes a term

∂Fi		∂(Unbi2 Unbi3 )	λ
∂λ		∂λ

to the free energy derivative and a term

Fi 1 ∂Fidλ

0 ∂λ

to the free energy difference. Summing over groups of atoms, one obtains the free energy contribution, or component, associated with each group. Unlike the total free energy, the free energy components are not state functions and thus are path-dependent [32,34]; i.e., they depend on the exact way UA is transformed into UB in a particular computation. Nevertheless, a number of studies have shown that when treated with care they can provide important insights into the microscopic interactions important for binding; see, e.g., Ref. 26 and references therein.

III. STANDARD BINDING FREE ENERGIES

The study of receptor–ligand binding is one of the most important applications of free energy simulations [35]. To study receptor–ligand binding theoretically, one must first partition the conformational space into ‘‘bound’’ and ‘‘unbound’’ states. There is no unique way to do this, but in practical situations there is often a natural choice. Thus, conformations where the ligand is within a well-defined binding pocket could be labeled ‘‘bound.’’ Because there is likely to be an energy barrier at the boundary of such a pocket,

182

Simonson

ligand conformations near the boundary will often have high energies and low statistical weights. Therefore, they will not contribute greatly to thermodynamic properties such as the binding constant, which will consequently be robust with respect to the exact definition of the pocket. In addition, when the binding of two similar ligands to a receptor is being compared, there will be some cancellation of the boundary region contributions of each ligand.

The equilibrium binding constant is

	Kb	ρRL		(21)
		ρRρL

where ρR, ρL, and ρRL are the concentrations (or number densities) of receptor, ligand, and complex and Kb has units of volume. The chemical potential of each species in solution is [18,35]

ρA	ZA
µA kT ln ρ0 kT ln			(22)
	Z0Vρ0

where A RL, R, or L; ρ0 is the standard state concentration, V the volume of the system, ZA the partition function of A in solution, and Z0 the partition function of the solution without A. The condition for equilibrium is

kT ln Kbρ0 ∆F b0	(23)
where ∆F 0b µ0RL µ0R µL0	is the standard binding free energy—the free energy to

bring two single molecules R and L together to form a complex RL when the concentrations of all species are fixed at ρ0.

To relate the standard binding free energy to free energies that can be obtained from simulations, we use

0		ZRLZ0Vρ0				QRLQ0Vρ0
∆F b kT ln				kT ln					(24)
		ZRZL				QRQL
		QRLρ0			QL
kT ln			kT ln
	QRQL0/V				QL0 Q0

where the second equality takes into account a cancellation of the velocity partition functions and QL0 is the configuration integral of the ligand alone (i.e., in the gas phase). The second term in Eq. (24) is the free energy to ‘‘annihilate’’ L in solution, i.e., the free energy to reversibly turn off its interactions with the surrounding solution, effectively transferring it to the gas phase. This operation can be done with the ligand either free to move or positionally constrained (fixed center of mass); the free energy is the same, by translational invariance of the solution. The first term is the free energy to ‘‘annihilate’’ the ligand in the binding site, with its center of mass fixed [29]. The standard concentration ρ0 appears explicitly here. Because the binding site does not have translational symmetry, this free energy takes the form of an average over all positions in the active site (see Ref. 29). To compute it, the ligand could in principle be fixed and annihilated in each of these positions, but this is impractical. A much better scheme does the annihilation in two steps [29,36]. First, the ligand’s interactions with its surroundings are reversibly turned off, and

Free Energy Calculations

183

at the same time a harmonic restraining potential is turned on, which confines the ligand to a region centered on the binding site and roughly equal to or slightly larger than it. This free energy does not depend on the standard concentration. Second, the free energy difference between the restrained and fixed ligands is calculated analytically; it has the form kT ln ρ0(2πkT/kh)3/2, where kh is the force constant for the harmonic restraint, and the standard concentration appears explicitly. This scheme is illustrated in Figure 4. For a moderately large ligand (i.e., significantly larger than a water molecule), it becomes necessary to also introduce rotational restraints; their contribution to the free energy is obtained analytically in close analogy to the positional restraint [30]. The use of translational and rotational restraints makes it possible to estimate the ‘‘cratic’’ contribution to the binding free energy [30]. For large ligands, the endpoints of the above ‘‘annihilation’’ processes require very careful sampling and a proper extrapolation of the free energy derivatives.

Figure 4 Thermodynamic pathway for the calculation of a standard binding free energy. In the first step, the interactions of the ligand (L) with its environment (receptor R solvent W) are gradually turned off; at the same time, one or more harmonic restraints are turned on, restricting the translation and possibly the rotation of the ligand. The restrained, ghost ligand (at right in figure) is strictly equivalent to the restrained ligand in the gas phase. The second step removes the harmonic restraints; the corresponding free energy has a simple analytical form (see text).