5.4. Umbrella Sampling Simulations and Free Energy Surfaces

In order to characterize the free energy surface for the exit of ligands 32, 2 and 41 from the FAK binding site, umbrella sampling simulations [40] of the three ligands in complex with FAK were carried out. These simulations were based on exit pathways from initial SMD simulations carried out using low forces in the expectation that these pathways would be more likely to be representative of the actual exit pathway. In the case of ligand 32, the umbrella sampling simulation was based on a separate SMD simulation with a force of 250 pN that had been carried out prior to the group of SMD simulations described above. In the case of ligands 2 and 41, the umbrella sampling simulations were based on the longest SMD exit simulations with a force of 350 pN taken from the main group of SMD simulations. In each case, windows for umbrella sampling were identified by first RMSD aligning the SMD trajectory with the reference structure along the peptide backbone, then choosing frames in which the center of mass of the ligand was 1 Å from the center of mass in any previous window. Each such center of mass was then used as the center of a harmonic umbrella potential for one of the windows. 30 windows were identified in this way for ligand 32 and 22 for ligand 2. In the case of ligand 41, an initial set of 40 windows was identified in this way, but the umbrella sampling simulation did not adequately sample near one of the free energy barriers, so an additional 6 windows had to be added centered on frames from the SMD simulation.

For each ligand, the chosen frames from the SMD simulations were each solvated in a rhombic dodecahedral box of TIP3 water [65] and sodium chloride to an ionic strength of 150 mM, and simulated in explicit solvent with the CHARMM 36 force field [66,67] and the particle mesh Ewald method [68]. This simulation consisted of heating to 300 K and equilibration under protocols similar to those used for the SMD simulations, followed by 20 ns of simulation for each window. Langevin dynamics also used to maintain constant temperature, with a damping constant of 10 ps−1. All simulations were also conducted at a constant pressure of 1 atm, with a barostat oscillation time of 100 fs and a decay time of 50 fs. The production simulations used a harmonic umbrella potential centered on the position of the center of mass in the initial frame taken from the SMD simulation and having the following form:

where UFF(x) is the potential given by the force field, rCM is the center of mass of the ligand after alignment of the protein backbone with the initial frame taken from the SMD simulation, and rCM,0 was the corresponding position taken from that frame. (These positions are shown in Figure 8.) The force constant k of the harmonic restraint was 1.0 kcal/mol Å2. The use of backbone alignment ensured that overall rotation and translation of the protein had no effect on the umbrella potential.

Positions of umbrella potential centers obtained from SMD simulation of ligand 32 exit, superposed on a close up view of the FAK binding site.

These simulations were carried out with NAMD [58,59] and the umbrella potential was imposed using the colvars module within NAMD. The first 2 ns of each trajectory was discarded. From the remaining 18 ns of each trajectory, a three-dimensional histogram of the position of the center of mass of the ligand was constructed with a bin size of 0.25 Å in each dimension. These histograms were then combined into a three-dimensional free energy surface for each ligand using the weighted histogram analysis method [69,70]. This method simultaneously corrected for use of the biasing potential and constructed the minimum variance estimate of the free energy surface by combining the histograms from each simulation using the following equations:

where G(r) and p0(r are the free energy surface and probability density in terms of the center of mass r of the ligand, hj(r) is the histogram obtained from simulation j, Uj(r)=k2rCM−rCM,j is the biasing potential applied to simulation j, and nj is the number of data points for simulation j. ZiZ0 is the ratio of the partition function for simulation i to that of the unbiased simulation. These equations were solved self-consistently by iterating until the values of lnZiZ0 for all simulations i had converged to less than 0.001.

Using VMD, each resulting three-dimensional free energy surface was then visualized as contour surfaces together with the protein. This made it possible to determine the location of free energy basins and saddle points corresponding to transition states in relation to the protein. To determine the locations, the isosurface of constant free energy was visualized for increasing free energy levels until two energy basins could connect with each other, as shown in Figure 1a. Each free energy surface was also used to produce a one-dimensional potential of mean force in terms of the distance of the ligand center of mass from its position in the reference structure. This was done by integrating p0(r), as calculated from Equation (6) over all bins corresponding to a given center of mass distance and converting to a free energy. Due to statistical error in the free energy surface, it is not possible to estimate the transition state free energy with a precision much higher than 0.2–0.3 kcal/mol.