Using the Fast Fourier Transform in Binding Free Energy Calculations
Trung Hai Nguyen Huan-Xiang Zhou David D. L. Minh
q-bio.BM cond-mat.stat-mech




Phase-matching-free parametric oscillators based on two dimensional semiconductors
    C. Conti^4,5
    December 30, 2023
==================================================================================




According to implicit ligand theory, the standard binding free energy is an exponential average of the binding potential of mean force (BPMF), an exponential average of the interaction energy between the ligand apo ensemble and a rigid receptor.
Here, we use the Fast Fourier Transform (FFT) to efficiently estimate BPMFs by calculating interaction energies as rigid ligand configurations from the apo ensemble are discretely translated across rigid receptor conformations.
Results for standard binding free energies between T4 lysozyme and 141 small organic molecules are in good agreement with previous alchemical calculations based on (1) a flexible complex (R ≈ 0.9 for 24 systems) and (2) flexible ligand with multiple rigid receptor configurations (R ≈ 0.8 for 141 systems).
While the FFT is routinely used for molecular docking, to our knowledge this is the first time that the algorithm has been used for rigorous binding free energy calculations.



Noncovalent Binding Free Energy, Implicit Ligand Theory, Fast Fourier Transform, Protein-Ligand, T4 Lysozyme











  iblabel[1]#1.
  

apsrev










§ INTRODUCTION
 





To accelerate absolute binding free energy calculations between proteins and small molecules, we developed a method based on the fast Fourier transform and then tested it, finding reasonable agreement with more expensive calculations.

The binding of small molecules to biological macromolecules plays a critical role in cellular processes including enzymatic reactions, signal transduction, and gene regulation.
Moreover, noncovalent associations with specific targets are essential to the mechanism of most drugs.
Molecular simulations to estimate standard binding free energies between small molecules and biological macromolecules, therefore, have become increasingly important for understanding biological mechanisms <cit.> and for computer-aided drug design <cit.>.

The most rigorous methods for estimating standard binding free energies are alchemical pathway methods <cit.>, 
in which conformations of both the ligand and receptor are sampled from a series of possibly nonphysical thermodynamic states along a pathway connecting two end states of interest.
These methods can be categorized into two flavors: 
relative and absolute binding free energy methods.
The former calculations attempt to calculate the relative binding free energy of one ligand compared to another. 
In these calculations, the pathway involves transforming one ligand into another.
The latter attempt to calculate binding free energies without comparison to a reference ligand.
In absolute binding free energy calculations, the pathway may involve turning off interactions between the ligand and receptor or physically separating them <cit.>.
Simulations may be done in explicit or implicit <cit.> solvent, with the latter increasing speed at the expense of accuracy <cit.>.
In an increasing number of publications, alchemical pathway methods have successfully calculated accurate protein-ligand binding free energies <cit.>.

If the binding free energy of many ligands is desired, however, 
sampling from a series of thermodynamic states for each receptor-ligand pair can be computationally expensive.
The ability to perform large-scale free energy calculations efficiently is desirable for the purposes of computer-aided drug discovery.
Absolute binding free energy calculations for a large library of diverse compounds may be used in lead discovery.
Relative binding free energy calculations between a congeneric series may be used in lead optimization.

A number of strategies are emerging to facilitate large-scale free energy calculations.
In λ-dynamics <cit.>, the alchemical parameter specifying the identity of a substituent is treated as a dynamical variable.
Its extension to multiple substituents at multiple locations on a scaffold, multi-site λ-dynamics <cit.>, enables many relative binding free energy estimates based on a single simulation. Another strategy involves pre-computing ensembles of the protein complexed with a reference molecule <cit.>, with a reference ligand <cit.>, or of the protein by itself, the apo ensemble <cit.>.
This approach is advantageous because it allows for receptor conformations to be exhaustively sampled once and used for a large set of ligands.
For a series of molecules that are similar to the reference molecule, a one-step perturbation is often sufficient to estimate their relative binding free energies.

Implicit ligand theory (ILT) <cit.> is a formalism that enables large-scale absolute binding free energy calculations based a pre-computed apo ensemble of the receptor.
The theoretical basis for ILT is similar to the basis for estimating standard binding free energies based on implicit solvent <cit.>.
While the implicit solvent formalism involves an integral over all solvent degrees of freedom to obtain a potential of mean force for solvation, 
ILT invokes an integral over all ligand degrees of freedom to obtain a binding potential of mean force (BPMF) - the binding free energy between a flexible ligand and rigid receptor.
The solvent potential of mean force is typically estimated by a continuum dielectric electrostatics model. 
In contrast, the BPMF has hitherto been estimated by an alchemical pathway method <cit.>.

