Royal Society Publishing

Grid-based steered thermodynamic integration accelerates the calculation of binding free energies

Philip W Fowler, Shantenu Jha, Peter V Coveney


The calculation of binding free energies is important in many condensed matter problems. Although formally exact computational methods have the potential to complement, add to, and even compete with experimental approaches, they are difficult to use and extremely time consuming. We describe a Grid-based approach for the calculation of relative binding free energies, which we call Steered Thermodynamic Integration calculations using Molecular Dynamics (STIMD), and its application to Src homology 2 (SH2) protein cell signalling domains. We show that the time taken to compute free energy differences using thermodynamic integration can be significantly reduced: potentially from weeks or months to days of wall-clock time. To be able to perform such accelerated calculations requires the ability to both run concurrently and control in realtime several parallel simulations on a computational Grid. We describe how the RealityGrid computational steering system, in conjunction with a scalable classical MD code, can be used to dramatically reduce the time to achieve a result. This is necessary to improve the adoption of this technique and further allows more detailed investigations into the accuracy and precision of thermodynamic integration. Initial results for the Src SH2 system are presented and compared to a reported experimental value. Finally, we discuss the significance of our approach.


1. Introduction

The ability to rapidly and accurately calculate free energies of binding is of widespread importance in condensed matter research. For example, understanding the interaction between a drug candidate and its target, usually a protein, is essential in both molecular biology and the pharmaceutical industry. The Gibbs free energy of binding (ΔGbind) characterizes the affinity of one molecular species for another and is a vital step towards understanding any biomolecular interaction.

However, in spite of recent improvements in the theory of and techniques for computing such free energies (Simonson et al. 2002), much remains to be done. The adoption of formally exact techniques for computing free energies, such as thermodynamic integration (Kirkwood 1935; Leach 2001), has been limited by the length of time required to reach a solution (Chipot & Pearlman 2002). This is a result of the substantial quantities of computing resources necessary to achieve adequate convergence. The aim of this paper is to describe a novel approach using a computational Grid that addresses this bottleneck. This will enable the increased adoption of thermodynamic integration and also permit more detailed investigations into the accuracy and precision of this technique. Although we discuss the calculation of a difference in the free energy of binding in the context of signal transduction (see §2), the methods described herein are equally applicable to other biological or physical condensed matter systems.

That the demand by scientists for computational resources will continue to outpace supply is universally acknowledged. Satisfying these demands is made harder by the need of many computational scientists for increasing amounts of tightly coupled, high-performance capability computing resources with a high quality of service (Harting et al. 2003). To meet this insatiable need requires creative solutions that go beyond Moore's law. Grid computing—distributed computing performed transparently across multiple administrative domains—is one possible solution to this problem (Foster & Kessleman 1999). This is reminiscent of high-performance distributed computing with the important additional feature that multiple administrative domains are transcended seamlessly, creating a single ‘meta-computer’ and thereby dramatically reducing the time to achieve solutions.

The effectiveness and impact of a computational solution to a problem is greatly enhanced by reducing the time taken to achieve a result to the same as or less than that achieved by experiment. We show in §3 that the combined use of computational Grids and the RealityGrid1 computational steering system (Chin et al. 2003; Fowler et al. 2004) is now making this possible. Details of how the model was constructed, the data analysed and the sources of error encountered are stated in §4. Our results are presented in §5 and they are discussed, along with plans for future work, and the implications of the Grid-based approach to computing differences in binding free energies, in §6.

2. Src homology 2 domain: peptide binding

A common form of intracellular signalling pathway incorporates protein tyrosine kinases. These are enzymes which phosphorylate tyrosine residues of specific intracellular proteins to indicate the presence of a chemical messenger bound to a receptor outside the cell. Other proteins containing Src homology 2 (SH2) domains recognize and bind these phosphorylated tyrosines in a specific manner and initiate a cascade of further protein–protein interactions, collectively referred to as tyrosine kinase-mediated signalling (TKM), leading ultimately to the desired effect, typically gene regulation (Shakespeare 2001).

