Binding affinity estimation from restrained umbrella sampling simulations

Theoretical foundation

Binding affinity is often quantified using the equilibrium dissociation constant (K_d), defined as:

where [P], [L] and [P:L] are the concentrations of protein, ligand and the protein–ligand complex, respectively. Computationally, the absolute binding free energy (ΔG°), which is the standard molar free energy of binding, is more convenient to calculate. The dissociation constant and the absolute binding free energy are related via

where R is the gas constant, T is the temperature and 1 M is 1 molar concentration. Various strategies have been used to estimate ΔG°, some of which were briefly discussed above. The methodology proposed here has a notable resemblance to the stratification strategy of refs. ^3,4. However, the two methods have major differences as will be discussed later in this section.

Absolute binding free energy or ΔG° is the free-energy change associated with moving the ligand from the bulk to the binding pocket (Supplementary Table 1). Within the formalism presented in this work, ΔG° is determined from the grid PMF G(x), where x is the position of the ligand mass center from the center of the binding pocket (Supplementary Table 1), G(x) is the PMF associated with the ligand position x. In practice, we need to bin the three-dimensional space and define the PMF at every bin or grid point as:

where p(x) is the probability of finding the ligand at bin x.

We define ΔG(x) = G(x) − G(0), where x = 0 (that is, the center of the binding pocket) is defined as the grid point associated with the lowest grid PMF. One can show:

in which the binding ‘pocket’ refers to all x ∈ V where the ligand is considered bound and ‘bulk’ refers to all x ∈ V where the ligand is not interacting with the protein. V here is a subset of space with a single protein in standard concentration (that is, 1 M). As ΔG(x) is the same everywhere in the bulk, we can simplify relation (4) as follows:

where V_B is the bulk volume per protein associated with the standard concentration, x_B is any grid point in the bulk and V_P is the binding pocket volume defined as:

Defining ΔG_V as the contribution of the difference between the volume of the binding pocket and the bulk to the binding free energy:

Combining equations (5) and (7), we have:

We can find the bulk volume (V_B) associated with the standard concentration for a single protein approximately as:

where N_A is Avogadro’s constant and L is the unit of volume (litres). We can now rewrite ΔG_V as:

in which ΔG_B is the bulk volume contribution and ΔG_P is the binding pocket contribution:

Determining both ΔG(x_B) and ΔG_P requires finding the grid PMF ΔG(x). ΔG(x_B) is the PMF difference between the binding pocket center and the bulk and ΔG_P also requires an estimate for ΔG(x) within the binding pocket. We therefore do not need to find ΔG(x) for all x if we have a good estimate for ΔG(x) within the binding pocket and in the bulk. Ideally, ΔG(x) for these points can be determined by pulling the ligand out of the binding pocket towards the bulk and using an enhanced sampling technique such as US to sample the space of a collective variable such as d, that is, the distance between the mass centers of the ligand and protein. ΔG(x) can be estimated for all sampled grid points x using this distance-based US simulation. Note that the collective variable used for biasing would be d, while the collective variable used for the PMF calculations would be the three-dimensional position vector of the mass center of ligand with respect to protein’s binding pocket center. One may estimate the grid PMF from the distance-based US simulations using a non-parametric reweighting algorithm as discussed in this section. ΔG(x) can also be used to estimate ΔG_P as defined in relation (11). There is often no need to strictly define the binding pocket as only low ΔG(x) values have non-negligible contribution to V_P and thus even if we include all sampled grid points, only those close to the binding pocket center have non-negligible contributions.

A practical issue with determining ΔG(x_B) is the convergence. The key obstacles for the sampling that slow down the convergence are the orientation of the ligand, and the conformational changes of the ligand and protein. Using an approach similar in spirit to the previously proposed stratification strategy^3,4,24, we can circumvent extensive sampling of these degrees of freedom. Let us first focus on the orientation of the ligand (Ω), defined using the orientation quaternion formalism. We can restrain Ω during the distance-based US simulations using a biasing potential ((frac12kvarOmega ^2) where a k is harmonic force constant) and later correct the free-energy difference based on the PMF associated with the Ω, which is different in the bulk (F(x_B, Ω)) and in the binding pocket (F(0, Ω)). More generally, for any grid point x, we may determine ΔG(x) based on the PMF associated with the Ω at x (F(x, Ω)) and 0 (F(0, Ω)):

Note that F(x, Ω) is the PMF associated with x and Ω, defined such that:

where c is an arbitrary constant. We therefore have:

We now define G_Ω(x) as the grid PMF of the restrained system (by Ω):

We also define U_Ω(x) as the average biasing potential at grid point x:

Now we have from relations (14), (15) and (16):

