{'fraction_non_alphanumeric': 0.03344780934540535, 'fraction_numerical': 0.015069205952752627, 'mean_word_length': 4.825876076148902, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Protein-Ligand Binding Potential of Mean Force Calculations with Hamiltonian Replica Exchange on Alchemical Interaction Grids

A binding potential of mean force (BPMF) is a free energy of noncovalent association in which one binding partner is flexible and the other is rigid. I have developed a method to calculate BPMFs for protein-ligand systems. The method is based on replica exchange sampling from multiple thermodynamic states at different temperatures and protein-ligand interaction strengths. Protein-ligand interactions are represented by interpolating precomputed electrostatic and van der Waals grids. Using a simple estimator for thermodynamic length, thermodynamic states are initialized at approximately equal intervals. The method is demonstrated on the Astex diverse set, a database of 85 protein-ligand complexes relevant to pharmacy or agriculture. Fifteen independent simulations of each complex were started using poses from crystallography, docking, or the lowest-energy pose observed in the other simulations. Benchmark simulations completed within three days on a single processor. Overall, protocols initialized using the thermodynamic length estimator were system-specific, robust, and led to approximately even replica exchange acceptance probabilities between neighboring states. In most systems, the standard deviation of the BPMF converges to within 5 kT. Even with low variance, however, the mean BPMF was sometimes dependent on starting conditions, implying inadequate sampling. Within the thermodynamic cycle, free energies estimated based on multiple intermediate states were more precise, and those estimated by single-step perturbation were less precise. The results demonstrate that the method is promising, but that ligand pose sampling and phase space overlap can sometimes prevent precise BPMF estimation. The software used to perform these calculations, Alchemical Grid Dock (AlGDock), is available under the open-source MIT license at https://github.com/ccbatiit/algdock/.


INTRODUCTION
Fast and accurate predictions of binding free energies between proteins and small organic ligands would have significant impact on designing drugs [1][2][3] and other modulators of biological processes. The clear relevance of protein-ligand binding affinity prediction in chemical biology and drug discovery has inspired a vast array of physics-based methods (for a broad review, see [4]), each with a different trade-off between computational accuracy and speed.
On one extreme, molecular docking focuses on speed. Docking algorithms are designed to quickly obtain plausible configurations of a protein-ligand complex. Scoring functions are then used to rank one configuration versus another. Docking programs are commonly assessed by their ability to redock ligands into crystallographic structures from which they have been removed. Comparative studies [5,6] and blinded exercises [7] consistently show that docking methods are adept at generating the native pose but are less competent at giving it the highest rank. In this context, it is not surprising that docking scores are poorly correlated with binding free energies [5,6,8,9].
In contrast, alchemical pathway methods are based on rigorous statistical mechanics. The methods involve sampling from a series of possibly nonphysical thermodynamic states in between end-states where the receptor and ligand are bound and unbound. In the unbound state, the receptor-ligand nonbonded interaction terms may be switched off or the species may be physically separated. In accordance to the established statistical mechanics of noncovalent binding [10], the receptor is usually allowed full flexibility. Unfortunately, sampling a fully flexible complex from multiple statistical distributions along an alchemical pathway generally requires substantial computing resources; it has been suggested that most published studies are not fully converged [11]! Even sampling the ligand binding pose is challenging [11]. To bypass this difficulty, most pathway calculations pursue relative binding free energies between similar molecules and assume that the binding mode does not change. In absolute binding free energy calculations, ligand sampling issues are usually alleviated by confining the molecule to a specific pose [12,13]. Within this restricted range of problems, alchemical pathway calculations are amassing a growing track record of accurate prediction (e.g. [14][15][16][17][18][19][20][21][22]). The success of these methods suggests that docking may be substantially improved by incorporating rigorous statistical mechanics.
Implicit ligand theory (ILT) [23], a recently derived statistical mechanics framework for noncovalent association, has the potential to inspire new methods that combine the speed of docking and the rigor of alchemical pathway methods. ILT formally separates receptor and ligand sampling into two distinct stages. Because the first stage of receptor sampling does not require a lig-and, receptor conformations can be sampled once and used with many different ligands. In contrast, conventional alchemical pathway methods require thorough receptor sampling for every receptor-ligand pair. The second stage is to calculate the binding potential of mean force (BPMF) -the binding free energy between a ligand and a rigid receptor configuration (Eq. 2) -for each receptor configuration. Finally, the standard binding free energy is an exponential average of BPMFs (Eq. 3).
While multiple BPMFs are required to estimate a standard binding free energy, the overall calculation should require less computer time than conventional methods because BPMFs are easier to estimate than binding free energies with a fully flexible receptor. The key reason for this speedup is that a BPMF calculation only requires sampling of the ligand, which usually has many fewer degrees of freedom than the complex. Furthermore, nonbonded interactions between a rigid receptor and ligand can be treated by interpolating precomputed threedimensional grids, a strategy first developed for docking [24]. Once the grid is stored, calculation time no longer depends on the size of the receptor. In contrast, conventional alchemical pathway methods require frequent force evaluation between flexible receptor atoms. As the number of pairwise interactions scales as O(N 2 ) with N receptor atoms (neglecting cutoffs), the relative efficiency of ILT-based methods will be more pronounced as receptor size increases.
In the first paper on ILT, BPMFs were estimated for simple host-guest systems [23]. Here, the focus is on precise BPMF estimation for protein-ligand systems. Mobley et al. [14], Ucisik et al. [25] previously computed binding free energies between simple ligands and a rigid protein, T4 lysozyme. In contrast, the present work involves more diverse systems. Mobley et al. [14] were equally rigorous, but did not implement specialized methods to significantly speed their calculations compared to flexible-receptor calculations. Ucisik et al. [25] developed a fast method based on a number of approximations. As in other BPMF [23] and standard binding free energy [22,26,27] calculations, Hamiltonian replica exchange is applied. The main methodological differences between the current and previous work are the adaptive initialization of thermodynamic states and the use of linearlyscaled interaction grids [24] with transformation-based smoothing [28]; these will be discussed in detail in the section on Theory and Methods. The method is demonstrated on the Astex diverse set [29], a curated database of 85 high-quality protein-ligand complexes of pharmaceutical or agrochemical interest.
The algorithm is implemented in a new software package, Alchemical Grid Dock (AlGDock), a python module based on the Molecular Modeling Toolkit (MMTK) 2.7.8 [30]. AlGDock is available under the open-source MIT license at https://github.com/ccbatiit/algdock/.

THEORY AND METHODS
This section reviews implicit ligand theory, details the algorithms in AlGDock, and describes the setup of the BPMF calculations on the Astex diverse set.

Implicit Ligand Theory
For the noncovalent association between a receptor R and ligand L to form a complex RL, R + L RL, the standard binding free energy is, where β = (k B T ) −1 is the inverse of Boltzmann's constant times the temperature, C • is the standard concentration (1 M = 1/1660Å 3 ), and C X is the equilibrium concentration of species X ∈ {R, L, RL}. Activities have been assumed to be unity, a reasonable approximation in the limit of low concentrations.
Coordinates of the complex, r RL , are partitioned into receptor (r R ) and ligand internal (r L ) and external (ξ L ) coordinates. Based on this partitioning, the interaction energy is defined as Ψ is an effective potential energy that includes the gas-phase potential energy U (r) and solvation free energy W (r) [10]. A BPMF is an exponential average of interaction energies over ligand coordinates in the binding site [23], where I ξ ≡ I(ξ L ) is an indicator function that takes values between 0 and 1 and specifies whether the receptor and ligand are bound or not. According to ILT, the standard binding free energy ∆G • is related to an exponential average of BPMFs over Boltzmann-distributed receptor configurations r R , where Ω = I ξ (ξ L )dξ L is the volume of the binding site.

Thermodynamic Cycle
A BPMF can also be expressed as a ratio of partition functions [23], In AlGDock, BPMFs are estimated by completing the thermodynamic cycle shown in Figure 1. Mathematical expressions for individual free energy differences within the cycle are given in Table I. The sum of all the free energy components in Table I is βB(r R ).
a&b Ligand desolvation f ab,L = − ln b&c Ligand heating f bc,L = − ln Expressions for free energy differences from the thermodynamic cycle in Figure 1. The target temperature TT = 300 K is the basis of βT = (kBTT ) −1 . Similarly, the high temperature TH = 600 K is the basis of βH = (kBTH ) −1 .
Over the course of the cycle, the receptor-ligand interaction strength is scaled and the temperature is varied. The temperature is varied because high-temperature states enhance transitions between local energetic minima. Because the protocol involves temperature changes, the notation is simplified by using a lower-case u to refer to the reduced potential energy [31], a log probability density that incorporates β. For example, at milestone c, the reduced potential energy is u(r RL ) = β H [U (r L ) + U (r R )]. Likewise, the reduced free energy difference between two milestones x and y is f xy . Be-cause converting from reduced to standard potential energies and free energies involves dividing by β, the units of these reduced quantities are specified as k B T .
For thermodynamic states that require I ξ , the ligand is confined to the binding site using a flat-bottom harmonic potential [23,26,27], where k = 10000 kJ/(mol nm 2 ) is the spring constant, d is the distance between the ligand center of mass and the center of the binding site, and d 0 = 6.0Å is the radius of the binding site. There is no restriction on ligand rotation. At milestone d, the reduced potential energy is, where Ψ g (r RL ) is the grid interaction energy. The grid interaction energy, is based on one electrostatic and two van der Waals grids. The electrostatic interaction energy Ψ P BSA (r RL ) is evaluated by multiplying atomic partial charges with the electrostatic potential. The electrostatic potential at each ligand position is obtained by trilinear interpolation of a precomputed grid. The grid is produced by solving the linear Poisson-Boltzmann equation around the minimized receptor molecule using APBS 1.4 [32] with sequential focusing. Coarse grids are at least 1.5 times larger than the range of the receptor molecule in each dimension. Fine grids have the same size as the van der Waals grids, and a spacing of 0.5Å. Coarse grids use multiple Debye-Huckel boundary conditions, and fine grids use coarse-grid solutions as boundary conditions. Both grids are solved with the following options: a quintic Bspline charge discretization, spline window width of 0.3, protein dielectric of 2.0, solvent dielectric of 80.0, solvent density of 10.0, solvent radius of 1.4Å, smoothed dielectric and ion-accessibility coefficients, and temperature of 300.0 K.
The van der Waals energy Ψ vdW (r RL ) is evaluated by an analogous grid-based procedure pioneered by Meng et al. [24]. In the AMBER [33] molecular mechanics force field, the receptor-ligand van der Waals interaction energy is represented as, Aij rij 12 − Bij rij 6 , a double sum over ligand atoms i ∈ 1, ..., N lig and receptor atoms j ∈ 1, ..., N rec . A ij and B ij are the repulsion and attraction parameters and r ij is the distance between atoms i and j. Molecular docking programs usually model these interactions by a geometric mean approximation [24], such that With this approximation, the receptor-ligand van der Waals interaction energy at a point k is, The values a k and b k are precomputed on a threedimensional grid with a spacing of 0.25Å. Grids span at least 10 + 1.5d maxÅ in each dimension surrounding the ligand center of mass in the crystallographic binding mode, where d max is the maximum distance from any ligand atom to its center of mass. For energy and force evaluations, a k and b k are approximated at actual ligand atom positions by trilinear interpolation. Because van der Waals potentials are highly nonlinear, straightforward trilinear interpolation of van der Waals grids provides energies that poorly match directly calculated energies [34]. Hence, van der Waals repulsive grid energies are calculated using a transformation, trilinear interpolation, and inverse transformation procedure [28]: If f (x i ) is the van der Waals attractive or repulsive energy at x i , a position vector on a grid corner, then g(x i ) = f (x i ) −1/m is evaluated for all grid corners. To estimate f (x) at a desired position vector x, g(x) is obtained by trilinear interpolation of g(x i ) for the eight x i surrounding x, and f (x) = g(x) −m . Previously, when equivalent powers were used for the van der Waals repulsive and attractive grids, m = 2 was found to dramatically reduce interpolation error [28]. Here, m = 2 is used for the repulsive grid and no transformation for the attractive grid; this combination was found to yield the most accurate van der Waals energies for docked poses.
Alchemical transformations that modulate the strength of interactions between atoms often face an "end-point catastrophe" where free energy changes are numerically unstable [1]. To circumvent this issue, a set of soft Lennard-Jones repulsive and electrostatic grids is introduced between milestones c&d. In these soft grids, the original grid value, v o , is replaced with v max tanh (v o /v max ). (Gallicchio and Levy [26] also used a hyperbolic tangent energy cap.) For the soft Lennard-Jones repulsive grid, v max = 10.0 kJ 1/2 mol −1/2 . To prevent ligands from being "pinned" by electrostatics, the soft electrostatic grid uses as maximum value such that the electrostatic energy is less than or equal to the soft Lennard-Jones repulsive energy for every heavy atom at every grid point. This is established by setting v max for the electrostatic grid to 10 times the minimum ratio of Lennard-Jones and electrostatic scaling factors. The reduced potential energy is switched according to the protocol, This protocol turns on the soft grids first, and then the unperturbed grids ( Figure 2). The potential is consistent with milestone c at α = 0 and milestone d at α = 1.

Sampling and Estimation
The sampling and estimation strategy in AlGDock is summarized as follows: 1. Ligand preparation: the ligand is minimized with 500 conjugate gradient steps. The temperature is ramped from 20 K to 300 K over 30 geometrically spaced simulations of 500 steps each.
2. Thermodynamic state initialization: starting from 50 seed configurations, simulations of 1000 steps were used to initialize each thermodynamic state between milestones b&c and c&d.
3. Hamiltonian replica exchange: simulations were conducted between milestones b&c and subsequently between c&d.
4. Postprocessing: samples from milestones b and d were postprocessed using the generalized Born/surface area (GBSA) water model [35]. For milestone e, grids were replaced by the rigid receptor configuration and pairwise interactions were directly computed.

5.
Estimation: Free energy differences that sum up to the BPMF are estimated.
During thermodynamic state initialization and Hamiltonian replica exchange, the No-U-Turn Sampler (NUTS) [36] is used to propagate the ligand from one configuration to the next. NUTS is a hybrid molecular dynamics (MD)/Markov chain Monte Carlo (MCMC) sampling technique [37] that eliminates the need to select the MD trajectory length between Monte Carlo trials. In other hybrid MD/MCMC techniques, the trajectory length is an adjustable parameter that must be carefully tuned. NUTS automatically truncates MD trajectories when they make a "U-Turn" that brings trajectory end points closer to one another, or when energy changes are too large. In the present implementation, the change in potential energy or total energy (potential and kinetic energy) is required to be less than 125 kJ/mol. The implementation also truncates trajectories so that the allotted number of MD integration steps is not exceeded.
Another adjustable parameter, the MD time step, is selected using dual averaging -an extension to NUTS [36]. Dual averaging adjusts the time step such that an acceptance statistic H t becomes close to δ. Parameters for dual averaging were δ = 0.6, γ = 0.05, t 0 = 10, κ = 0.75, and µ = ln(10∆t). Because adjusting the time step disturbs the equilibrium distribution, dual averaging is only used during initialization.
To accelerate transitions between binding poses, external coordinate Markov chain Monte Carlo moves are attempted for thermodynamic states between milestones c&d for α < 0.01. The external coordinate move consists of: 1. Random rotation. A random quaternion is converted into a matrix that is used to rotate the molecule about its center of mass.
2. Random translation. The magnitude of translation in each dimension is drawn from a Gaussian distribution with a standard deviation of 0.6Å.
The move is accepted or rejected according to the Metropolis criterion. Moves are not attempted for α > 0.01 due to low acceptance probability.

Thermodynamic state initialization
For milestones b&c, the first thermodynamic state (k = 0) is at 300 K and the ligand is steadily warmed to 600 K. The first thermodynamic state between milestones c&d depends on whether there is a fully-bound pose available in the binding site. If a pose is available, then the first state is fully bound at 300 K. Otherwise, the first state is fully unbound at 600 K. In the latter situation, 50 randomly selected configurations from milestone c are placed in the binding site at 453 random center-ofmass positions (0.5 positions perÅ 3 ) and rotated with 100 different random orientations.
After state k is initialized, state k + 1 is initialized as follows: 1. Parameter selection: Parameters for state k + 1 are selected using a new algorithm designed to separate states at approximately even intervals in thermodynamic length.
Thermodynamic length is a metric of the distance on the manifold of thermodynamic states [38]. For a sequence of states, the statistical error in free energy calculations is minimized and the replica exchange frequency is nearly maximized when the intermediate states are equidistant in thermodynamic length [39]. Suppose that thermodynamic parameters (e.g. temperature, pressure, grid scaling strengths) are specified by a vector λ with components λ i . Let γ ≡ γ(α) describe the dependence of λ on the variable α, such that γ(0) is the initial and γ(1) is the final thermodynamic state. For microscopic systems, the thermodynamic length is defined by the path integral [39,40], Given a parameter vector λ, the reduced potential is the effective potential energy and T λ is the temperature. The normalized log probability of observing a configuration x is l λ ( dx is the partition function. These quantities are used to define elements of the Fisher information matrix, where σ 2 λ is the covariance in state λ and ∂ i denotes a partial derivative with respect to λ i . For a protocol in which only one parameter λ i varies with α, the length is, L = 1 0 Numerical estimates of L are most accurate when samples are drawn from many states between 0 < α < 1 [40]. Such exhaustive sampling, however, is unavailable during initialization. A simple approximation for the thermodynamic length when one parameter changes is, L = ∆λ i σ 0 ∂ i l λ , where ∆λ i is the total change in the value of the parameter λ i and σ 0 ∂ i l λ is a standard deviation in the initial state. Thus, if one desires L to be approximately constant between different intermediate stages in a protocol, then the change in parameter should be inversely proportional to σ 0 ∂ i l λ , where s is an adjustable parameter, the thermodynamic speed.
Between milestones b&c, the parameter that varies is the temperature T. As the log probability of a ligand configuration is where s bc = 0.1. Between milestones c&d, the log probability of r RL is l with s cd = 0.2.
If the targeted value of the parameter is exceeded (e.g. temperature increased above 600 K), then the targeted value is used.
2. Seed selection: 50 configurations from state k are resampled as starting seeds for simulations in state k + 1.
Configurations are drawn from state k with weights is the reduced potential energy in state k. In the limit of infinite sampling of state k, resampled configurations would be Boltzmann-distributed in state k + 1. With imperfect sampling of state k, resampled configurations approximate the Boltzmann distribution in state k + 1.
3. Sampling: Simulations of 1000 steps are run from each seed.
NUTS is run with dual averaging to adjust the time step. Each simulation of 1000 steps has its own estimate of the ideal time step. The median time step, restricted to between 0.25 and 2.5 fs, is subsequently used in Hamiltonian replica exchange.

Verification:
The mean replica exchange probability, p acc , is used to verify that thermodynamic states are not too distinct nor similar.
It is estimated by taking the sample mean of p acc (Eq. 17) for every pair of initial samples (at the same time index) from states k and k + 1. If p acc is estimated to be too low (below 0.3), then parameters for state k + 1 are reselected with a smaller increment (thermodynamic speed is adjusted by a factor of 4/5) and simulations are repeated. If it is too high (above 0.99), then state k is removed.

Hamiltonian replica exchange
Initialization is followed by Hamiltonian replica exchange [23,26,27,41] between milestones b&c (15 cycles) and c&d (18 cycles). Replica exchange is a Markov chain Monte Carlo move that swaps the configurations of a pair of simulations at different thermodynamic states. (Equivalently, it may be regarded as swapping the states.) Consider the thermodynamic states a and b with reduced energies u a and u b , respectively. If x is the original configuration in state a and y the original configuration in state b, then the acceptance probability, preserves the Boltzmann distribution in both states. Typical replica exchange protocols attempt exchanges between pairs of neighboring thermodynamic states, but this restriction is unnecessary. As replica exchange is a type of Gibbs sampling [42], an arbitrary number of attempts can be made between arbitrary pairs of states. In the present simulations, each replica exchange cycle consists of 1000 iterations with 50 NUTS steps, 20 external coordinate MCMC moves (if appropriate), and a sweep of replica exchange. Each sweep of replica exchange includes attempts to swap configurations between pairs of states that are 1, 2, ..., min(5, K) states apart, where K is the total number of thermodynamic states in the direction.
The statistical inefficiency g is estimated from the integrated autocorrelation time of replica exchange indices, as described [43]. Multiple samples (up to 3 between milestones b&c and 25 between milestones c&d) were saved for every g replica exchange sweeps.

Estimation
BPMFs were estimated according to, This expression replaces the sum f bc,R + f cd,o with, which can be evaluated without determining the receptor internal energy U (r R ). The receptor desolvation free energy is estimated by the difference, f ab,R = β T (U (r R ) − U(r R )). Other free energy differences are estimated based on equilibrated samples from replica exchange between milestones b&c or c&d. Equilibration is determined by calculating the mean and standard deviation of the energy in milestones b or d, respectively. Systems are considered to be equilibrated once the mean energy is within a standard deviation of the mean energy in the last cycle. f ab,R and f de are estimated by free energy perturbation [44] using configurations drawn from milestones a and d, respectively. f bc,L and f cd are estimated by the multistate Bennett acceptance ratio [31], which uses potential energies from every replica.

Astex diverse set BPMF calculations
The Astex diverse set is a subset of the protein data bank carefully curated for quality and diversity [29]. The members of the set will be referred to be their protein data bank identifiers, e.g. 1s3v. Each system in the set was set up using AmberTools 14 [45]. Proteins and ions (from protein.mol2) were parameterized using the AMBER ff14SB force field. Zinc parameters were also used for iron, mercury, and magnesium. Non-standard amino acid residues were converted to their standard counterparts. Using the python wrapper to the Open-Eye OEChem Toolkit, all crystallographic ligands were parameterized with the Generalized Amber Force Field [46] and AM1BCC charges [47,48]. The main ligands from ligand.mol were regenerated to lose memory of the crystallographic pose. Noncrystallographic atoms in each receptor were minimized for 2500 steps in gas and 2500 steps in GBSA solvent using NAMD 2.9 [49].
For each system, 15 independent simulations were performed between milestones b&c. Subsequently, simulations were continued between milestones c&d from three different starting points: crystallography, docking, or the lowest-energy pose (according to u g ) sampled in all the other simulations. In most situations, starting from docked poses will be the most practical approach, but the others are a useful comparison.
Docked poses were obtained using the anchor and grow algorithm in UCSF DOCK 6 [50] to redock ligands into their respective proteins. First, sphgen was run with default parameters: dotlim = 0.0, radmax = 4.0Å, and radmin = 1.4Å. DOCK was run with 20000 maximum orientations, using internal energy with an exponent of 12, a flexible ligand, an a minimum anchor size of 40. Pruning was performed with clustering, 1000 maximum orientations, a clustering cutoff of 1000, and conformer score cutoff of 25.0. A bump filter was used with max bumps anchor = 12 and max bumps growth = 12. Final conformations were clustered with a root mean square deviation (RMSD) threshold of 2.0Å. For two systems, 1t46 and 1n46, no docked poses in the binding site were obtained.
After the starting poses were minimized for 1000 conjugate gradient steps in u g , the lowest-energy pose was used to initialize milestone d. After the thermodynamic states were initialized, docked poses within the binding site were also used as starting points for each state in replica exchange. The lowest-energy pose was used in milestone d. Higher-energy poses were used in intermediate replicas to fill all available thermodynamic states. If there were more states than docked poses, the lowest-energy pose was duplicated. For simulations starting from a single pose, replicas were not reinitialized at the beginning of replica exchange.
Calculations were performed on the Open Science Grid [51], Duke shared computing resources, or the Minh group computing cluster at IIT. Each node on the Minh group cluster includes dual AMD Opteron 6376 2.30GHz 16MB cache sixteen-core processors and 64GB DDR3 1600 ECC registered memory. Benchmark calculations were run on the Minh group cluster.
To quantify the precision of the results, we consider the root second moment about a reference value x r , for N samples indexed by x n . When x r is the mean, then Equation 20 is the standard deviation, which will simply be referred to as σ[x]. We also consider the root second moment about the minimum, σ[x, x min ], and about the sample mean at the end of a simulation, σ[x,x end ].

Protocol length
Overall, the present method for thermodynamic state initialization robustly yields protocols that are adapted to specific systems. The large range in the number of thermodynamic states, N states , provides evidence of system-specific adaptation (Figure 3). Between milestones b&c, the average number of states,N states ranges from 35.1 to 93.6. Between milestones c&d for simulations initialized from the crystallographic pose, thē N states ranges from 84.6 to 190.1. For most systems, the standard deviation of the number of states, σ[N states ] is small relative to theN states , indicating that the procedure is robust. The robustness and apparent systemspecificity of the protocols implies that a universal thermodynamic protocol may strongly deviate from optimality, e.g. p acc may be inconsistent across replicas.
In some systems, however, N states is not robust; N states is dependent on the starting pose and/or σ[N states ] is large. Protocols will vary when the initialization stage samples from distinct local minima and estimates different thermodynamic lengths. Between milestones b&c, the σ[N states ] is larger than 2 in three sets of simulations: 1tz8 (6.1), 1r1h (7.3), and 1vcj (13.6). The larger variance in protocol lengths occurs for different reasons (Figure 4). In 1r1h simulations, the protocols require variable number of states to reach 400 K, after which the slope is similar. For 1tz8 and 1vcj, there are two classes of protocols, indicative of sampling from distinct local minima. With 1tz8, the two protocols are obtained with approximately equal probability. With 1vcj, the short protocol was observed in only 2 of 15 simulations. Protocols between milestones c&d are more variable than those between milestones b&c. Interaction grids introduce more possibilities for trapping in local minima. For all three starting points, 1l7f has the largest variance in N states . The protocols are most similar for the first 40 thermodynamic states and the last 10, but take different paths in between ( Figure 5). It is worth noting, however, that high variation in protocols does not necessarily lead to inaccurate results, as will be explained in the following sections.

Time steps
In the vast majority of initialization processes, dual averaging converged to a time step around 1.5 fs ( Figure  6 and Table II). This is true both between milestones b&c and c&d, suggesting that the grid does not significantly affect the time step in the vast majority of cases. As a notable exception, the shortest time steps appear around α = 0.5 between milestones c&d in simulations of 2bm2 and 1n46. In 1n46, there is a steady decrease in time step away from the end points. Because simulations of 1n46 did not start with a docked pose, they are more likely to occupy metastable minima in the grid. Near α = 0, the grid has little influence. Near α = 1, simulations are dominated by a smaller number of more stable poses. The longest time steps are observed with 1jd0, where most simulations converged to time steps between 1.8 and 2 fs. The ligand in 1jd0 is notable for being small, relatively rigid, and possessing few hydrogen atoms (5). The small number of hydrogen atoms is relevant because unconstrained hydrogen vibrations often limit the maximal time step of molecular dynamics integrators.

Replica exchange acceptance probabilities
A good replica exchange protocol has a reasonable p acc between all neighboring states. The key benefit of replica exchange is to spread sampled configurations across a range of different thermodynamic states. A bottleneck in the exchange of configurations across the pair of states can eliminate this benefit of replica exchange; the simulations effectively become two independent sets of simulations. Low exchange probabilities are also indicative of poor phase space overlap, which can limit the convergence of free energy estimates [52].
Because of its importance to efficient simulation, Al-GDock includes an estimate of p acc to verify new thermodynamic states during the initialization step. (The notationp acc is used to refer to an estimator for p acc ). In some simulations, however, the phase space explored during replica exchange can be distinct from that explored during initialization, leading to a substantial change in the observed p acc . Thus, thermodynamic state verification was followed up by estimating p acc based on equi- Estimates of p acc indicate that the initialization protocol does not lead to any simulations with replica exchange bottlenecks (Figure 7 and Table III). Between milestones b&c,p acc estimated during replica exchange are high and have low variance. This outcome is consistent with achieving the goal of nearly equal thermodynamic length between adjacent thermodynamic states. It also implies that the configuration spaces explored during initialization and replica exchange are largely the same. Between milestones c&d, the replica exchange rates are also high but there is a larger variance. While most p acc are between 0.7 and 1.0, there are a few simulations wherep acc is much lower; the lowest observed values are during simulations of 2bm2 (0.21) and 1gm8 (0.28). In these simulations, the acceptance probability dips around α = 0.75, but are high for most other values of α. These drops indicate that different conformations are explored in replica exchange compared to during thermodynamic state initialization. However, even these lowp acc are not low enough to be considered a bottleneck; they are expected to allow configurations to pass through the ther-  The observedp acc also underscore the point that large variation in protocols is not necessarily problematic. While simulations of 1tz8, 1r1h, and 1vcj have a relatively large variation in N states between milestones b&c, allp acc estimated from replica exchange are greater than 0.9 ( Figure 8). As is particularly clear for low temperature states of 1r1h, longer protocols have a higherp acc . Another bump inp acc occurs between the last pair of thermodynamic states, where an even interval in ther- modynamic length would necessitate going higher than the temperature of interest.

Simulation time
Due to the heterogeneity of computing resources, variability of queue times, and the diverse nature of the systems, calculations took a variable length of time to complete. Some calculations completed in one day, and all were done within a week.
Total CPU time of benchmark simulations spanned a large range, from 13.6 to 71.1 hours (Figure 9). In all simulations, the majority of time was spent in replica exchange between milestones c&d, and the second largest fraction in replica exchange between milestones b&c. Initialization and analysis between milestones c&d consumes a small but noticeable fraction of simulation time, whereas the analogous steps between milestones b&c consume a nearly negligible fraction of simulation time.
There is a trend relating the calculation time affiliated with milestones a to c and milestones c to e, but the correlation is imperfect.

Sampling
At milestone c, sampled configurations are uniformly distributed in a sphere. As α increases from 0 to 1 between milestones c&d, the phase space of the ligand is gradually restricted. First, the soft grids prevent ligand atoms from overlapping with receptor atoms. At intermediate values of α, the ligand may assume several poses (e.g. Figure 10). Finally, at milestone d, the ligand usually but does not always sample from a single minimum.
In many simulations, the phase space sampled at higher α is a subset of that sampled at lower α (e.g. Figure 11). In others, however, the conformational minima for α ≈ 0.5, where the soft grids are at full strength, are entirely distinct from the minima for α ≈ 1 (e.g. Figure 10). When there are shifts in important phase space, the thermodynamic states at α ≈ 0.5 are less beneficial to sampling from milestone d. Nonetheless, phase space overlap between adjacent states is a sufficient condition for precise free energy estimates. At milestone d, most simulations sample from a single minimum (Figure 12), but there are several other situations. These include, 1. Sampling several alternate poses that are (e.g. Figure 12(a)) or are not (e.g. Figure 12(d)) clearly separated; 2. Sampling a clearly defined but sparsely populated alternative pose (e.g. Figure 12(b)); 3. Sampling alternate poses that share a common warhead position, but have a floppier tail (e.g. Figure 12(g)).

Native pose identification
A large fraction of simulations include the crystallographic pose among samples from milestone d ( Figure  13). Unsurprisingly, it is most often observed in simulations starting from the crystallographic pose. After equilibration, all the simulations still sample a pose with a heavy-atom RMSD less than 0.75Å. In contrast, a native pose (defined as having a heavy-atom RMSD less than 2.0Å) is observed in 85.6% of simulations starting from the lowest energy pose and 77.7% of simulations starting from docked poses. According to u g , the lowest-energy observed pose is a native pose in 61/85 = 71.8% of systems. This implies that Boltzmann sampling was able to start at another pose and find the native pose in 13.8% of simulations. Starting docked poses include a native pose in 75/85 = 88.2%, and rank the native pose with the lowest u g (and lowest UCSF DOCK 6 grid score) in 48/85 = 56.5% of systems. This implies that the sampling procedure was able to bring the native pose into milestone d, via Monte Carlo moves or replica exchange, in 21.2% of simulations.
Overall, the GBSA force field is better at identifying the native pose than the grid energy u g . According to the grid energy, the lowest energy configuration has a RMSD less than 2.0Å in 78.5% of simulation starting from the crystallographic pose, 68.9% of simulations starting from the lowest energy pose, and 55.6% of simulations starting from docked poses. According to the GBSA force field, the analogous fractions are 87.1%, 75.1%, and 61.5% of simulations, respectively.

Standard deviation of free energy estimates
In most systems, the standard deviation of the BPMF converges to within 5 k B T ( Figure 14 and Table IV). Within the thermodynamic cycle, free energies estimated based on multiple intermediate states (f bc and f cd ) were more precise, and those estimated by single-step perturbation (f ab and f de ) were less precise. Simulating from multiple intermediate states promotes phase space overlap between adjacent states, leading to more precise free FIG. 11. Samples from evenly spaced thermodynamic states between milestones c&d, taken from a representative simulation of 1kzk starting from docked poses. The protein structure structure is shown with ribbons and the crystallographic ligand pose is shown with a thick licorice representation and purple carbon atoms. The same illustration scheme is used in Figures 10 and 12. Figures 11, 10, and 12 were generated with VMD [53].  (Table V) to improve the precision of free energy estimates. While σ[f, f end ] usually decreases monotonically with sample size, this trend is not universally true ( Figure 15). In some systems there are temporary jumps in σ[f, f end ] due to instability in determining the equilibrated cycle. More importantly, in other systems, σ[f, f end ] appears flat, with little or no change with increased sampling.
In some systems, the variance is low but free energy estimates are dependent on the starting poses. When results from all simulations between milestones c&d are grouped together, the standard deviation in many systems no longer converges to within 5 k B T . Because converged thermodynamic quantities should not depend on initial conditions, this inconsistency implies that some simulations are not fully converged.
When free energy estimates contradict, it is not completely clear which result is most consistent with the force field. False convergence is quantified based on the minimum value, using σ[x, x min ]. A system is considered to be falsely converged when σ[f ] is less than 5 k B T but σ[f, f min ] is greater than 15 k B T . Based on this definition, false convergence occurs in a nontrivial fraction of systems (Table IV).
Ultimately, poor convergence and false convergence are caused by incomplete sampling. The most straightforward explanation for incomplete sampling is that a simulation starts and remains trapped in a local min-imum distinct from the global minimum energy. Indeed, between milestones c&d, most falsely converged calculations starting from the crystal structure (4/5) and docked poses (7/11) have no starting pose within 1.0 A of the lowest-energy observed pose. However, sampling the lowest-energy structure is not a sufficient condition for complete sampling, as other poses may have a lower free energy.

DISCUSSION
The selection of intermediate states between two thermodynamic milestones of interest is a ubiquitous problem in molecular simulation. The usual approach is an iterative trial-and-error process starting from a naive protocol, checking for issues such as replica exchange bottlenecks, and inserting and removing states as necessary. I have developed a simple and robust approach to initialize a series of thermodynamic states based on only a single adjustable parameter, the thermodynamic speed. Due to this automated procedure, I was able to run a large number of simulations on a diverse array of protein-ligand complexes. In the vast majority of cases, the protocols yielded consistent replica exchange rates across neighboring thermodynamic states without further fine-tuning. The described approach to trailblazing thermodynamic state space may find use in other classes of simulations.
Replica exchange calculations in this present study include more thermodynamic states than most published molecular simulations.
Conventional wisdom about replica exchange is that an optimal number of replicas will maximize efficiency. With too few replicas, there is limited phase space overlap between neighbors and exchange rates are vanishingly small. With too many replicas, metrics of replica exchange efficiency, such as a mean round-trip time, diminish. However, in a recent study involving extensive simulation of several distinct processes, it was found that if there are no bottlenecks, the number of states has little impact on the convergence of free energy estimates (the manuscript is under preparation). Hence, I chose to include a large number of states to allow for the possibility that later sampling will explore different regions of phase space and reduce the exchange rate between neighboring states. While useful, the described thermodynamic state initialization process remains imperfect. Replica exchange is particularly beneficial when the important phase space of a thermodynamic states is a subset of the important phase space of another. The present procedure is limited to the variation of a single thermodynamic parameter between milestones, and can do little to promote the subset relationship. Future improvements could accomodate varying multiple parameters (e.g. separate parameters for the temperture, van der Waals grids, and electrostatic grids) in between thermodynamic states of interest, while maintaining a phase space subset relationship. While a direct comparison is difficult, the described procedure appears less successful than the best docking programs at identifying the native pose. Hartshorn et al. [29] found that with a larger site (10Å), the top-ranked GOLD solution is within 2Å of the experimental binding mode in 80.4% of calculations. In contrast, only 56.5% of the present UCSF DOCK 6 calculations were successful, suggesting that setup and parameterization, and perhaps the docking algorithm, could be fine-tuned. Native pose identification based on AlGDock performs similarly to UCSF DOCK 6, suggesting that replica exchange sampling from the Boltzmann distribution does not significantly increase the diversity of poses. As evidenced by the variable performance as a function of starting pose(s), improving Boltzmann sampling schemes should improve native pose identification.
Nonetheless, the sampling and estimation approach is able to attain precise BPMF estimates in the majority of systems. As for the others, it is unsurprising that the same approach does not work equally well in all cases. With a diverse set, ligands and receptors will exhibit distinct levels of complexity and may be driven to bind by different factors. Indeed, the inadequate sampling that led to convergence issues for several systems is probably not unique to BPMF calculations [11]. Because of the relative simplicity of the calculations, I have been able to conduct more replicates from more distinct starting points than other alchemical standard binding free energy calculations. The sheer number of simulations has allowed us to expose problems that may not be clear from fewer calculations. Our results suggest that Hamiltonian replica exchange with molecular dynamics is not always an adequate approach for sampling distinct ligand binding poses. Specialized sampling methods will likely be required in a more broadly effective protein-ligand BPMF estimation procedure. It is also clear that for some systems, there may be limited phase space overlap between different implicit solvent models. More precise BPMF estimation can be attained by introducing intermediate states or implementing more similar force fields for milestones a&b and d&e.

CONCLUSIONS
I have developed a method to initialize thermodynamic states and estimate BPMFs for protein-ligand systems.  The least precise free energy estimates between adjacent milestones. Counts are based on all simulations or the imprecise subset, when the standard deviation of the total BPMF is greater than 5.0.
The largest sources of imprecision are found to be ligand pose sampling and phase space overlap between implicit solvent models. Intriguingly, I found that starting from the lowest-energy pose is not a sufficient condition for converged BPMF estimates. Future studies may build on the current work to compare BPMFs with molecular docking scores as classifiers of binding. Another useful direction will be to assess the convergence of standard binding free energy estimates upon increased receptor conformational sampling. ACKNOWLEDGMENTS I thank Michael Shirts (University of Virginia), John Chodera (MSKCC), and Trung Hai Nguyen (IIT) for helpful discussions. Peter Eastman (Simbios) and John Chodera assisted with implementing an early version of the code in OpenMM. Michael Sherman (Simbios) made an invaluable suggestion of using precomputed grids. Rob Gardner, Lincoln Bryant, and Balamurugan Desinghu (Open Science Grid) and Tom Milledge (DSCR) provided assistance with computing resources. David Beratan (Duke) was a supportive postdoctoral advisor. OpenEye Scientific Software, Inc. and UCSF provided academic licenses to their software. At the beginning of this project, I was a postdoctoral scholar at Duke, supported by NIH 2P50 GM-067082-06-10. I also spent a month as an OpenMM visiting scholar at Simbios.