The elucidation of the crystal structure of the Src SH2 domain by Waksman et al. (1992) initiated a significant research effort directed towards understanding how these protein domains bind and recognize specific peptide sequences. Inhibition of specific SH2 domains is expected to lead to the control of specific TKM pathways through the activation or deactivation of the genes that they regulate. Precisely because these pathways are so generic, the possibilities to influence a wide range of ailments, from cancer to immunological disorders, are large (Waksman et al. 2004). Despite this effort, there has been limited success in developing drug candidates (Shakespeare 2001).

It has proved extremely difficult to discover molecules that specifically inhibit a SH2 domain (Ladbury & Arold 2000), making a detailed thermodynamic study of SH2-mediated binding and recognition both interesting and useful. NMR solution studies have shown that changes in chemical shift upon binding the peptide extend over large regions of the Src SH2 domain (Ladbury & Williams 2004). This suggests that the binding has a large dynamic component and may explain why simple empirically based computational models, while fast, are insufficiently accurate for SH2 domains (Henriques et al. 2000; Ladbury & Williams 2004). Formally exact computational methods, such as thermodynamic integration, include the dynamics of the system as well as the effect of solvent molecules and ions, and therefore, are promising methods of investigation.

Experimental techniques, for example isothermal titration calorimetry (ITC; Lubman & Waksman 2003), while accurate, can only measure macroscopic thermodynamic properties including, in addition to the Gibbs free energy of binding (ΔGbind), the enthalpy (ΔHbind) and entropy (ΔSbind) of binding. It is possible, through site-directed mutagenesis studies, to examine the contribution to the free energy that different residues make, but this is a tedious process and is limited in the insight that it can derive (it usually assumes a simplistic additive effect from each residue). In addition to calculating the free energy, thermodynamic integration (and other related computational techniques) has the potential to relate the strength of binding to microscopic, molecular behaviour: for example, how the binding of a peptide affects the extent of the solvation shells (Simonson et al. 2002). Detailed computational studies and experimental measurements complement one another with experiment informing computation and vice versa. Unfortunately, thermodynamic integration is computationally very expensive and, for the type of system described, has a far longer ‘turnaround’ time than experiment. We shall discuss this in more detail shortly, but first we shall outline briefly the theory on which thermodynamic integration is based to provide the context to describe our Grid-based solution.

Thermodynamic integration is an established computational method for calculating free energy differences (Chipot & Pearlman 2002) within both the academic (Hansson et al. 2002) and pharmaceutical (Jorgensen 2004) communities. To illustrate the basis of the technique, consider a thermodynamic cycle as shown in figure 1. The difference in free energy between molecules A and B binding to a given SH2 protein domain is given byEmbedded Image(2.1)where we have made use of the fact that ΔG is a thermodynamic state function. Next, we assume that ΔG is a continuous function of a parameter, λEmbedded Image(2.2)where λ represents in some sense the extent to which the system has changed from A to B. At intermediate values of λ, the system is in an unphysical state. The concept of changing one molecule into another leads to this procedure being referred to as an alchemical transmutation. Fortunately, computing ΔG1 and ΔG2, while not trivial, is at least computationally tractable. Leach (2001) contains a good description of the theory of thermodynamic integration which we will only outline here. We first substitute in the standard thermodynamic resultEmbedded Image(2.3)where Q(λ) is the partition function for our system. Further writing for Q(λ) in terms of the Hamiltonian and thereby the potential energy, U, gives us the final resultEmbedded Image(2.4)The angular brackets represent the ensemble averaging of the free energy derivative. This integral is numerically evaluated over a set of intermediate values, or ‘windows’, of λ. When selecting the set of λ-values, care has to be taken to balance minimizing the number of simulations against accurately mapping the free energy derivative function. Each simulation must be run for a sufficient duration to ensure first, that it is equilibrated and second, that the ensemble average within equation (2.4) has properly converged (see §5 for more detail).

Figure 1

The thermodynamic cycle used to compute the difference in binding free energy between peptide A and the SH2 domain and peptide B and the SH2 domain. A plus indicates the species are bound and a minus that they are unbound.