where, the free energy of grid point x from the center 0 (ΔG(x)) is calculated based on its equivalent free energy (ΔG_Ω(x)) in a system biased by a harmonic restraint on Ω and a correction term ΔU_Ω(x). For x = x_B:

To determine the above ensemble averages, we need to determine the PMF along Ω for the bound and unbound ligand and calculate the ensemble averages analytically using relation (16). ΔG_Ω(x_B) can be determined from PMF calculations, where the distance between the protein and ligand is varied and the orientation of the ligand is restrained (distance-based BEUS with restrained orientation). We note that:

where we assume ΔU_Ω(x) is negligible for x within the binding pocket. In other words, (langle {mathrme^ – fracfrac12kvarOmega ^2RT} rangle _boldsymbolx approx langle {mathrme^ – fracfrac12kvarOmega ^2RT} rangle _boldsymbol0) for x close to 0.

In brief, if we choose to restrain the orientation, our absolute binding-free-energy estimate includes the following terms (using relations (8) and (17)):

F(x_B, Ω) can be calculated numerically from orientation angle distribution of a free ligand: (F(boldsymbolx_mathrmB,varOmega ) = – RTln pleft( varOmega right),) where p(Ω) is determined from the distribution of Euler angles ((pleft( phi ,theta ,psi right) = frac18uppi ^2sin theta), where 0 ≤ ϕ, ψ ≤ 2π and 0 ≤ θ ≤ π) given that:

(langle {mathrme^ – fracfrac12kvarOmega ^2RT} rangle _{mathrmbulk}) can then be calculated using relation (16) with numerically estimated F(x_B, Ω) and the k value used in the simulations. (Fleft( {boldsymbolx_mathrmB,varOmega } right) = – RTln p(varOmega )) was numerically estimated by discretizing each of the 3 Euler angles with a bin width of 1° and a total of 360 × 360 × 180 bins to estimate p(Ω) from p(ϕ, θ, ψ). F(0, Ω) can be determined approximately using orientation-based US simulations of bound ligand. F(0, Ω) can then be used to estimate (langle {mathrme^{ – fracfrac12kvarOmega ^2RT}} rangle _mathrmpocket) using relation (16).

The above strategy can be extended to other degrees of freedom for which unbiased sampling may hinder the convergence. Most notably, the internal conformational changes of the ligand and that of the protein may also play a crucial role in slowing down the convergence. In the following, we show how one can restrain not only the orientation of the ligand but also the RMSD of the ligand (denoted here by r) in distance-based US simulations (along d) to speed up convergence. In this case, the grid PMF difference ΔG(x) is calculated based on ΔG_Ω,r(x), the grid PMF of a system whose Ω and r are both restrained:

Using a similar strategy as in relation (14), we have:

which results in:

Here we have defined G_Ω,r(x) as:

where k′ is the harmonic force constant associated with the r based on biasing potential ((frac12k^prime r^2)). We also define U_r(x) similar to U_Ω(x) in relation (15) except for using r instead of Ω. (U_varOmega^r(boldsymbolx)) is also defined similar to U_Ω(x) except for the additional restraint on r:

Finally, we have:

In brief, if we choose to restrain both the orientation and RMSD, our absolute binding-free-energy estimate includes the following terms:

Here we are using an approximation similar to that in relation (19):

Using relations (20) and (28), we can generalize the stratification strategy to include three restraints on arbitrary collective variables α, β and γ:

where:

Isothermal titration calorimetry of hFGF1 with heparin hexasaccharide

ITC data were obtained using MicroCal iTC 200 (Malvern) with microcal origin software. The change in heat during the biomolecular interaction was measured by titrating heparin (loaded in the syringe) into the hFGF1 solution in the calorimetric cell. Both the protein and the heparin samples were prepared in the buffer containing 10 mM phosphate buffer with 100 mM NaCl at pH 7.2 and were degassed before loading. The protein-to-heparin ratio was maintained at 1:10 with the protein concentration being 100 μM and the heparin concentration being 1 mM. A total of 30 injections were conducted with a constant temperature of 25 °C and stirring speed of 300 rpm. One set of sites binding model was used for the ITC binding curve⁶⁸. The standard binding free energy ΔG° was determined from dissociation constant via relation (2) at T = 25 °C. The experiment was repeated three times with the same sample and the results obtained were very similar to each other. The mean and standard deviation were reported for both K_d and ΔG°.

All-atom MD simulations