Recently, <cit.> demonstrated the feasibility of using ILT <cit.> (ILT) to estimate absolute binding free energies based on multiple rigid conformations of a protein.
Receptor conformations were first drawn from alchemical binding free energy calculations between T4 lysozyme and 6 different ligands.
The BPMF was then calculated between 141 ligands and each of the receptor conformations.
In accordance with ILT, the standard binding free energy was calculated as an exponential average of the BPMFs.
ILT-based calculations closely reproduced (with Pearson's R = 0.9 and RMSE = 1.59 kcal/mol) results obtained by YANK <cit.>, a program based on alchemical pathway calculations with a flexible receptor.

Here, we perform similar calculations using an alternate way to estimate BPMFs. 
<cit.> computed BPMFs using the program Alchemical Grid Dock <cit.> (AlGDock), which implements an alchemical pathway method in which the receptor is held rigid.
While these calculations were much faster than simulations with a flexible receptor, they still required sampling from multiple thermodynamic states. 
To understand the alternate approach, it is helpful to consider the definition of the BPMF between a flexible ligand and a rigid receptor <cit.>,

    B(r_R) = - β ^ -1ln< I(ζ_L) exp[ - βΨ(r_RL) ]  > _ζ_L, r_L,

where β = (k_B T)^-1.
Internal coordinates of the complex, r_RL, include internal coordinates of the receptor and ligand, r_R and r_L, respectively, and external degrees of freedom of the ligand, ζ_L.
I(ζ_L) is an indicator function between 0 and 1 that specifies whether the receptor and ligand are bound.
Ψ(r_RL) is the effective interaction energy between the receptor and ligand, the difference in energy between the solvated complex and solvated receptor and solvated ligand.
The bracket < ·>_ζ_L, r_L denotes the ensemble average over the external and conformational degrees of freedom of the ligand in the apo ensemble.
In this paper, the apo ensemble of a ligand (or receptor) refers to the probability density of the molecule in the unbound state, alone in implicit solvent.
Equation <ref>, therefore, specifies that BPMFs can be calculated by sampling ligand internal coordinates r_L from the apo ensemble, sampling external degrees of freedom from the distribution proportional to I(ζ), calculating interaction energies, and estimating B(r_R) with a sample mean,

    B̂(r_R) = -β^-1ln1/N∑_n=1^Nexp[ -βΨ(r_RL,n) ].


The fast Fourier transform (FFT) facilitates the use of Equation <ref> providing an efficient way to compute interaction energies as the ligand is translated across the receptor.
If interaction energies are calculated at M translational positions per dimension, 
the complexity of direct calculations is O(M^6).
<cit.> pioneered an alternate procedure based on the FFT 
in which interaction energies are calculated from the cross-correlation between three-dimensional grids that represent each binding partner.
Their grids were very simple - points were 0 outside the protein, 1 on its surface, and a constant value for the interior - but more recent studies have incorporated elements of a molecular mechanics nonbonded interaction energy, including electrostatics <cit.>, van der Waals interactions <cit.>, or both <cit.>.
After mapping binding partners onto the grids, a discrete Fourier transform is performed, the two grids are convoluted, and an inverse Fourier transform is performed, yielding interaction energies for relative translations specified by grid point positions.
Using the FFT, this algorithm is of order M^3 ln(M^3) or less!
Due to this speedup, docking algorithms based on the FFT are routinely used for protein-protein docking <cit.> and have been applied to docking fragments to proteins <cit.>. Qin and Zhou <cit.> have implemented an FFT-based method for calculating an exponential average of the intermolecular interaction energy (similar to Equation <ref>) and determining the chemical potential of a tracer protein in a crowded macromolecular solution.

In this study, we develop a similar FFT-based method for estimating BPMFs and standard binding free energies between T4 lysozyme L99A and the same set of ligands as in <cit.>.
We will refer to our approach as FFTΔG.
To our knowledge, this is the first time that the FFT has been used for rigorous binding free energy calculations.









§ METHODOLOGY








Our calculations were performed on the same systems as in the prior study <cit.>: T4 lysozyme L99A complexed with 141 ligands with experimentally measured activities.
In a thermal denaturation shift assay, 69 of the ligands were determined to be active and 71 inactive <cit.>.
Another of the ligands, iodobenzene, was determined to be active by isothermal titration calorimetry.
We used the same AMBER <cit.> topology and coordinate files for the ligands and receptor as before <cit.>.
BPMFs were calculated between these ligands and the same T4 lysozyme conformations, extracted from alchemical pathway simulations for six ligands: 1-methylpyrrole, benzene, p-xylene, phenol, N-hexylbenzene, and (±)-camphor.



§ LIGAND SAMPLING


To sample ligand conformations from the apo ensemble, we performed parallel tempering <cit.> in the gas phase with a script based on OpenMM 6.3.1 <cit.>. Langevin dynamics simulations were performed at eight geometrically spaced temperatures between 300 and 600 K for 1 ns at each temperature using a time step of 2 fs.
Exchanges between nearest neighbors were attempted every 1 ps.
Ligand snapshots from the simulation at 300 K were saved every 1 ps, 
resulting in the total of 1000 snapshots per ligand. 
In the apo state, the energy is isotropic with respect to translation and rotation.
To sample from this state, therefore, each snapshot was randomly rotated about its centroid.
Translational sampling was based on the FFT.



§ FFT-BASED BPMF ESTIMATION



BPMFs between each ligand and receptor snapshot were estimated based on the thermodynamic cycle shown in Figure <ref>.
Solvent contributions to the BPMF are (1) the desolvation of the rigid receptor and flexible ligand from state (A) to (B), and (2) the solvation of the complex from state (C) to (D).
The gas phase BPMF is the free energy difference between states (B) and (C), 
estimated by using the FFT to calculate interaction energies for a set of discrete translations of the ligand relative to the receptor.



Interaction energies were calculated using grid-based terms for electrostatic and repulsive and attractive van der Waals interactions.
The general form for these terms is <cit.>,


    Ψ_gen = ∑_i ∈ Lγ_i  ∑_j ∈ Rγ_j V(r_ij),


where the inner sum is over all receptor atoms and the outer one over all the ligand atoms. V(r_ij) is a function of the distance r_ij between atoms i and j.
For the electrostatic energy,

    Ψ_ELE = ∑_i ∈ L q_i  ∑_j ∈ R q_j 1/r_ij,

where q_i and q_j are atomic partial charges.
As partial charges were specified in electronic charge units and r_ij in Å, Ψ_ELE was multiplied by 332.05 such that the energy is in units of kcal/mol.
The van der Waals interaction energy is based on the Lennard-Jones potential using the geometric (opposed to arithmetic) mean to combine the well depth of two atoms as ϵ_ij = (ϵ_iiϵ_jj)^1/2 and radii as σ _ij = (σ_iiσ_jj)^1/2.
It includes a repulsive term,

    Ψ_LJr = ∑ _ i ∈ Lϵ ^ 1/2 _ iiσ ^ 6 _ ii∑ _ j ∈ Rϵ ^ 1/2 _ jjσ ^6 _jj/ r ^ 12 _ ij

and an attractive term,

    Ψ_LJa = - ∑ _ i ∈ L 2 ϵ ^ 1/2 _ iiσ ^ 3 _ ii∑ _ j ∈ Rϵ ^ 1/2 _ jjσ ^3 _jj/ r ^ 6 _ ij.

Combining all three terms, the total interaction energy is Ψ_ELE + Ψ_LJr + Ψ_LJa.

Each of these terms was calculated by discretizing Equation <ref> onto a pair of three-dimensional grids: one for the receptor and one for the ligand.
Values on the receptor grid are based on the inner sum of Equation <ref>,

    G_R (n) = ∑_j ∈ Rγ _ j V ( r_j - n),

where n is the position vector of a grid point.
To discretize the ligand grid G_L (n), 
first, the ligand was placed near the corner (i, j, k) = (0,0,0) such that its lowest x, y and z coordinates are 2 grid spacings away from the planes i=0, j=0 and k=0, respectively.
Then the generalized “charges” γ _ i of the ligand are distributed into its ten nearest grid points as described by <cit.>.
With these discretizations, the interaction energy between the protein, in its original position, and the ligand, translated by the vector m, is approximated by the cross-correlation function,

    Ψ_gen(m) ≈ C(m) ≡∑ _ n G_R (n) G_L (n + m) .

where the approximation is due to discretization. According to the cross-correlation theorem, it is possible to evaluate the cross-correlation function C for every position on the grid by,

    C = FT^-1[ FT(G_R) ·FT(G_L)^* ],

where FT and FT^-1 denotes the forward and inverse discrete Fourier transforms, respectively, and the dot symbol · denotes the element-wise product of matrices.
While direct evaluation of grid-based interaction energies for all points has a complexity of O(M^6), the complexity of FFT-based estimates scales as M^3 ln(M^3) or better, where M is number of translations per dimension.

FFT-based interaction energy calculations were performed for two grid sizes.
In one set of calculations, we used a smaller cubic grid of 16 Å along each dimension, centered around the binding site (as defined in <cit.>).
For a subset of 21 ligands in which experimental binding free energies are available, 
we also performed calculations with a larger cubic grid of 62 Å along each dimension, which encompasses the whole receptor surface.

Interaction energies computed directly and by the FFT-based method were compared for three grid spacings: 0.125 Å, 0.25 Å, and 0.5 Å.
Direct calculations based on equations <ref>, <ref>, and <ref> were performed using a custom python script.
Based on the comparison, grid spacings of 0.25 Å and 0.125 Å were used for the small grid and a spacing of 0.25 Å was used for the larger grid.

While the FFT evaluates interaction energies at all grid points, not all of the points are useful for BPMF estimation. 
First, because the FFT is periodic, an interaction energy can be evaluated for an unphysical situation where a ligand is split between opposite sides of a receptor grid.
To filter out configurations corresponding to these unphysical translations, 
a ligand box was defined as a rectangular region containing all ligand coordinates and two additional grid spacings.
Interaction energies corresponding to translations where the ligand box is split across the system were discarded.

Second, the receptor grid given by Equation <ref> implies that when an atom is very close to a grid point n, the repulsive van der Waals term will cause G_R (n) to be an extremely large positive number. This can lead to a floating point overflow in the FFT calculation. 
To avoid extremely large grid values, we set V(r) in Equation <ref> to zero for r values less than the van der Waals parameter 2^-1/6σ _ ii. (Note that setting V(r) to zero for translations that result in steric clash is arbitrary and does not affects the final result because we did not use the interaction energy values given by FFT calculation for these translations. Instead, we assumed that they are infinite.)
To keep track of translational positions that give rise to a steric clash, we used occupancy grids, defined the same way for both the receptor and ligand as,

    G_occ (n) = 
    
    1    if ∃ i: r_i - n < 2^-1/6σ_ii, 
    
    0    otherwise.

The FFT correlation function of the ligand and receptor occupancy grids gives positive values whenever there is at least one pair of atoms that sterically clash. It gives zero otherwise. 
Interaction energies for translations that have steric clashes were assumed to be infinite. 
Interaction energies (including infinite values) were collected for all sampled translations, random rotations, and conformations of the ligand to estimate the BPMF using Equation <ref>.

Solvent contributions to the BPMF were estimated based on a single-step perturbation,

    B̂_j,k(r_R) = -β^-1ln1/N∑_n=1^Nexp[ -βΔ U_j,k(r_RL,n) ],

where B̂_j,k(r_R) is the free energy difference between the sampled state j and target state k and Δ U_j,k(r_RL,n) is the difference in the potential energy of the complex in state k versus state j.
In the target state, the receptor, ligand, or complex, was solvated in 


  * Onufriev-Bashford-Case (OBC2) generalized Born / surface area implicit solvent <cit.>, implemented in OpenMM 6.3.1 <cit.>; or

  * Poisson-Boltzmann / surface area (PBSA) implicit solvent, implemented in the sander program from AmberTools <cit.>. 

For desolvation of the rigid receptor and flexible ligand, the sampled state was (B) and target state was (A).
Since the FFT approach does not directly sample from either state (C) or (D), we performed sampling importance resampling from the conformations in state (B). 
That is, 100 conformations were drawn from state (B) with weight proportional to e^-βΨ, where Ψ is the total interaction energy. The weight, e^-βΨ, is the ratio of probabilities for observing the conformation in the bound state versus the apo state.
After drawing from state (C) using sampling importance resampling, the free energy difference is estimated using Equation <ref> with state (C) as the sampled state and (D) as the target state.



§ STANDARD BINDING FREE ENERGY ESTIMATION


According to implicit ligand theory <cit.>, the standard binding free energy Δ G^∘ is given by an exponential average of BPMFs over the apo ensemble of receptor conformations,


    Δ G^∘ = - β ^-1ln< exp[ - β B(r_R) ] > _ r_R - βln( Ω C^∘/8 π^2) ,


where Ω = ∫ I(ζ_L) dζ_L = 8π^2 V_site. 
The 8π^2 comes from an integral over the rotational degrees of freedom on the ligand. Because there are no restraints on this rotation, the integral is over the full range of Euler angles.
V_site is the volume of the binding site, defined as the rectangular region where the ligand box does not go outside of the receptor grid (see above for the definition of the “ligand box”). Since each FFT translation yields a slightly different box size, the average volume of all rectangular regions is used to calculate V_site. C^∘ is the standard state concentration.

Standard binding free energies were estimated using a weighted sample mean of the BPMFs of different receptor snapshots,


    ΔĜ^∘ =
    -β^-1ln∑_r_R W(r_R) exp[ -β B(r_R) ] + β^-1ln( Ω C^∘/8 π^2),


where W(r_R) is the normalized weight (∑_r_R W(r_R) = 1) associated with each receptor configuration r_R in the apo ensemble.
<cit.> explained in detail why we used the weighted mean and how to obtain the weights. Briefly, receptor conformations were not sampled from the apo state but from a series of alchemical states in the presence of a ligand using the program YANK <cit.>. Therefore we had to reweigh these conformations to the apo ensemble. To this end we used the multistate Bennett Acceptance Ratio <cit.> to obtain weights for all the receptor snapshots in the apo ensemble. However, the set of receptor conformations selected for binding free energy calculations is just a small subset of all the snapshots from YANK simulations. Therefore, there is no obvious way to assign weights for snapshots in the selected subset given the weights of all snapshots.
As in <cit.>, we combined MBAR <cit.> weights using six weighting schemes:



  * Each snapshot is assigned its own MBAR weight; each YANK simulation has equal weight.

  * Each snapshot is assigned its own MBAR weight; each apo state has equal weight.

  * Each snapshot is assigned the cumulative MBAR weight of the thermodynamic state it represents; each YANK simulation has equal weight.

  * Each snapshot is assigned the cumulative MBAR weight of the thermodynamic state it represents; each apo state has equal weight.

  * Each snapshot is assigned the MBAR weight of its neighbors; each YANK simulation has equal weight.

  * Each snapshot is assigned the MBAR weight of its neighbors; each apo state has equal weight.




§ CORRELATION AND ERROR STATISTICS


To quantify the agreement between different data sets, we used the Pearson's R, the root mean square error (RMSE), and adjusted RMSE (aRMSE).
The RMSE between two series of data points {x_1, x_2, ..., x_N} and {y_1, y_2, ..., y_N} is,

    ϵ = √(1/N∑_n=1^N[ x_n - y_n ]^2 ).

The aRMSE is <cit.>,

    ϵ = √(1/N∑_n=1^N[ x_n - y_n - (x̅ - y̅) ]^2 ),

where the x̅ and y̅ are the sample mean of x and y, respectively. The aRMSE accounts for systematic deviation between the series and is useful for assessing whether relative binding free energies are accurate.



§ RESULTS AND DISCUSSION






§ FFT-BASED INTERACTION ENERGY ACCURACY DEPENDS ON GRID SPACING


Distributing generalized ligand charges onto a grid results in a discretization error that is more pronounced for a larger grid spacing (Figure <ref>).
For a grid spacing of 0.125 Å, the FFT interaction energies are very close to the direct computation values and discretization error is essentially nonexistent.
Increasing grid spacing to 0.25, 0.5, and 0.8 Å leads to a gradual reduction of the Pearson's R to 0.89 and increase of the root mean square error to 0.63 kcal/mol (Figure S1 in the supplementary material).

While the accuracy at larger spacing may still be acceptable, increased spacing also reduces the number of samples. 
Therefore, for BPMF calculations with the smaller grid with 16 Å edges, we used grid spacings of 0.25 Å and 0.125 Å. For the larger grid with 62 Å edges, we reduced computational expense by only using a grid spacing of 0.25 Å.

Our free energy calculations employ a finer grid spacing than Qin and Zhou <cit.> used to estimate the excess chemical potential of proteins in a crowded solution. In these studies, grid spacings between 0.15 Å and 0.75 Å were all sufficient to calculate the fraction of proteins that do not clash with a crowding molecule <cit.>. However, chemical potential estimates were found to suffer from a systematic bias that increases with the grid spacing, and is approximately 0.75 kcal/mol at 0.6 Å<cit.>. To improve accuracy while retaining the speed of a larger spacing, <cit.> corrected for these discretization errors by adjusting atomic radii and partial charges; this strategy may be pursued in future work on FFT-based protein-ligand interaction energies.







§ FFT-BASED INTERACTION ENERGIES FAVOR THE BINDING SITE AND PROTEIN SURFACE


A cross-section of the interaction energy grid between benzene and the whole protein (Figure <ref>) shows the lowest values in the binding site and near the protein surface and many infinite values within the protein.
Low values in the binding site are expected, but low values on the surface may indicate either alternative binding sites or an artifact of the force field.
The existence of alternative sites was suggested by <cit.>, who predicted that several ligands including benzene bind to multiple sites on T4 lysozyme and that each site contributes to the binding affinity.
Interaction energies for surface sites may be comparable to the main binding site because the gas phase interaction energy does not account for solvation.
To elaborate, placing the small hydrophobic molecule benzene in the main binding site would be favorable because it would be almost entirely desolvated, but this effect is not captured in our treatment of the interaction energy.
The infinite interaction energies result from steric clashes between the ligand and receptor.
The existence of many translational vectors with steric clashes is common in FFT-based docking, evident from even the first paper in the area <cit.>.






§ BPMF CALCULATIONS APPARENTLY CONVERGE


Ligand conformational and binding pose sampling appears to be thorough.
Parallel tempering and random rotation leads to ligands being uniformly oriented in three-dimensional space (Figure <ref>a). 
When ligand conformations are translated across the receptor binding pocket, 
a large number of translations result in steric clashes and infinite interaction energies.
However, due to the diversity of sampled ligand conformations and orientations, a substantial number of poses have finite interaction energies (Figure <ref>b).





Ligand sampling is sufficient for BPMF estimates to apparently converge (Figure <ref>).
As the number of ligand conformations increases, the standard deviation of the BPMF decreases and, for the ten arbitrarily chosen systems in the figure, starts to level off at around 1000 configurations.
At this point, the standard deviations of the BPMFs range from 0.04 to 0.3 kcal/mol.




The main caveat of this analysis is that exponential averages such as Equation <ref> are known to suffer from convergence issues due to conformational sampling limitations and finite sample bias <cit.>.
As ligand conformations important to the holo ensemble may differ from those in the apo ensemble, 
many ligand conformations may be required to sample from the energetic minima in the holo ensemble.
Moreover, a large number of samples may need to be sampled from within each minima in order to observe rare events that have the largest contribution to the exponential average.
The established approach to address these issues of configuration space overlap and finite sample bias is to compute free energy differences between similar thermodynamic ensembles in an alchemical pathway.
Comparison to the alchemical pathway method in AlGDock allows us to better evaluate the accuracy of the FFT-based single-step perturbation.

FFT-based and AlGDock BPMF estimates are highly correlated but there is an error of 3.31 ± 0.21 kcal/mol  (Figure <ref>). As seen in Figure <ref>b, much of the deviation from AlGDock stems from using a single-step perturbation instead of an alchemical pathway to estimate the gas-phase BPMF. Some error also results from the complex solvation free energy, but the ligand solvation free energy estimates are very consistent.





§ FFTΔG IS CONSISTENT WITH YANK


For 24 ligands binding to T4 lysozyme, using the FFT-based minimum interaction energy as an approximation to the BPMF yields results that are highly correlated with the binding free energy from YANK, but with large error (Figure <ref>).
In the context of ILT, using the minimum interaction energy instead of the exponential average is referred to as the dominant state approximation <cit.>.
When the minimum interaction energy from all receptor structures is used as a free energy estimate, the correlation with free energies is high, with a Pearson's R of 0.90, but the slope is much greater than one and RMSE is large, at 12.70 kcal/mol.
Using an exponential average of BPMFs estimated by the dominant state approximation retains the high correlation with a reduction in RMSE.



The full FFTΔG procedure of using the Equation <ref> yields binding free energy estimates that are highly correlated with YANK and have much smaller error (Figure <ref> and Table <ref>). 
The RMSE for the full set of ligands ranges between 2.07 and 2.31 kcal/mol and Pearson's R between 0.87 and 0.90 (with the exact value depending on the weighting scheme).
The slope of the linear regression line is close to 1 (≈ 1.2).
When considering only the 11 active ligands, the RMSE is only slightly different from the full data set but the Pearson's R much lower and has greater uncertainty.
This large uncertainty in the Pearson's R likely results from the relatively small spread of free energies among active ligands.
On the other hand, due to a large spread in YANK free energies of inactive ligands, the subset of 13 inactive ligands achieves even higher correlation (Pearson's R ≈ 0.96) than the full set of 24 ligands.








A notable outlier is N-hexylbenzene, whose free energy is estimated to be around -2.5 and -8 kcal/mol by FFTΔG and YANK, respectively. N-hexylbenzene is the largest active ligand and possesses a long and flexible hydrocarbon chain. As it is difficult to sample the bound conformation of this chain from a simulation of the ligand by itself, its binding free energy is overestimated.

The performance of different weighting schemes are fairly similar to one another.
Overall, weighting scheme (d) gives the lowest RMSE for the full set of ligands and sets containing only active or inactive ligands (see Tab. <ref>). 
However, the differences among weighting schemes in the present work are not significant due to the sizable error bars of the RMSE and Pearson's R.
To be consistent with <cit.>, 
who found that scheme (c) led to the greatest consistency between AlGDock and YANK results, 
we will hitherto use weighting scheme (c).

In recapitulating YANK, FFTΔG is less accurate than AlGDock.
AlGDock yields RMSEs in the range 1.49 to 1.76 kcal/mol <cit.>. 
This reduction in accuracy may be linked to ligand being sampled only from the apo ensemble in the FFT approach as opposed to a thorough sampling from a series of alchemical states in AlGDock method <cit.>. 
With both methods, the slope is slightly greater than one, likely because we used the same receptor conformations and these conformations did not include those relevant to the weakest-binding ligands. 
This makes the weakest-binding ligands appear even weaker, leading to an increased slope.



§ CONSISTENCY WITH YANK REQUIRES A DIVERSE SET OF RECEPTOR SNAPSHOTS


The consistency between FFTΔG and YANK is sensitive to the choice of receptor snapshots.
If only snapshots from a single alchemical simulation are used, 
there is generally a weaker correlation and larger RMSE (and uncertainty in the RMSE) compared to YANK (Table <ref> and Figure S3).
However, there are exceptions. When only snapshots from alchemical simulations of T4 lysozyme in complex with p-xylene or (±)-camphor are used, FFTΔG gives rather strong correlation and relatively low RMSE compared to YANK. 
These two ligands are the second largest active and the largest inactive ligands, respectively. 
The quality of these calculations can be attributed to the fact that larger ligands tend to open up the binding pocket more widely and hence improve the FFT translational sampling of ligands. 
However, N-hexylbenzene, the largest active ligand with a long carbon chain, gives the worst FFT-based free energy estimates. 
YANK simulations with this ligand may have induced a binding pocket shape that is unfavorable for binding to most of other ligands.
These trends are consistent with our previous results <cit.>, although N-hexylbenzene snapshots did not perform as poorly when using AlGDock.

In contrast to <cit.>, using only snapshots from alchemical simulations of T4 lysozyme complexes with four active ligands (Tab. <ref>) results in worse consistency with respect to YANK than using all snapshots. This suggests that FFTΔG may require a larger set of receptor snapshots to obtain converged binding free energy estimates, likely due to the large number of steric clashes.
 





§ FFTΔG CORRELATES WITH ALGDOCK


For 141 ligands, binding free energy estimates based on FFTΔG and AlGDock are highly correlated (Figure <ref> and Figure S3 in the supplementary material).
The RMSE for 140 ligands (excluding one ligand with high free energies) is 2.24 (0.13) kcal/mol and Pearson's R is 0.82 (0.04). 
The subsets of 69 active and 71 inactive ligands maintain essentially the same level of consistency with AlGDock results with RMSE/Pearson's R at 0.70 (0.08) / 2.23 (0.18) kcal/mol for active ligands and 0.86 (0.04) / 2.24 (0.19) kcal/mol for inactive ligands.
For the full set of 141 ligands, the performance is similar (Figure S3 in the supplementary material).
In the majority of deviations between the two estimates, FFTΔG has a higher value.
The likely cause of these deviations is that the FFT-based procedure does not sample poses with sufficiently low interaction energies.



Using a larger grid spacing of 0.25 Å does not appear to have a significant effect on the accuracy of FFTΔG with respect to YANK, but reduces accuracy with respect to AlGDock (Figures S2 and S4 in the supplementary material). With respect to YANK, the Pearson's R is not much lower and RMSE not much higher. However, accuracy is more sensitive to the weighting scheme. On the other hand, for the larger data set of 141 ligands, the Pearson's R and RMSE with respect to AlGDock is 0.73 and 3.28 kcal/mol, respectively; the RMSE is slightly higher than with 0.125 Å grid spacing. Because interaction energies are fairly accurate at 0.25 Å grid spacing, the likely cause of the deviation is reduced sampling of translational positions.

The overall consistency with respect to YANK and AlGDock suggests that it is feasible to estimate protein-ligand binding free energies using FFTΔG.



§ FFTΔG CONVERGENCE REQUIRES MORE RECEPTOR SAMPLING THAN ALGDOCK


Compared to AlGDock, FFTΔG requires a greater number of receptor snapshots in order for Δ G^∘ estimates to converge. If receptor snapshots are randomly selected, AlGDock requires about 200 snapshots in order for the free energy to converge <cit.>. In contrast, the Pearson's R and RMSE values do not level off until about 400 snapshots processed by FFTΔG (red curves in Figure <ref>). When receptor snapshots are selected in the order of increasing the DOCK 6 scores, the convergence of Pearson's R and RMSE is improved significantly (green curves in Figure <ref>). They all level off after about 100 snapshots. Nevertheless, the convergence of FFT method is still much slower that AlGDock <cit.>, which needs only a few snapshots to recover the best possible correlation and RMSE when snapshots were selected from lowest to highest docking scores. The best convergence can be achieved by selecting receptor snapshots in the order of increasing BPMF as shown by blue curves in Figure <ref>. These curves serve as an ideal unobtainable reference point because a priori, we do not know which snapshots in the set give lowest BPMFs unless we carry out the calculation for all of them. The convergence analysis done here and in <cit.> suggest that it is better to use docking scores to sort out receptor snapshots before performing BPMF calculations. Convergence is likely worse with FFTΔG than with AlGDock because the ligand has no flexibility and must truly fit as a lock-and-key into the selected receptor conformations to obtain accurate free energies.





§ FFTΔG IS CONSISTENT WITH EXPERIMENT


Binding free energy estimates based on FFTΔG are moderately consistent with experimental results (Figure <ref>), with similar but slightly worse performance than AlGDock. 
Isothermal titration calorimetry has been used to measure binding free energies between T4 lysozyme L99A and 21 active ligands <cit.>.
Calculated binding free energies for iodine-containing ligands, as in <cit.>, were found to be highly overestimated.
Therefore, iodobenzene was excluded from the present analysis (but is included in the supplementary material).
The correlation of FFTΔG estimates and measured values is dependent on the solvent model.
Comparing the results from the OBC2 model and experiment, the RMSE is 1.29 ± 0.14 kcal/mol, aRMSE is 1.08 ± 0.08 kcal/mol, and Pearson's R is 0.49 ± 0.17. 
When the PBSA implicit solvent model is used, the RMSE is 3.42 ± 0.42 kcal/mol, aRMSE is 1.99 ± 0.28 kcal/mol, and Pearson's R is 0.52 ± 0.19, which has slightly worse error but comparable correlation. 
For a point of comparison, <cit.> used the PBSA implicit solvent model in AlGDock to attain a RMSE of 2.81 ± 0.32, aRMSE of 1.35 ± 0.27, and Pearson's R of 0.65 ± 0.05.
In both the present results and in <cit.>, the reason for the poor RMSE is that the PBSA model causes a positive shift in the estimated binding free energies.



Using alternate estimators for the binding free energy increases the error but does not significantly reduce the correlation with experiment. 
The error is largest using the dominant state approximation, the minimum interaction energy, as the BPMF.
Error is reduced by taking the exponential average of the minimum interaction energy, and is least with FFTΔG.
These results are consistent with comparisons to YANK, where we found that these approximations also increase the error with respect to YANK calculations.
They are also consistent with <cit.>, who found that using the exponential average is no better than the minimum energy for T4 lysozyme.
However, they found that in other, more complex systems, an exponential average was superior to the minimum interaction energy for reproducing experimental free energies.

We hypothesized that agreement with experiment could be improved by using a larger grid (62 Å on each side) to account for binding to alternative sites. This grid size encompasses the whole surface of the receptor. However, in spite of being more computationally demanding, the FFT calculations with large grid size do not lead to improved agreement with experiment (Figure <ref>). 
When the OBC2 solvent model is used, the Pearson's R, RMSE, and aRMSE with respect to experiment are 0.43 ± 0.18), 2.30 ± 0.19 kcal/mol, and 1.04 ± 0.12 kcal/mol, respectively.
When PBSA solvent model is used, the Pearson's R slightly increases to 0.53 ± 0.13 but the RMSE becomes slightly larger, 3.50 ± 0.24 kcal/mol, and the aRMSE is 1.08 ± 0.16 kcal/mol.
The lack of improvement in agreement with experiment suggests that weak binding of these ligands outside of the L99A cavity may not make a significant contribution to their binding free energy.





§ FFTΔG IS A SLIGHTLY BETTER BINARY CLASSIFIER THAN MOLECULAR DOCKING


Binding free energy methods can also be used to predict whether a molecule is active or inactive against a target.
Here we use receiver operating characteristic (ROC) <cit.> curves, the area under the ROC curve (AUC), and the area under the semi-log ROC curve (AUlC) to assess the ability of FFT free energy estimates to discern active from inactive molecules. 
A ROC curve illustrates the fraction of true positives versus the fraction of false positives as the threshold separating two categories is changed. 
An ideal ROC consists of a vertical line from (0,0) to (0,1), and then a horizontal line from (0,1) to (1,1), meaning that all active molecules are more highly ranked than any inactive molecules. 
The AUC ranges from 0 for completely incorrect to 0.5 for random to 1 for completely correct classification. The intent of the AUlC <cit.> metric is to emphasize top-ranked molecules, which are more likely to be pursued in subsequent experiments and calculations. For a random classifier, the AUlC is 0.14.







These different metrics show that DOCK 6 is essentially random, FFTΔG is slightly better, and AlGDock is the best binary classifier (Figure <ref> and Table <ref>).
With the OBC2 model, ROC curves of AlGDock and FFT methods (based on the minimum interaction energy or BPMF) are about the same and slightly better than DOCK 6. With the PBSA implicit solvent model, AlGDock has better, but FFTΔG has worse ROC curves. The relatively poor binary classification performance of FFTΔG with the PBSA implicit solvent model is consistent with its relatively poor correlation with experimental Δ G^∘ measurements.

As observed in <cit.>, computing a binding free energy opposed to the minimum or mean interaction energy has no benefit to binary classification of active and inactive ligands to T4 lysozyme L99A. As the ligands for T4 lysozyme are small and relatively rigid, the entropic contribution to binding may not change much from ligand to ligand. 



§ FFTΔG IS FASTER THAN ALGDOCK FOR COARSE AND SMALL GRIDS


While FFTΔG is less accurate than AlGDock, its main potential benefit is speed. To compare the computational cost of FFTΔG and AlGDock, we carried out BPMF calculations for a small and a large ligand. Both methods basically consist of two steps: sampling and post-processing with implicit solvent. The post-processing step happens in the same way for both methods and expected to consume more or less the same amount of CPU time. Therefore we only compare timing for the sampling step. For the purposes of this comparison, we sampled the binding pose using AlGDock <cit.> based on Hamiltonian replica exchange as previously described <cit.> and  in FFTΔG using random rotation and the cross-correlation.



Our benchmark calculations (Table <ref>) show that FFTΔG is faster than AlGDock for coarse (0.25 Å spacing) and small grids (16 Å on each edge) and is less sensitive to the ligand size. For a finer grid with 0.125 Å spacing, the speed of FFTΔG is comparable or slower than AlGDock. It is considerably slower for a large grid, but the comparison is unfair because the AlGDock calculation is restricted to smaller binding site opposed to the entire protein surface. One apparent advantage of the FFT approach is that it is less sensitive to the size of the ligand. That is, the larger ligand, N-hexylbenzene, requires approximately twice the amount of computer time with AlGDock, but takes around the same amount of time with FFTΔG. The reason that the calculations take about the same amount of time is that they depend on the number of ligand grid points, not the number of ligand atoms.

A comparison between FFTΔG and YANK is less meaningful because YANK calculations were performed on graphical processing units. As mentioned in <cit.>, each YANK calculation took about one week on a single graphical processing unit.



§ CONCLUSIONS




We have demonstrated that the FFT may be used to estimate standard binding free energies based on implicit ligand theory.
To our knowledge, this is the first time that the FFT has been used to calculate noncovalent binding free energies.
For the binding of T4 lysozyme to 141 small molecules, FFTΔG is less accurate than AlGDock at reproducing YANK, less correlated with experiment, and less capable of classifying active and inactive molecules.
With small grids, however, FFT-based sampling is considerably faster than AlGDock and less sensitive to ligand size.

Given the benefits and limitations of FFTΔG, the approach is mostly likely to be useful for larger chemical libraries than AlGDock.
Virtual screening can be performed by first using FFTΔG on a large library. 
Subsequently, AlGDock may be used for a refined library with the greatest predicted affinity. 
Finally, a flexible-receptor technique such as YANK may be used on yet a smaller library prior to experimental validation.

Reflecting a niche of FFT-based docking, another possible use of FFTΔG may be for computing binding free energies between relatively inflexible fragments and a protein. The FFT is an efficient approach for translating the ligand across an entire protein surface. ILT provides a rigorous formalism for combining interaction energies from multiple sites into a single binding free energy.

Due to its relative insensitivity to the number of ligand atoms, FFTΔG is likely to be useful for larger ligands than AlGDock.
However, FFTΔG may not be efficient if the large ligands are too flexible.
Ultimately, another niche of FFTΔG may be the same as the present niche of FFT-based molecular docking: docking of folded protein domains to each other.
In these cases, the ligands are large but relatively ordered; it may be feasible to recapitulate their flexibility with molecular dynamics simulations.

As a closing comment, we would like to point out a possible way to improve the scoring function of FFT-based docking methods. 
The current paradigm in FFT-based docking, as in other docking strategies, is to search for the configuration with the lowest interaction energy between binding partners. 
Once this position is obtained, all other interaction energies are ignored and the binding strength is scored solely based on the lowest value.
We suggest an alternative approach that requires essentially no additional computational expense:
following the approach of our present study, one may calculate an exponential average of the interaction energy and use the BPMF to score the strength of binding.



§ ACKNOWLEDGMENTS




We thank Paula Bianca Viana Pinheiro for participating in this project as an undergraduate summer intern through the Brazil Scientific Mobility Program.
We thank OpenEye scientific software for providing a free academic license to their software.
Computer resources were provided by the Open Science Grid <cit.>.
This research was supported by the National Institutes of Health (R15GM114781 to DDLM and R35GM118091 to HXZ).

Additional Supporting Information may be found in the online version of this article.
These include figures comparing  
FFT-based vs direct interaction energies with a grid spacing of 0.8 Å, 
binding free energies for 24 ligands estimated using YANK and FFTΔG with a  grid spacing of 0.25 Å,
binding free energies for 141 ligands estimated using AlGDock and FFTΔG with grid spacings of 0.125 and 0.25 Å,
and FFT free energy estimates with experiment for 21 ligands (including iodobenzene) for different grid sizes.








69
natexlab#1#1bibnamefont#1#1bibfnamefont#1#1citenamefont#1#1url<#>1urlprefixURL 


[Dadarlat and Skeel(2011)]Dadarlat2011
authorV. M. Dadarlat and
  authorR. D. Skeel,
  journalBiophys. J. volume100,
  pages469 (year2011).

[Cong et al.(2015)Cong, Campomanes,
  Kless, and Schapitz]Cong2015
authorX. Cong,
  authorP. Campomanes,
  authorA. Kless, and
  authorI. Schapitz,
  journalPLoS ONE volume10,
  pagese0135998 (year2015).

[Lee et al.(2015)Lee, Seok, and
  Im]Lee2015a
authorH. S. Lee,
  authorC. Seok, and
  authorW. Im, journalJ.
  Chem. Theory Comput. volume11, pages1255
  (year2015).

[Calabro et al.(2016)Calabro, Woods,
  Powlesland, Mey, Mulholland, and Michel]Calabro2016
authorG. Calabro,
  authorC. J. Woods,
  authorF. Powlesland,
  authorA. S. J. S. Mey,
  authorA. J. Mulholland,
  and authorJ. Michel,
  journalJ. Phys. Chem. B volume120,
  pages5340 (year2016).

[Li and Nam(2017)]Li2017
authorY. Li and
  authorK. Nam,
  journalChem. Sci. volume8,
  pages3453 (year2017).

[Jorgensen(2004)]Jorgensen2004
authorW. L. Jorgensen,
  journalScience volume303,
  pages1813 (year2004).

[Michel and Essex(2010)]Michel2010
authorJ. Michel and
  authorJ. W. Essex,
  journalJ. Comput.-Aided Mol. Des. volume24,
  pages639 (year2010).

[Chodera et al.(2011)Chodera, Mobley,
  Shirts, Dixon, Branson, and Pande]Chodera2011
authorJ. D. Chodera,
  authorD. L. Mobley,
  authorM. R. Shirts,
  authorR. W. Dixon,
  authorK. Branson, and
  authorV. S. Pande,
  journalCurr. Opin. Struct. Biol. volume21,
  pages150 (year2011).

[Mobley and Klimovich(2012)]Mobley2012
authorD. L. Mobley and
  authorP. V. Klimovich,
  journalJ. Chem. Phys. volume137,
  pages230901 (year2012).

[Parenti and Rastelli(2012)]Parenti2012
authorM. D. Parenti and
  authorG. Rastelli,
  journalBiotechnol. Adv. volume30,
  pages244 (year2012).

[Gilson et al.(1997)Gilson, Given, Bush,
  and McCammon]Gilson1997
authorM. K. Gilson,
  authorJ. A. Given,
  authorB. L. Bush, and
  authorJ. A. McCammon,
  journalBiophys. J. volume72,
  pages1047 (year1997).

[Deng and Roux(2009)]Deng2009
authorY. Deng and
  authorB. Roux, journalJ.
  Phys. Chem. B volume113, pages2234
  (year2009).

[Bekker et al.(2017)Bekker, Kamiya,
  Araki, Fukuda, Okuno, and Nakamura]Bekker2017
authorG.-j. Bekker,
  authorN. Kamiya,
  authorM. Araki,
  authorI. Fukuda,
  authorY. Okuno, and
  authorH. Nakamura,
  journalJ. Chem. Theory Comput. volume13,
  pages2389 (year2017).

[Heinzelmann et al.(2017)Heinzelmann,
  Henriksen, and Gilson]Heinzelmann2017
authorG. Heinzelmann,
  authorN. M. Henriksen,
  and authorM. K.
  Gilson, journalJ. Chem. Theory Comput.
  volume13, pages3260 (year2017).

[Gallicchio et al.(2010)Gallicchio,
  Lapelosa, and Levy]Gallicchio2010
authorE. Gallicchio,
  authorM. Lapelosa, and
  authorR. M. Levy,
  journalJ. Chem. Theory Comput. volume6,
  pages2961 (year2010).

[Wang et al.(2013)Wang, Chodera, Yang,
  and Shirts]Wang2013
authorK. Wang,
  authorJ. D. Chodera,
  authorY. Yang, and
  authorM. R. Shirts,
  journalJ. Comput.-Aided Mol. Des. volume27,
  pages989 (year2013).

[Michel and Essex(2008)]Michel2008
authorJ. Michel and
  authorJ. W. Essex,
  journalJ. Med. Chem. volume51,
  pages6654 (year2008).

[Boyce et al.(2009)Boyce, Mobley,
  Rocklin, Graves, Dill, and Shoichet]Boyce2009
authorS. E. Boyce,
  authorD. L. Mobley,
  authorG. J. Rocklin,
  authorA. P. Graves,
  authorK. A. Dill, and
  authorB. K. Shoichet,
  journalJ. Mol. Biol. volume394,
  pages747 (year2009).

[Ge and Roux(2010)]Ge2010
authorX. Ge and
  authorB. Roux, journalJ.
  Phys. Chem. B volume114, pages9525
  (year2010).

[Wang et al.(2012)Wang, Berne, and
  Friesner]Wang2012
authorL. Wang,
  authorB. J. Berne, and
  authorR. A. Friesner,
  journalProc. Natl. Acad. Sci. USA volume109,
  pages1937 (year2012).

[Zhu et al.(2013)Zhu, Travis, and
  Elcock]Zhu2013
authorS. Zhu,
  authorS. M. Travis,
  and authorA. H.
  Elcock, journalJ. Chem. Theory Comput.
  volume9, pages3151 (year2013).

[Wang et al.(2015)Wang, Wu, Deng, Kim,
  Pierce, Krilov, Lupyan, Robinson, Dahlgren, Greenwood et al.]Wang2015
authorL. Wang,
  authorY. Wu,
  authorY. Deng,
  authorB. Kim,
  authorL. Pierce,
  authorG. Krilov,
  authorD. Lupyan,
  authorS. Robinson,
  authorM. K. Dahlgren,
  authorJ. Greenwood,
  et al., journalJ. Am. Chem. Soc.
  volume137, pages2695 (year2015).

[Aldeghi et al.(2016)Aldeghi, Heifetz,
  Bodkin, Knapp, and Biggin]Aldeghi2016
authorM. Aldeghi,
  authorA. Heifetz,
  authorM. J. Bodkin,
  authorS. Knapp, and
  authorP. C. Biggin,
  journalChem. Sci. volume7,
  pages207 (year2016).

[Aldeghi et al.(2017)Aldeghi, Heifetz,
  Bodkin, Knapp, and Biggin]Aldeghi2017
authorM. Aldeghi,
  authorA. Heifetz,
  authorM. J. Bodkin,
  authorS. Knapp, and
  authorP. C. Biggin,
  journalJ. Am. Chem. Soc. volume139,
  pages946 (year2017).

[Wan et al.(2017)Wan, Bhati, Zasada,
  Wall, Green, Bamborough, and Coveney]Wan2017
authorS. Wan,
  authorA. P. Bhati,
  authorS. J. Zasada,
  authorI. Wall,
  authorD. Green,
  authorP. Bamborough,
  and authorP. V.
  Coveney, journalJ. Chem. Theory Comput.
  volume13, pages784 (year2017).

[Kong and Brooks III(1996)]Kong1996
authorX. Kong and
  authorC. L. Brooks III,
  journalJ. Chem. Phys. volume105,
  pages2414 (year1996).

[Knight and Brooks III(2011)]Knight2012
authorJ. L. Knight and
  authorC. L. Brooks III,
  journalJ. Chem. Theory Comput. volume7,
  pages2728 (year2011).

[Hayes et al.(2017)Hayes, Armacost,
  Vilseck, and Brooks III]Hayes2017
authorR. L. Hayes,
  authorK. A. Armacost,
  authorJ. Z. Vilseck,
  and authorC. L. Brooks
  III, journalJ. Phys. Chem. B volume121,
  pages3626 (year2017).

[Mark et al.(1995)Mark, Xu, Liu, and Van
  Gunsteren]Mark1995
authorA. E. Mark,
  authorY. Xu,
  authorH. Liu, and
  authorW. F. Van Gunsteren,
  journalActa Biochim. Pol. volume42,
  pages525 (year1995).

[Oostenbrink and van Gunsteren(2005)]Oostenbrink2005
authorC. Oostenbrink and
  authorW. F. van Gunsteren,
  journalProc. Natl. Acad. Sci. USA volume102,
  pages6750 (year2005).

[Raman et al.(2012)Raman, Vanommeslaeghe,
  and Mackerell]Raman2012
authorE. P. Raman,
  authorK. Vanommeslaeghe,
  and authorA. D.
  Mackerell, journalJ. Chem. Theory Comput.
  volume8, pages3513 (year2012).

[Minh(2012)]Minh2012
authorD. D. L. Minh,
  journalJ. Chem. Phys. volume137,
  pages104106 (year2012).

[Xie et al.(2017)Xie, Nguyen, and
  Minh]Xie2017
authorB. Xie,
  authorT. H. Nguyen,
  and authorD. D. L.
  Minh, journalJ. Chem. Theory Comput.
  volume13, pages2930 (year2017).

[Minh(2015)]Minh2015
authorD. D. L. Minh,
  journalarXiv pages1507.03703v1
  (year2015).

[Katchalski-Katzir
  et al.(1992)Katchalski-Katzir, Shariv, Eisenstein, Friesem,
  Aflalo, and Vakser]Katchalski-Katzir1992
authorE. Katchalski-Katzir,
  authorI. Shariv,
  authorM. Eisenstein,
  authorA. A. Friesem,
  authorC. Aflalo, and
  authorI. A. Vakser,
  journalProc. Natl. Acad. Sci. USA volume89,
  pages2195 (year1992).

[Gabb et al.(1997)Gabb, Jackson, and
  Sternberg]Gabb1997
authorH. A. Gabb,
  authorR. M. Jackson,
  and authorM. J. E.
  Sternberg, journalJ. Mol. Biol.
  volume272, pages106 (year1997).

[Bliznyuk and Gready(1999)]Bliznyuk1999
authorA. A. Bliznyuk and
  authorJ. E. Gready,
  journalJ. Comput. Chem. volume20,
  pages983 (year1999).

[Qin and Zhou(2013)]Qin2013
authorS. Qin and
  authorH.-X. Zhou,
  journalJ. Chem. Theory Comput. volume9,
  pages4633 (year2013).

[Qin and Zhou(2014)]Qin2014
authorS. Qin and
  authorH.-X. Zhou,
  journalJ. Chem. Theory Comput. volume10,
  pages2824 (year2014).

[Ritchie and Kemp(2000)]Ritchie2000
authorD. W. Ritchie and
  authorG. J. Kemp,
  journalProteins: Struct., Funct., Genet.
  volume39, pages178 (year2000).

[Mandell et al.(2001)Mandell, Roberts,
  Pique, Kotlovyi, Mitchell, Nelson, Tsigelny, and Ten Eyck]Mandell2001
authorJ. G. Mandell,
  authorV. A. Roberts,
  authorM. E. Pique,
  authorV. Kotlovyi,
  authorJ. C. Mitchell,
  authorE. Nelson,
  authorI. Tsigelny, and
  authorL. F. Ten Eyck,
  journalProtein Eng. volume14,
  pages105 (year2001).

[Chen et al.(2003)Chen, Li, and
  Weng]Chen2003
authorR. Chen,
  authorL. Li, and
  authorZ. Weng,
  journalProteins: Struct., Funct., Genet.
  volume52, pages80 (year2003).

[Eisenstein and
  Katchalski-Katzir(2004)]Eisenstein2004
authorM. Eisenstein and
  authorE. Katchalski-Katzir,
  journalC. R. Biol. volume327,
  pages409 (year2004).

[Kozakov et al.(2006)Kozakov, Brenke,
  Comeau, and Vajda]Kozakov2006
authorD. Kozakov,
  authorR. Brenke,
  authorS. R. Comeau,
  and authorS. Vajda,
  journalProteins: Struct., Funct., Bioinf.
  volume406, pages392 (year2006).

[Moal and Bates(2010)]Moal2010
authorI. H. Moal and
  authorP. A. Bates,
  journalInt. J. Mol. Sci. volume11,
  pages3623 (year2010).

[Brenke et al.(2009)Brenke, Kozakov,
  Chuang, Beglov, Hall, Landon, Mattos, and Vajda]Brenke2009
authorR. Brenke,
  authorD. Kozakov,
  authorG.-Y. Chuang,
  authorD. Beglov,
  authorD. Hall,
  authorM. R. Landon,
  authorC. Mattos, and
  authorS. Vajda,
  journalBioinformatics volume25,
  pages621 (year2009).

[Ngan et al.(2012)Ngan, Bohnuud,
  Mottarella, Beglov, Villar, Hall, Kozakov, and Vajda]Ngan2012
authorC. H. Ngan,
  authorT. Bohnuud,
  authorS. E. Mottarella,
  authorD. Beglov,
  authorE. A. Villar,
  authorD. R. Hall,
  authorD. Kozakov, and
  authorS. Vajda,
  journalNucleic Acids Res. volume40,
  pagesW271 (year2012).

[Morton et al.(1995)Morton, Baase, and
  Matthews]Morton1995
authorA. Morton,
  authorW. A. Baase, and
  authorB. W. Matthews,
  journalBiochemistry volume34,
  pages8564 (year1995).

[Morton and Matthews(1995)]Morton1995a
authorA. Morton and
  authorB. W. Matthews,
  journalBiochemistry volume34,
  pages8576 (year1995).

[Su et al.(2001)Su, Lorber, Weston,
  Baase, Matthews, and Shoichet]Su2001
authorA. I. Su,
  authorD. M. Lorber,
  authorG. S. Weston,
  authorW. A. Baase,
  authorB. W. Matthews,
  and authorB. K.
  Shoichet, journalProteins: Struct., Funct., Genet.
  volume42, pages279 (year2001).

[Wei et al.(2002)Wei, Baase, Weaver,
  Matthews, and Shoichet]Wei2002
authorB. Q. Wei,
  authorW. Baase,
  authorL. Weaver,
  authorB. Matthews, and
  authorB. K. Shoichet,
  journalJ. Mol. Biol. volume322,
  pages339 (year2002).

[Graves et al.(2005)Graves, Brenk, and
  Shoichet]Graves2005
authorA. P. Graves,
  authorR. Brenk, and
  authorB. K. Shoichet,
  journalJ. Med. Chem. volume48,
  pages3714 (year2005).

[Mobley et al.(2007)Mobley, Graves,
  Chodera, McReynolds, Shoichet, and Dill]Mobley2007
authorD. L. Mobley,
  authorA. P. Graves,
  authorJ. D. Chodera,
  authorA. C. McReynolds,
  authorB. K. Shoichet,
  and authorK. A. Dill,
  journalJ. Mol. Biol. volume371,
  pages1118 (year2007).

[Graves et al.(2008)Graves, Shivakumar,
  Boyce, Jacobson, Case, and Shoichet]Graves2008
authorA. P. Graves,
  authorD. M. Shivakumar,
  authorS. E. Boyce,
  authorM. P. Jacobson,
  authorD. A. Case, and
  authorB. K. Shoichet,
  journalJ. Mol. Biol. volume377,
  pages914 (year2008).

[Case et al.(2017)Case, Cerutti, T.E.
  Cheatham, Darden, Duke, Giese, Gohlke, Goetz, Greene, Homeyer
  et al.]Case2017
authorD. Case,
  authorD. Cerutti,
  authorI. T.E. Cheatham,
  authorT. Darden,
  authorR. Duke,
  authorT. Giese,
  authorH. Gohlke,
  authorA. Goetz,
  authorD. Greene,
  authorN. Homeyer,
  et al., titleAMBER 2017
  (year2017).

[Swendsen and Wang(1986)]Swendsen1986
authorR. H. Swendsen and
  authorJ.-S. Wang,
  journalPhys. Rev. Lett. volume57,
  pages2607 (year1986).

[Sugita and Okamoto(1999)]Sugita1999
authorY. Sugita and
  authorY. Okamoto,
  journalChem. Phys. Lett. volume314,
  pages141 (year1999).

[Eastman and Pande(2010)]Eastman2010
authorP. Eastman and
  authorV. S. Pande,
  journalComput. Sci. Eng. volume12,
  pages34 (year2010).

[Eastman et al.(2013)Eastman, Friedrichs,
  and Chodera]Eastman2013
authorP. Eastman,
  authorM. Friedrichs,
  and authorJ. Chodera,
  journalJ. Chem. Theory Comput. volume9,
  pages461 (year2013).

[Pattabiraman et al.(1985)Pattabiraman,
  Levitt, Ferrin, and Langridge]Pattabiraman1985
authorN. Pattabiraman,
  authorM. Levitt,
  authorT. E. Ferrin,
  and
  authorR. Langridge,
  journalJ. Comput. Chem. volume6,
  pages432 (year1985).

[Meng et al.(1992)Meng, Shoichet, and
  Kuntz]Meng1992
authorE. C. Meng,
  authorB. K. Shoichet,
  and authorI. D. Kuntz,
  journalJ. Comput. Chem. volume13,
  pages505 (year1992).

[Onufriev et al.(2004)Onufriev, Bashford,
  and Case]Onufriev2004
authorA. Onufriev,
  authorD. Bashford, and
  authorD. A. Case,
  journalProteins: Struct., Funct., Bioinf.
  volume55, pages383 (year2004).

[Shirts and Chodera(2008)]Shirts2008
authorM. R. Shirts and
  authorJ. D. Chodera,
  journalJ. Chem. Phys. volume129,
  pages124105 (year2008).

[Wood et al.(1991)Wood, Muhlbauer, and
  Thompson]Wood1991
authorR. H. Wood,
  authorW. C. F. Muhlbauer,
  and authorP. T.
  Thompson, journalJ. Phys. Chem.
  volume95, pages6670 (year1991).

[Zuckerman and Woolf(2004)]Zuckerman2004
authorD. M. Zuckerman
  and authorT. B. Woolf,
  journalJ. Stat. Phys. volume114,
  pages1303 (year2004).

[Nunes-Alves and Arantes(2014)]Nunes-Alves2014
authorA. Nunes-Alves and
  authorG. M. Arantes,
  journalJ. Chem. Inf. Model. volume54,
  pages2309 (year2014).

[Swets et al.(2000)Swets, Dawes, and
  Monahan]Swets2000
authorJ. A. Swets,
  authorR. M. Dawes, and
  authorJ. Monahan,
  journalSci. Am. volume283,
  pages82 (year2000).

[Mysinger and Shoichet(2010)]Mysinger2010
authorM. M. Mysinger and
  authorB. K. Shoichet,
  journalJ. Chem. Inf. Model. volume50,
  pages1561 (year2010).

[Pordes et al.(2007)Pordes, Petravick,
  Kramer, Olson, Livny, Roy, Avery, Blackburn, Wenaus, Würthwein
  et al.]Pordes2007
authorR. Pordes,
  authorD. Petravick,
  authorB. Kramer,
  authorD. Olson,
  authorM. Livny,
  authorA. Roy,
  authorP. Avery,
  authorK. Blackburn,
  authorT. Wenaus,
  authorF. Würthwein,
  et al., journalJ. Phys.: Conf. Ser.
  volume78, pages012057
  (year2007).






