The amount of computational resource needed to compute a single thermodynamic integration is usually very large, although this varies depending upon the number of atoms involved in the alchemical mutation and the size of both the bound and unbound systems. Even if each simulation is run on a tightly coupled parallel architecture then it may be weeks to months before a result is obtained: for examples of complicated studies see the articles by Archontis et al. (1998), Michielin & Karplus (2002) and Yang et al. (2004). By contrast, a single ITC experiment, assuming one has prepared protein and peptide solutions, only takes a day or two to measure a binding free energy. This disparity in turnaround time, along with the cost and difficulty of performing the computations, has prevented the wide-spread adoption of thermodynamic integration within academia and the pharmaceutical industry (Chipot & Pearlman 2002). We describe in §3 a Grid-based approach which dramatically speeds-up such calculations and makes them easier to perform.

3. Thermodynamic integration on the Grid: motivation, design and implementation

In this section, we discuss the need for a Grid-based solution to compute free energies and describe a Grid workflow that reduces the time to solution. The implementation of the Grid workflow using the RealityGrid computational steering system (ReG CSS; Chin et al. 2003; Pickles et al. in press) is described before we discuss our experiences of developing this approach and demonstrating the calculation of a single difference in binding free energy. This Grid-based approach to computing free energies using thermodynamic integration provides an example of the use of both Grid computing and the ReG CSS. For a further example concerned with defect dynamics in gyroid liquid crystals see González-Segredo & Coveney (2004) and Pickles et al. (2004).

The solution must meet the following requirements. It must be able to launch and handle multiple, long-running, tightly coupled parallel simulations, each of which may behave differently. For example, simulations differing only in the value of a single parameter may take different lengths of time to converge and therefore consume varying amounts of computational time. The scientist must also be able to interact with and steer each simulation to, for example, ensure the effective use of resources.

(a) Motivation

Although in principle running many simulations concurrently in different administrative domains is possible without the use of Grid technology, in practice this would require an enormous effort to cope with the heterogeneity of the different computers used. The aim of Grid middleware is to shield the user from these complexities and leave the scientist free to interact with the simulations as if they were running locally. The current versions of Grid middleware are some way from achieving this vision and the Grid computing infrastructure itself is not yet in a robust and reliable state (Chin & Coveney 2004). Despite this, the potential offered by Grid computing is large.

There are many theoretical approaches to computing free energies (Chipot & Pearlman 2002). Thermodynamic integration is especially amenable to a Grid-based approach as it naturally partitions the problem into several independent simulations, each requiring tightly coupled parallel resources. As each simulation can be distributed around the Grid, it is easily possible to be using at any given instant in time, more high-performance computing (HPC) resources than are available in a single administrative domain. As a result the wall-clock time taken to reach a solution can be dramatically reduced.

The peptide-SH2 protein system is not as large as some of the problems we and others are currently studying (Wan et al. 2005). In addition to the scientific interest inherent in the binding of peptides to SH2 domains, its intermediate size makes it a good candidate to demonstrate proof-of-principle for the acceleration of exact calculations of binding affinities using computational Grids.

(b) Grid workflow

Consider the process of calculating a difference in binding free energy using thermodynamic integration on a Grid from the perspective of a scientist. First, a single initial simulation, usually at a low value of λ, is launched on a HPC resource somewhere on the Grid. The scientist monitors the simulation and assesses when to spawn a second simulation with the next value of λ. Although there is no obvious parameter to directly monitor to assess when the simulation has converged, it is prudent to use, at a minimum, the convergence of the running average of the free energy derivative. When the scientist decides to spawn a new simulation, a checkpoint is taken and the new simulation with the next value of λ is started on another, or the same, HPC resource. The original parent simulation continues for a suitable duration accumulating data. Meanwhile, the newly spawned simulation is monitored in exactly the same manner as just described to assess when the third λ simulation can be spawned. This process of spawning, or cloning, is carried out as often as available resources permit and ensures that the equilibration time for each λ simulation is minimized. This is an important consideration when studying alchemical mutations involving large numbers of atoms. The scientist monitors and controls each of the simulations, the number of which is constrained by the resources available within the Grid. For example, at any time a single λ simulation may be migrated onto a different computer in a separate administrative domain. This is a good example of the flexibility of this approach and is necessary if, for example, a HPC resource has scheduled maintenance or the scientist wishes to take advantage of a faster or cheaper machine. He or she continues the process until all simulations have converged and the desired amount of data have been collected. Finally, the data are analysed and the value of ΔΔGAB is calculated. By comparison, a workflow for typical thermodynamic integration performed using molecular dynamics (MD) is serial in nature (Donnini & Juffer 2003; Wan et al. 2004; Yang et al. 2004); each λ simulation is run to completion before the next is launched, often on the same machine. The acceleration over the conventional approach stems from the twin abilities of computational Grids to supply increased resources and to allow the easy interaction with multiple simulations, wherever they are running.