For our bound state, we utilized the X-ray crystal structure of the dimeric hFGF1 combination with heparin hexasaccharide (PDB 2AXM; resolution, 3.0 Å)⁶⁹, and for our apo state, we used the X-ray crystal structure of monomeric hFGF1 (PDB 1RG8; resolution, 1.1 Å)⁷⁰. The NAMD 2.13 (ref. 71) was used to run MD simulations. Using a conjugate gradient, we energy-minimized the system for 10,000 steps. We next relaxed the systems using stepwise restrained MD simulations (for 1 ns) using CHARMM-GUI⁷². All production runs were done in an NPT (constant N, number of atoms; P, pressure; T, temperature) ensemble after the first NVT (constant N, number of atoms; V, volume; T, temperature) relaxation. Simulations were done at 300 K with a 2 fs time step and a 0.5 ps⁻¹ damping coefficient using a Langevin integrator. Nosé–Hoover–Langevin pistons were used to maintain 1 atm pressure⁷². Long-range electrostatic interactions were estimated using the particle mesh Ewald approach. The initial runs were done for 15 ns, followed by the production run on the Anton 2 supercomputer (Pittsburgh Supercomputing Center) for 4.8 μs with a 2.5 fs time step.

MD simulations of free heparin hexasaccharide

Heparin hexasaccharide⁶⁹ was simulated in a rectangular water box without the protein. The system was set up as described previously in the ‘All-atom MD simulations’ section. The final conformation after relaxation was then used as the starting conformation for 10 production runs for 40 ns each. The total simulation time was around 400 ns.

SMD simulations

The final conformations of the hFGF1–heparin⁷³, apo hFGF1 (ref. 73) and free heparin hexasaccharide equilibrium simulations were used to generate starting conformations for the non-equilibrium pulling simulations. Four collective variables⁷⁴ were used for the SMD simulations⁷⁵: (1) distance between the heavy-atom center of mass of heparin and that of the protein (d); (2) the orientation angle of heparin with respect to the protein (Ω) defined using the orientation quanternion formalism; (3) RMSD of the protein (r_P); (4) RMSD of heparin (r_L). Six independent sets of simulations were performed. The distance-based SMD simulation was run for 9.5 ns, while the orientation-based SMD simulation was run for 8 ns. The distance-based SMD simulation was used to pull the heparin away from the protein by approximately 30 Å (10 Å → 40 Å) with a force constant of 100 kcal (mol Ų)⁻¹. The orientation angle was also restrained in this simulation with a force constant of 0.5 kcal (mol degree²)⁻¹ to stay close to its initial orientation in the bound state. The orientation-based SMD simulation was used to rotate the bound heparin locally with respect to the protein (0° → 73°) with a force constant of 100 kcal (mol degree²)⁻¹. Four RMSD-based SMD simulations were run for 10 ns each using a force constant of 50 kcal (mol Ų)⁻¹: (1) to change the RMSD of the bound protein (0.5 Å→2 Å) (the RMSD of heparin was restrained in this simulation with a force constant of 1 kcal (mol Ų)⁻¹); (2) to change the RMSD of the bound heparin (1.5 Å → 4 Å); (3) to change the RMSD of the unbound protein (0.8 Å → 3.2 Å); (4) to change the RMSD of the free heparin (1.5 Å → 5.5 Å).

BEUS simulations

BEUS^53,76,77, which is a variation of the US simulation method, was performed to estimate grid PMF (Supplementary Fig. 1). Four independent sets of distance (d)-based BEUS simulations were performed, with no restraints, restraint on orientation angle of heparin with respect to the protein (Ω), restraint on RMSD of the ligand (r_L) and RMSD of the protein (r_P), and restraints on Ω, r_L and r_P. Two sets of BEUS simulations were also performed using the Ω collective variable, one with and one without a restraint on r_L and r_P. In addition, two sets of BEUS simulations were performed using the r_P collective variable (bound protein with restraint on r_L; unbound protein) and two sets were performed using the r_L collective variable (bound ligand; free ligand). Selected SMD conformations were assigned to individual BEUS windows with equal spacing in each one of these BEUS simulations. The distance-based BEUS simulations ran for 10 ns with 31 replicas/windows and the orientation-based simulations ran for 10 ns with 30 replicas/windows. The RMSD-based BEUS simulations ran for 10 ns with 12 replicas/windows. The force constant used for ligand–protein distance (d) in distance-based BEUS was 2 kcal (mol Ų)⁻¹ while the orientation was restrained as in SMD simulations using a force constant of 0.5 kcal (mol degree²)⁻¹. For orientation-based BEUS simulations, the force constant for the ligand orientation angle (as in SMD simulations) was set to 0.5 kcal (mol degree²)⁻¹. The force constant used for r_L and r_P in all cases was 1 kcal (mol Ų)⁻¹. See Supplementary Fig. 1 for a schematic representation of these simulations.

Free-energy calculations using non-parametric reweighting