(c) Implementation

The ReG CSS is central to our ability to perform the Grid workflow just described and provides the functionality to carry out the necessary steps in a generic, extensible and modular way. Interfacing the ReG CSS with an existing classical MD code automatically implemented the functionality required to perform a thermodynamic integration calculation on a Grid. Onus is placed on the application developer, often the scientist, to interface the application code with the ReG CSS. We chose to interface the ReG CSS with the classical MD package NAMD2.5 (Kalé et al. 1999). This proved to take longer than expected, and we will discuss the reasons for this in §3d. We shall not discuss the design, architecture or scope of the ReG CSS as these are described elsewhere in this issue (Pickles et al. in press). The ReG CSS provides the scientist with the ability to launch, checkpoint, restart and migrate instances of a simulation code on a Grid and to monitor and steer individual variables within each simulation. It is often referred to as ‘upper-level’ middleware, to distinguish it from ‘lower-level’ middleware like Globus2 or Unicore.3

The scientist launches and migrates individual simulations using the supplied launcher application. A separate application, the steering client, allows the scientist to monitor and control all the simulations he or she has running. This client usually takes the form of a conventional application running on a desktop or laptop computer however, to provide flexibility, there exist steering clients that can run on a personal digital assistant (Kalawsky et al. in press) or can be accessed through a web portal. The steering client uses the RealityGrid steering Grid service via a registry to communicate with all the active simulations.

(d) Experience

It was a non-trivial task to interface the ReG CSS with NAMD, primarily due to difficulties with using globus sockets to handle the communication between NAMD and VMD (Humphrey et al. 1996), its attendant visualization package. Many of the difficulties were finally solved by using TCP/IP sockets for this communication. We have interfaced other codes with the ReG CSS in less time, for example a lattice-Boltzmann code (Chin et al. 2003; Harting et al. 2004) and therefore in our experience, this was an unusually difficult application to interface with the ReG CSS.

We demonstrated this Grid-based method for computing differences in binding free energies at the UK e-Science All Hands Meeting 2004 in Nottingham, UK, using the workflow previously described (Fowler et al. 2004). Figure 2 shows how the UK National Grid Service (NGS)4 and the US TeraGrid5 were federated using the US–UK network link, UKLight.6 This combined Grid was used to run 24 λ simulations of varying durations and sizes in 41 h, simulating on average 0.3 ns for each λ value. Halfway through the calculation there were nine simulations running on five different machines in both the US and UK using a total of 140 processors. By contrast, a practice calculation of the same degree of rigour took four weeks to complete using up to half the capacity of a 36 processor SGI Altix machine at University College London.

Figure 2

The computational Grid and the network used at the UK e-Science AHM Meeting 2004.

Many difficulties were encountered during the course of the experiment. These ranged from the simple (running out of disc quota) via the annoying (the occasional failure of software) to the catastrophic (an entire node of the NGS crashing). As an indication of the effort required to mount this pioneering calculation, up to six people were working on solving various problems at any one time with a further person controlling all the active simulations. Securing acceptance of our UK e-Science certificates by all the HPC resources comprising our Grid was a key obstacle and will likely remain so until projects of this kind become commonplace. The lack of specialized, interactive queues on the HPC resources comprising this Grid meant that it was always a risk that launched simulations would get stuck in a batch queue. Fortunately, due to the relatively small numbers of processors (8–24) required for each λ simulation, we were usually able to avoid this situation. If this Grid-based approach is to be used more widely or on larger systems then HPC queuing policies will need to change to allow instantaneous access, albeit at a higher cost, to the scientist. If these changes cannot be made, routine use of such interactive Grid-based approaches will remain difficult.

4. Methods

In this section, we discuss the motivation for investigating the difference in binding affinities between our two chosen peptides before we list the details of the construction of the model and the MD simulations. Since the analysis of the data is complicated, we also briefly describe the method to compute the difference in binding free energies.

(a) The model

We chose to compute the difference in the free energy of binding between peptides of sequence PQpYEpYIPI (pYEpYI) and PQpYIpYVPI (pYIpYV, single amino acid notation; pY, phosphorylated tyrosine; Alberts et al. 1994) and the Src SH2 domain (see figure 3). The ITC experimental measurements of the free energy of binding for both peptides by Lubman & Waksman (2003) show that both peptides have similar free energies of binding (ΔΔG=−1.0±0.2 kcal mol−1), one binding predominantly entropically, the other predominantly enthalpically. This change in behaviour results from the substitution of two amino acids, including the introduction of electric charge into the system through the glutamic acid to isoleucine mutation. As we shall see, the large number of atoms being alchemically transmuted and the introduction of charge makes the precise determination of the free energy of binding difficult. In ongoing work, we are studying the binding of simpler peptides to the Src SH2 domain.

Figure 3

On the left an Src SH2 domain is shown bound to a peptide. As this is shown in the ‘cartoon’ representation, the peptide is only visible as a blue tube running from the left α-helix to the right α-helix. On the right both peptides are shown without any hydrogens. This figure was produced with VMD (Humphrey et al. 1996).

The structures of both peptides bound to the SH2 domain were also solved by Lubman & Waksman (2003; PDB codes: 1NZL and 1NZV). The structure of pYEpYI bound to the SH2 domain (1NZL) was taken to be the λ=0 state and an alchemical peptide with dual topology was constructed. Both unbound and bound systems were solvated with TIP3P water resulting in system sizes of 5025 and 18 651 atoms, respectively. Each system was subjected to a combined minimization, warming and equilibration protocol. PME was used to calculate the electrostatic forces and periodic boundary conditions were employed (Darden et al. 1993). The van der Waals forces were cut off at 12 Å with a switching function applied from 10 Å. SHAKE was used to facilitate a 2 fs timestep (Ryckaert et al. 1977) and all simulations were run in an NPT ensemble using a Langevin thermostat and Berendsen barostat. To match experimental conditions, the temperature was fixed at 298 K and the pressure at 1 atm.

Existing CHARMM27 parameters for phosphotyrosine were used (Feng et al. 1996; MacKerell et al. 1998) in conjunction with the NAMD2.5 parallel MD package (Kalé et al. 1999). The thermodynamic integration functionality of NAMD2.5 was employed using the dual topology formulation and a linear scaling of the two Hamiltonian endpoints (Dixit & Chipot 2001). The spacing and number (11) of λ-values were chosen to minimize the number of simulations required while mapping as accurately as possible the free energy derivative, which itself was written to disc every 40 fs.

(b) Analysis of data

Configurations generated by a MD simulation are correlated, thus successive measurements of the free energy derivative are not independent. We used the method of statistical inefficiency (Leach 2001; Yang et al. 2004) to estimate the correlation time (τ) for each λ simulation. Values ranging from 25 to 75 ps were observed and we chose a conservative value of 100 ps for all λ simulations. This is larger than the correlation time observed by Yang et al. (2004), however, this is to be expected since our alchemical transformation has a larger number of degrees of freedom. To verify that the correlation time was correctly estimated, ΔΔG was calculated with correlation times of 200 and 400 ps. No significant variation in results was found, save the expected increase in the error as the number of independent measurements decreased.

The values of the free energy derivative were then calculated from these N independent measurements for each λ simulation. Figure 4 plots Embedded Image (see equation (2.4)) for all values of λ for both the bound and unbound systems. The individual alchemical free energy changes, ΔG1 and ΔG2, for the bound and unbound systems, respectively, were computed by integrating the resultant function according to equation (2.4).

Figure 4

The integrand from equation (2.4) for both the bound and unbound systems (6 ns duration of data included—see table 1). Note error bars are smaller than the size of the points.

View this table:
Table 1