Once the BEUS simulations described above were converged, a non-parametric reweighting method^76,78, which is somewhat similar to the multi-state Bennett acceptance ratio method⁷⁹, was used to construct the PMF. In this method⁷⁶, each sampled configuration will be assigned a weight, which can be used to construct the PMF in terms of a desired collective variable. Suppose that a system is biased (for instance, within a BEUS scheme) using N different biasing potentials U_i(r), where i = 1, …, N, and r represents all atomic coordinates. Typically, U_i(r) is a harmonic potential defined in terms of a collective variable with varying centers for different i. Assuming an equal number of sampled configurations from each of the N generated trajectories, we can combine them in a single set of samples r_k (irrespective of which bias was used to generate each sample r_k) and determine the weight of each sample (w_k) as:

where, c is the normalization constant such that (mathop sumnolimits_k w_k = 1) and both w_k and perturbed free energies F_i are determined iteratively using the above equation and the following:

Converged w_k values can be used to construct any ensemble averages including any PMF (for example, G(ζ) PMF of the atomic system in the collective variable space(ζ))) in terms of not only the collective variable used for biasing but also any other collective variables that are sufficiently sampled. One may use a weighted histogram method to construct the PMF as follows:

To estimate the uncertainty of any of PMF calculations described above, one may use bootstrapping. Here, we have used a block Bayesian bootstrapping technique⁷⁷, where 100 alternative datasets are resampled from the existing dataset and the same non-parametric reweighting algorithm and the same PMF calculation is repeated for each set to generate 100 alternative PMFs. The standard deviation of the PMF at any point along the reaction coordinate provides an estimate for the error.

Alchemical FEP simulations

We used the BFEE2 (ref. 27) package to estimate the absolute free energy of binding in silico for an alchemical or geometrical route with multiple subprocesses and geometric constraints. Alchemical FEP simulations were performed to calculate the absolute binding free energy for the interaction of hFGF1 with heparin hexasaccharide. We used a double annihilation protocol⁸⁰, wherein the heparin hexasaccharide is annihilated in both the free and bound states. The final conformations of the hFGF1–heparin complex⁷³ and free heparin hexasaccharide equilibrium simulations (discussed previously in the ‘All-atom MD simulations’ section) were used to generate starting conformations for the bound hFGF1–heparin and free heparin FEP simulations respectively. For the alchemical route, four separate simulations are performed: (1) coupling the restraints of seven collective variables in the bound state; (2) decoupling the ligand alchemically in the bound state; (3) coupling the ligand alchemically in the unbound state; (4) decoupling the conformational restraints in the unbound state. The FEP simulations 1 and 3 were performed bidirectionally using 200 λ-windows (λ is the coupling parameter associated with the FEP that could vary between 0 and 1). Each λ-window included a 0.5 ns of equilibration and 5.0 ns of averaging for both the unbound and bound states, for a total of 2.3 μs (Supplementary Table 3). The decoupling FEP simulations 2 and 4 were also performed bidirectionally, each one for 51 ns. All FEP simulations were performed using the NAMD 2.13 (ref. 71) simulation package with the CHARMM36m all-atom additive force field, using the protocol discussed previously for the equilibrium simulations. We used the state-of-the-art BFEE2 (ref. 27) method to make input files and analyze the FEP simulations.

Binding-free-energy calculations using geometrical route

The extended ABF technique with an umbrella integration estimator was used to calculate the free-energy change along the coarse variables required to characterize reversible heparin–hFGF1 binding^3,24,27. We used the software BFEE2 (ref. 27) to generate the input files for these simulations. In the geometrical route, these collective variables are often subjected to restrictions, and the amount of reversible work required to impose each constraint is determined by a sequence of very accurate PMF simulations. The collective variables used here are the RMSDs of the two proteins’ backbone distances from the reference, native conformation, the three Euler angles (Θ, Φ and Ψ) that describe their relative orientation and the polar (θ) and azimuth angles (φ) that describe their relative position^27,81. The geometrical path consists of a sequence of separate PMF computations performed sequentially with the gradual inclusion of restrictions (RMSD, Θ, Φ, Ψ, θ and φ), as shown in Supplementary Table 2. Each geometric collective variable (RMSD, Θ, Φ, Ψ, θ, φ and r = (1/β) ln(S*I*C°); β = (𝑘_BT)⁻¹, with k_B the Boltzmann constant and T the temperature; C° denotes the standard concentration of 1 M. I*, which stands for the separation term, and S*, which stands for the surface term, indicate the percentage of a sphere with radius r*, centered at the binding site of the reference protein, that is, accessible to its partner) simulation was run with 10 replicas per restriction, and each replica simulation included 20 ns (RMSD, Θ, Φ, Ψ, θ and φ) of simulation time (r collective variables simulations were run for 40 ns for each replica), for a total of 1.6 μs. The BFEE2 (ref. 27) Gui was used to analyze the final ABF simulation data.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.