The variation of the computed values for ΔΔGAB with the duration of data included in the calculation (using extrapolation method 3 of §4b)

Figure 4 clearly illustrates the well-known ‘endpoint’ problem: as λ=0 and 1 are approached, the van der Waals forces exerted by the ‘vanishing’ alchemical atoms on the remainder of the system are scaled by λ. This introduces a discontinuity by preventing the surrounding solvent molecules sampling the space excluded by the van der Waals repulsive term until an endpoint, λ=0 or 1, is reached. Using a single rather than dual topology description might ameliorate this problem, but in our opinion the mutation is too complex for it to help here. A further approach (which we have adopted) is to dynamically modify the width of the windows as the endpoints are approached. Non-linear scaling of the Hamiltonians can reduce the problem but such an approach is limited in scope. The use of a ‘soft-core’ van der Waals potential (Zacharias et al. 1994), i.e. one that smoothly reduces to zero with λ, would be expected to solve this problem, but its implementation in a package like NAMD2.5 is not simple. Instead, we must extrapolate the values of the free energy derivative at the endpoints.

Extrapolation introduces a significant source of error into our calculation of the difference in binding free energy. We use three separate extrapolation methods of increasing sophistication so that we may understand the magnitude of this error. The first method (method 1) uses no extrapolation: we simply integrate the function for 0.02≤λ≤0.98. This is equivalent to assuming the function is zero for 0.02≤λ≤0 and 0.98≤λ≤1. The second method (method 2) assumes that the values of the function at λ={0, 1} can be approximated by those at λ={0.02, 0.98}. Any differences in the gradients at the λ endpoints are therefore ignored. The final method (method 3) splits the function into its electrostatic and van der Waals components. For the electrostatic component, the values of the function at λ={0, 1} are linearly extrapolated using the last three datapoints, weighted by their standard errors. The values of the van der Waals component at λ={0, 1} are extrapolated in a similar way by fitting a function of the form −3/4+b, where a and b are fitting coefficients. This functional form has been shown to be analytically correct by Simonson (1993) and is an established technique (Archontis et al. 1998).

The alchemical free energy derivatives are integrated using the trapezium rule to compute the alchemical free energy differences (ΔGA and ΔGB). N-values of ΔGA and ΔGB are separately calculated using single sets of N independent measurements of the free energy derivatives. Each 100 ps of data contributes an extra measurement to the set of alchemical free energy changes. The standard errors for ΔGA and ΔGB are computed from these N-values.

(c) Sources of error

When computing a difference in free energy using thermodynamic integration, it is essential to minimize and quantify all the errors present. In addition to the extrapolation error mentioned above, the two main sources of error are the imprecise description of the Hamiltonian and the incomplete exploration of phase space at each value of λ. The latter arises, for example, due to simulations being of insufficient length to achieve convergence of the ensemble average (Simonson et al. 2002).

To both mitigate and investigate the error resulting from unconverged simulations, we massively extended all the λ simulations performed during the course of the UK e-Science All Hands Meeting 2004 by at least an order of magnitude in duration. Each λ simulation was extended for a minimum of 7 ns, with some simulations being extended by up to 25 ns, making a total of 0.235 μs of MD trajectory, split evenly between the bound and unbound systems. The protein was stable for all simulations, for example, the average root mean square deviation between the alpha carbons of the X-ray crystallographic structure and the simulation trajectory for the bound λ=0.92 simulation was 1.42 Å with a standard deviation of 0.17 Å.

We emphasize that this is extremely long for any thermodynamic integration calculation and would not have been possible without the extensive use of HPC. Using the TeraGrid node at NCSA, a state-of-the-art supercomputer, still required 31 500 CPU hours. This permits the investigation and analysis of the behaviour of the binding free energy as a function of the amount of data used in the calculation.

It remains difficult to assess whether a λ thermodynamic integration simulation has equilibrated since the free energy derivative is very small compared to the fluctuations in the total system energy. In a promising method recently introduced by Yang et al. (2004), the border between the equilibrating and equilibrated region is determined by applying a statistical test for normality as the dataset is penetrated in reverse order. It assumes that the running average is stationary and normally distributed when the simulation has equilibrated. The first point at which the probability of the data being drawn from a Gaussian falls below a chosen level of significance is taken as the border between the two regions. However, the application of the Wilk–Shapiro test for normality (Shapiro & Wilk 1965; Thode 2002) to our data did not lead to a clear identification of equilibrating regions. Unlike Yang et al. (2004), we used a conventional Lennard–Jones 6–12 van der Waals potential which is acknowledged to slow the equilibration of the system as the scaled alchemical atoms perturb the motions of other atoms (Yang et al. 2004). We cannot yet conclude whether our simulations are not equilibrated or whether the technique is inappropriate for our chosen method. We assume that the λ simulations may not all be equilibrated and therefore we may not have minimized this source of error. This will be investigated in §5 by examining the sensitivity of ΔΔGAB, the amount of data used in its calculation.

5. Results

Computed values for ΔΔGAB are shown in table 1 as a function of the amount of data included in the calculation when extrapolation method 3 is used (see §4 for definition). The first 0.2 ns of each λ simulation was excluded from this calculation to remove any significant initial transients. During the UK e-Science All Hands Meeting 2004, we collected data for approximately 0.2 ns at each λ value, based upon which our initial and unsophisticated analysis suggested that the value for ΔΔGAB was in the region of −9 to −12 kcal mol−1. This is significantly different to the experimental value of −1.0±0.2 kcal mol−1. As more data were included in the calculation, the accuracy relative to the experimental value improved. The change in the computed value of ΔΔGAB, as each additional nanosecond of data was included, decreases until successive values are statistically indistinguishable.

The same analysis was repeated using the other methods for extrapolating values of the free energy derivative at λ={0, 1}. The values of ΔΔGAB are plotted in figure 5a. Improvement in the accuracy and precision of the computed values for ΔΔGAB is observed for all methods of extrapolation, the majority of the improvement occurring within the first 2 ns. The computed value for the difference in binding free energies varies significantly with the chosen method of extrapolation, and therefore our errors likely underestimate the contribution from extrapolating values of the free energy derivative at λ={0, 1}. Method 3, the extrapolation method with the most theoretical grounding and therefore expected to be the most accurate, is bracketed by the other two simpler methods. We might hypothesize that if we have mapped the free energy derivative sufficiently well, then the final result will not be sensitive to the method of extrapolation used. This is not the case and therefore this calculation would benefit from the additional measurements of the free energy derivative at new values of λ. Running the entire simulation in reverse, i.e. from λ=0.98 to 0.02 would give a further estimate of the systematic error and provide another check on the consistency of the approach.

Figure 5

The variation of the computed values for ΔΔGAB with both the extrapolation method used (§4b) and the duration of data included in the calculation. The errors are calculated as described at the end of §4b. (a) Dataset ordered in forward direction, (b) dataset ordered in temporally reversed direction.

By analogy with the approach taken by Yang et al. (2004), we have repeated this analysis on a temporally reversed dataset. This is illuminating as it necessarily excludes the first, and therefore presumably the least converged, part of each dataset. The dataset used is not the same as the forward dataset. The last 7 ns from all the simulations, irrespective of their durations, were taken and hence for the longer simulations, a significant portion of the data is not included. The variation of ΔΔG with the quantity of data can be seen in figure 5b. The more rapid convergence of ΔΔG and the agreement, to within errors, of the value computed using the forward ordered dataset suggests that the majority of the equilibration occurs within approximately the first 2 ns, although this is unlikely to be true for all λs.

The behaviour of the errors is plotted in figure 6 and, perhaps contrary to expectations, the computed errors in the value of ΔΔGAB do not monotonically reduce as more data are included in the calculation. This is a result of the individual λ simulations insufficiently sampling phase space; short λ simulations do not measure the fluctuations in the free energy derivative and consequently underestimate the error. This is an acknowledged weakness of thermodynamic integration, but it is unusual to observe the effect over time-scales, as long as nanoseconds.

Figure 6

The variation of the computed errors in ΔΔGAB for the forward dataset with both the extrapolation method used and the duration of data included in the calculation.

Table 2 decomposes the free energy differences and therefore the value of ΔΔGAB into its electrostatic and van der Waals components. As expected, due to the destruction of the charge, both ΔG1 and ΔG2 are dominated by the electrostatic component and are in excess of +100 kcal mol−1 as we are alchemically transforming a hydrophilic peptide (pYEpYI) to a hydrophobic peptide (pYIpYV). Neither component of the unbound system varies significantly with the quantity of data used; however, the van der Waals component of the bound system varies by 4.8 kcal mol−1, thereby accounting for the majority of the variation in the computed value of ΔΔGAB. This suggests that the bound system is taking longer to equilibrate, which is not surprising as in this case the process is dominated by the relatively slow local rearrangements around the alchemical atoms.

View this table:
Table 2

Variation of the electrostatic (elec) and van der Waals (vdW) components of ΔΔGAB with the duration of data

6. Conclusion

Our Grid-based method dramatically accelerates the calculation of differences of binding free energies and was demonstrated at the UK e-Science All Hands Meeting 2004 (Fowler et al. 2004). The reduction in the time taken to calculate a result will hopefully encourage the adoption of this useful technique within condensed matter research. Computational techniques for computing differences in binding free energies yield more insight into a binding event than that of experiment and, if sufficiently rapid, can complement, or compete with, experiment. We emphasize that, like the technique of thermodynamic integration, our Grid-based approach can be applied to many different fields and is not restricted to protein–ligand interactions. We are currently using a somewhat similar Grid-based approach to study a very different problem, namely, understanding the translocation of DNA through protein nanopores.

The reduced time to compute differences in binding free energies permits detailed investigation of the technique of thermodynamic integration itself. For example, the ability to collect more data allowed us to examine the variation of the computed value of ΔΔGAB with the amount of data used in its calculation and, therefore, we can assess if the simulations have equilibrated and are converged. Using large quantities of HPC, we massively extended all the individual λ simulations resulting in a total of 0.235 μs of classical MD trajectory, significantly more data than most thermodynamic integration calculations. Our computed value for ΔΔGAB gets more precise and accurate (relative to the known experimental result) as more data are included in the calculation, although our computed value (−4.9±0.2 kcal mol−1) remains significantly different to that measured by experiment (−1.0±0.2 kcal mol−1). In forthcoming work we show that for simpler mutations than the one shown here, we are able to achieve agreement with measured affinities to within experimental error. The computed error underestimates the uncertainties introduced by the method for extrapolating values of the free energy derivatives at λ={0,1} and overestimates the convergence of the simulations, although we have shown the latter is reduced when we include more data. Interestingly, the calculation underestimates the error in the difference in binding free energy if simulations that are too short are used. This serves as a reminder of the importance of ensuring that simulations are of sufficient duration.

Approximations made to the Hamiltonian used to describe the system are also a source of error. For example, constraints (e.g. SHAKE, as used here) are often introduced into thermodynamic integration calculations to speed up the simulations, and/or only a small number of atoms within a fixed radius of the binding site are properly modelled. Both approaches introduce error into any calculation. Our Grid-based method will allow some of these sources of error to be examined more rigorously than previously possible.

Understanding the binding of small molecules and peptides to SH2 domains has important clinical applications and we have developed a Grid-based method to accelerate the calculations of differences in binding free energy. The difference in binding free energies between two more similar peptides will be treated in future work. We hope that this Grid-based approach, in combination with experiment, will lead to a better understanding of the binding of peptides to the Src SH2 domain.


Discussions with Dr Shunzhou Wan are gratefully acknowledged and we thank Prof. Gabriel Waksman and his group, especially Sebastien Geroult, for providing insight into the SH2 problem. We thank our RealityGrid collaborators, in particular the SVE group at Manchester University, including Dr Stephen Pickles, Dr Andrew Porter, Dr Robin Pinning and Rob Haines, for developing the Steering Library and associated steering clients. We are grateful to ESPRC for funding much of this research through RealityGrid grant GR/R67699 and for providing a Ph.D. studentship for PWF. Our work was partially supported by the National Science Foundation under NRAC grant MCA04N014, and utilised computer resources at the Pittsburgh Supercomputer Center, the National Computational Science Alliance and the TeraGrid. Finally, this work used both the UKLight and UK National Grid Service.



View Abstract