Massive-Scale Binding Free Energy Simulations of HIV Integrase Complexes Using Asynchronous Replica Exchange Framework Implemented on the IBM WCG Distributed Network
To perform massive-scale replica exchange molecular dynamics (REMD) simulations for calculating binding free energies of protein-ligand complexes, we implemented the asynchronous replica exchange (AsyncRE) framework of the binding energy distribution analysis method (BEDAM) in implicit solvent on the IBM World Community Grid (WCG) and optimized the simulation parameters to reduce the overhead and improve the prediction power of the WCG AsyncRE simulations. We also performed the first massive-scale binding free energy calculations using the WCG distributed computing gird and 301 ligands from the SAMPL4 challenge for large-scale binding free energy predictions of HIV-1 integrase complexes. In total there are ~10 thousand simulated complexes, ~1 million replicas, and ~2000 microseconds of aggregated MD simulations. Running AsyncRE MD simulations on the WCG requires accepting a tradeoff between the number of replicas that can be run (breadth) and the number of full RE cycles that can be completed per replica (depth). As compared with synchronous Replica Exchange (SyncRE) running on tightly coupled clusters like XSEDE, on the WCG many more replicas can be launched simultaneously on heterogeneous distributed hardware, but each full RE cycle requires more overhead. We compared the WCG results with that from AutoDock and more advanced RE simulations including the use of flattening potentials to accelerate sampling of selected degrees of freedom of ligands and/or receptors related to slow dynamics due to high energy barriers. We propose a suitable strategy of RE simulations to refine high throughput docking results which can be matched to corresponding computing resources: from HPC clusters, to small or median-size distributed campus grids, and finally to massivescale computing networks including millions of CPUs like the resources available on the WCG.
Introduction<!>Binding energy distribution analysis method (BEDAM)<!>WCG AsyncRE implementation of BEDAM method<!>Metrics to evaluate the efficiency of the REMD algorithm<!>BCD toy model<!>HIV-1 integrase<!>Optimization of simulation parameters<!>Overhead analysis due to the heterogeneity of WCG simulations<!>Binding free energy predictions of whole SAMPL4 library<!>Binding free energy predictions from different docking poses<!>Conclusion
<p>There are three critical components for accurate binding free energy prediction using molecular dynamics simulations which are of importance for structure-based drug design 1–7 in the early stage of computer aided drug discovery: the statistical theory and computational approximations for binding free energy calculations; the force field functions and parameters for describing the physical systems involved such as the receptor, the ligand, and the solvent; the sampling methods for exploring relevant conformational space. In this report we target the last component, how to combine the replica exchange sampling methods with our binding energy distribution analysis method (BEDAM)8 for calculating absolute binding free energy and optimize the simulation strategy for a distributed resource like the World Community Grid.</p><p>Computational methods for calculating binding free energy1–7 generally perform individual molecular dynamics (MD) simulations at many intermediate thermodynamic states besides the free and fully coupled states. The MD aggregate times, however, are typically limited to the order of microseconds 9–12 even using high performance computing (HPC) resources from XSEDE or specialized CPU/GPU computing units13–15. Developing more advanced conformational sampling methods in the context of generalized ensembles16–31 such as parallel replica exchange (RE) or parallel tempering methods is one way to accelerate the conformational sampling and overcome the timescale challenge due to high free energy barriers resulting in slow dynamics of biomolecular complexes. However, conventional RE methods are implemented with synchronous exchanges18–22,32,33 and are designed for homogeneous environments such as HPC clusters that require the allocation and maintenance of necessary resources for all replicas during the entire simulation and are intolerant to the failure of any individual replica simulation. Those limitations prevent the traditional SyncRE approach from being a feasible solution for new RE application simulations requiring hundreds to thousands of replicas.34–36</p><p>On the other hand the available computing units are not limited to high-end HPC clusters, there exist massively distributed computing units such as the IBM World Community Grid (WCG), a volunteer grid consists of more than 0.7 million members distributed all over the world and 3.0 million heterogeneous computing units including personal or public workstations, laptops and mobile devices, installed with different operating systems such as Linux, Mac OS, and Windows. Those distributed computing grids are highly dynamic and heterogeneous due to the volunteer nature of joined members, the diversity of the computing units, and random pause or termination of running jobs. The implementation of conventional RE methods for those distributed computing grids is difficult and also much less efficient since the slowest computing unit determines the efficiency of the whole RE simulation. There exists previous attempts to develop algorithms better suited for heterogeneous computing grids. Serial tempering (ST)37–39 or simulated tempering only carries out a single thread of an MC/MD simulation in position space, and updates of the thermodynamic state (such as temperature) of the system are performed periodically. ST methods can perform simulations on a single computing unit but requires a pre-estimation of free energy weights at different thermodynamics states and their values are iteratively adjusted to equalize state populations visited.37–39 Similarly serial replica exchange40 also performs periods of MD simulations in a single replica but the selections of jumping to other thermodynamic states need corresponding estimated potential energy distribution functions at those states accumulated from time series of previous simulations. Other varieties of serial replica exchange such as virtual replica exchange41 distributed replica sampling42 and simulated tempering distributed replica sampling41 also need to estimate similar potential distribution functions in other states during the simulations, which make the massive-scale simulations of complex systems less applicable.</p><p>In recent work we proposed an asynchronous parallel replica exchange (AsyncRE) methodology43 and corresponding python software package44 to utilize massive heterogeneous computing grids without pre-estimation of those free energy weights or potential distribution functions. This AsyncRE framework removes the synchronizing feature of the standard implementations of parallel replica exchange. This asynchronization feature make it possible to optimize the usage of heterogeneous HPC clusters and dynamically distributed computing networks. The main idea of AsyncRE is to use a client-master communication architecture: the local master server assigns all replicas to one of two states, "running" or "swapping" in a dynamic way; "running" replicas are submitted for execution of a predetermined amount of simulation time (an MD cycle or period) on a remote client through ssh or BOINC transportation mechanisms and thermodynamic parameters of "swapping" replicas are periodically exchanged on the local server independently of the "running" replicas; a "running" replica becomes a "swapping" replica once its remote MD cycle is completed and simulation results are returned to the server. This client-master communication mechanism does not require a direct network link between compute nodes due to the fact that all exchanges are managed and conducted by the master process, and the asynchronous exchange feature does not require a static homogeneous set of compute nodes. This AsyncRE framework allows the loss of compute nodes or replica simulation failures and is able to scale to very large numbers of replicas and take advantage of heterogeneous and dynamically distributed computational resources, including XSEDE high performance clusters, university grid networks, and world-wide distributed networks contributed by volunteer computing units such as IBM WCG. Our previously developed Python package44 (https://github.com/ComputationalBiophysicsCollaborative/AsyncRE) combines the AsyncRE framework for job preparations and exchanges with the BOINC server for job assignments and collections from clients using the Python interface/wrap to separate the interactions between the MD engine and the BOINC client/server. In this manuscript, we report the C++ implementation which allows more direct communications between the IMPACT MD engine, the BOINC Client/Server and the AsyncRE framework for massive scale simulations involving hundreds of thousands of replicas.</p><p>The binding energy distribution analysis method (BEDAM)8 developed in our group computes absolute binding free energies for receptor-ligand systems in implicit solvent using a single alchemical decoupling path. It is among several implementations1–7 of a statistical mechanics theory of molecular association equilibria based on atomistic molecular dynamics simulations1,6,45 that are performed to calculate the free energy differences between neighboring thermodynamic states along predefined paths connecting the fully decoupled (apo-protein + free ligand) state to the fully coupled (binding complex) state. Those paths can be realized as in a physical way such as the potential of mean force (PMF) method46 constructing the free energy profile along the distance between the ligand and the receptor, or in an alchemical way such as double decoupling methods (decouple the ligand respectively from the receptor site and the solution environments to vacuum47 from explicit solvent) and single decoupling as in the BEDAM method (decoupling the ligand directly from the receptor site to the solution environment in implicit solvent8) using a parameter (λb) to define intermediate thermodynamic states and rescale (or weakening) the binding energy (including the intermolecular interactions and implicit solvation energies) between the ligand and the receptor. In the early implementation of BEDAM method in the IMPACT48 molecular simulation package, parallel Hamiltonian replica exchange molecular dynamics (HREMD) sampling (λ hopping)26,49 in the synchronous nearest-neighbor replica exchange scheme (SyncRE) was employed to realize the conformational diffusion along the alchemical path connecting the bound (λb = 1) to unbound states (λb = 0) and accelerate the sampling of the intermolecular (external) degrees of freedom between the ligand and the receptor by weakening the binding energies of intermediate states. Recently we have implemented the BEDAM method in the AsyncRE framework43,44 to enable asynchronous replica exchange simulations using heterogeneous HPC clusters and campus distributed computing grids. In this report, we focus on the combination of BEDAM and AsyncRE framework for massive-scale simulations of protein-ligand complexes on the IBM WCG.</p><p>Due to the rescaling (weakening) feature of binding energy using the parameter λb the BEDAM method is capable of carrying out extensive intermolecular conformational sampling of the relative position and orientation of the ligand with respect to the receptor, an advantage of BEDAM8,50–53 over existing free energy perturbation (FEP) and absolute binding free energies protocols in explicit solvent. BEDAM is particularly well suited to refine the prediction results from docking methods and to prioritize the hits for further lead optimization using more expensive free energy simulation methods. Namely docking methods perform the initial search for strong binders using simplified conformational sampling and energy scoring functions in order to perform high throughput virtual screening of large libraries of small molecules (on the order of millions). Then BEDAM simulations can be started from the initial complexes of the top 102 to 103 candidates predicted by docking methods and BEDAM can be used to refine binding poses and provide more accurate estimations of binding free energies which can be used to increase enrichment factors. The top predicted results of BEDAM (on the order of 101 to 102 ) can be passed to carry out explicit solvent FEP simulations with more accurate atomistic simulations of protein-ligand complexes and more detailed physical potentials and solvation environments.1–7
Hence BEDAM targets a niche between low resolution docking methods and high resolution FEP methods in explicit solvent.1,6,45
with the major goal of refining the predicted results from docking methods and significantly reducing the number of candidates to perform further lead optimization using FEP methods in explicit solvent.</p><p>Although post docking analysis using pharmacophore model filters and visual checks can improve initial docking predictions,54 in this report we focus on an alternative BEDAM refinement procedure starting from initial AutoDock predictions using the SAMPL4 library for the HIV-1 integrase (HIV-1 IN).53–55 The SAMPL4 library is a low diversity library of ~300 ligands designed to target the LEDGF allosteric site of HIV-1 integrase, which was the result of a prior lead identification campaign.55 The protocols described here have been implemented for the FightAIDS@Home Phase 2 project, using the physics based BEDAM model and AsyncRE computational method running on the IBM WCG volunteer grid to score the top candidates from the Phase 1 screening using the AutoDock protocol. Massive binding free energy calculations using replica exchange on heterogeneous distributed volunteer computing networks like WCG is challenging in part because of the need for optimization of the parameters which control how the AsyncRE simulation runs on the massive grid to utilize efficiently the resources in a highly dynamic and heterogeneous environment due to the fact that the wall clock time to finish a MD simulation at client side is very diverse. This diversity is not only because of the diversity from computing units (in CPU speeds and device types such as desktop, laptop and mobile) but also based on the fact that the simulation can be paused, stopped and aborted in unpredictable ways. Furthermore, there also exist different kinds of time overheads including the delay of returning results on the client side and others on the server side such as generating and transferring of input files, queuing of individual jobs, and exchanges of waiting replicas. Previous attempts to utilize distributed computing networks to calculate binding free energies were limited to a few examples56–59 using short MD simulations without exchanges to build Markov models utilizing cloud computing. These studies have mostly been focused on protein folding and large conformational changes.41,60–63 To our knowledge no systematic studies for optimizing the simulation parameters and designing general procedures for massivescale simulations of protein-ligand binding have been reported. We fill this gap using as an example the SAMPL4 library for the HIV-1 integrase, one of several important viral enzymes64 during the life cycle of the HIV-1 virus. Integrase is responsible for integrating the viral genome into the host genome with the help of the human LEDGF protein (a transcription factor) linking the HIV-1 integrase to the human chromosome.65 The SAMPL4 library is a focused library developed for lead optimization targeting the LEDGF binding site of HIV-1 IN and preventing the LEDGF protein from binding to stop the insertion of the HIV genome into the human genome, consisting of ~300 ligands sharing several common molecular scaffolds such as benzoic acid, benzodioxole acid, and benzodioxine acid.53,55 Since the SAMPL4 library is not very diverse this presents a particular challenge to binding affinity predictions based on docking alone.</p><!><p>BEDAM is based on the statistical mechanics theory of molecular association6 and the Widom potential distribution theorem.66 The ligand is decoupled from the complex environment with full intermolecular interactions (bound state) to an uncoupled (unbound) state in the implicit solvent environment through a predefined alchemical path which is realized by introducing intermediate states with a parameter (λb) rescaling (or weakening) the binding energy (including the intermolecular interactions and implicit solvation energies) between the ligand and the receptor. The BEDAM method8 calculates the absolute binding free energy ΔGb0 between a receptor P and a ligand L by employing a λb-dependent effective potential energy function as follows:
(1)VÎ»b(r)=V0(r)+Î»bu(r)
where V0(r)=V(rP)+V(rL) is the total potential energy of the complex including both the receptor P and the ligand L at the uncoupled (unbound) state in the implicit solvent, and u(r)=u(rP,rL)=V(r)−V0(r)=V(rP,rL)−V(rP)−V(rL) is defined as the binding energy for each complex conformation (r=(rP,rL)), corresponding to the difference between the total effective potential energies V(r) with implicit solvation effects67–69 of the fully coupled (bound) and decoupled (unbound) states of the complex with the same fixed internal conformations. The standard free energy of binding for this system can be calculated as6,8
(2)ΔGb0=−kBTlnC0Vsite+ΔGb=−kBTlnC0Vsite−kBT∫dup0(u)e−βu,
with the first ideal term on the right −kBTlnC0Vsite=−kBTlnVsiteV0 representing the entropic change when moving the free ligand in the volume ( V0 ) of standard concentration into the volume of the binding site ( Vsite ) and the second one ΔGb=−kBT∫dup0(u)e−βu denotes the total free energy change when turning on the interaction energies between the receptor and the ligand. Moreover kB is the Boltzmann constant, T is temperature, C0 (=1M) is the the standard concentration of ligand molecules. and p0(u) is the probability distribution of binding energy (u(r)) collected in an appropriate decoupled ensemble of conformations in which the ligand is present in the binding site but the interactions between the receptor and the ligand are turned off but both molecules are present in the solvent continuum. To obtain a good estimation of p0(u) for biomolecules with slow dynamics, a series of discretized λb values between 0 and 1 are introduced to represent intermediate thermodynamic states and depict an alchemical path decoupling the ligand from the complex environment with full interactions in the receptor site (λb = 1) to a pure solvent state without protein interactions (λb = 0). In practice the ΔGb can be calculated by reweighting analysis of all snapshot binding energies of all thermodynamic states using the unbinned UWHAM70 or similar methods.71–73</p><p>Comparing with the docking methods with empirical docking score functions, one of the advantages of the BEDAM method is the ability to capture the entropic changes for binding processes8,50,51. In the BEDAM method, the binding process can be decomposed into two separate steps: 1) the ligand and the receptor reorganizes its conformational (internal and external) ensembles in the unbound state to match corresponding distributions in the bound complex; 2) the receptor-ligand interactions are created in the binding site with no entropic change involved due to the fact of no change in the configurational ensemble of the binding partners. The free energy change for the second stage can be calculated as the average binding energy at the fully coupled state (λb = 1), represented by ΔGIIo=<u>RL. We should emphasize that ΔGIIo not only includes the interaction energies between the ligand and the receptor but also includes the solvation contribution (the difference of implicit solvation energy between the bound and unbound state). In contrast, there are both entropic and enthalpic changes in the first stage due to the reorganization of conformational ensembles, denoted as the reorganization free energy ΔGIo=ΔGreorgo. Hence the standard binding free energy can be described as
(3)ΔGb0=ΔGIo+ΔGIIo=ΔGreorgo+<u>RL.</p><p>We should mention that using of implicit solvation model is crucial to achieve the decomposition of binding free energy and the single-path decoupling of the BEDAM method. Namely the ligand is only decoupled from the coupled complex environment to the free state in the solution and the additional decoupling from the solution to the vacuum as in the double decoupling method is not necessary since the solvation energies can be evaluated from the ligand and receptor conformations through the AGBNP2 solvation model.67–69 Implicit solvent models could also lead to faster convergence of free energies than explicit models due to lack of solvent fluctuations although simulations with implicit solvation are more difficult to parallelize due to the relatively small number of atomic sites and the complexity of algorithms to evaluate the Born radii.67–69 For some systems which require explicit representation of water mediated interactions around the binding site, simulations using implicit models may not fully capture the ligand binding effects and become less accurate than that using explicit solvation models as in the double decouple method.47</p><!><p>The early implementation of the BEDAM method was accomplished in the SyncRE scheme through the MPI implementation in the IMPACT48 simulation package targeting homogeneous HPC clusters where all CPUs involved in an RE simulation have the same speed but with a limited number (tens) of replicas due to the necessity of synchronization of replica exchange. To overcome this synchronization limitation and enable large-scale simulations (thousands) using heterogeneous HPC clusters and distributed computing resources, recently we proposed the AsyncRE framework43and developed a corresponding python package44 where only the replicas in the local waiting list can participate in the exchange process using an algorithm for the sampling of the state permutation space which does not require the prior identification of neighboring states and attempts to exchange two replicas are randomly picked but follow the same Metropolis criterion and repeated many (M3 to M5 where M is the total number of replicas on the exchange queue) times to reach the infinite swap limit.43,74,75 This Python package (https://github.com/ComputationalBiophysicsCollaborative/AsyncRE) has been installed to perform AsyncRE simulations on XSEDE high performance resources, and on BOINC distributed computing networks at Temple University, Brooklyn College at CUNY. 43,44,76 The Python interface/wrap separate the direct communication between the MD engine and the BOINC server and there is also no direct data transferring between the MD engine and the BOINC client for the screensaver visualization. In this way we can minimize the tasks for code development of new implementations using different MD engines, which is more suitable for setting up small or median size computing grids up to ~ 103 computing units. For massive computational grids such as the IBM WCG grid involving hundreds of thousands of computing units, in this report we introduce the C++ implementation (See Fig. 1 for the corresponding workflow for Figh-tAIDS@Home Phase 2 refinement) which allows more directly communications between the IMPACT MD engine and the BOINC Client, and between the BOINC Server and the AsyncRE framework for massive scale simulations involving hundreds of thousands of replicas, although it shares the same framework of AsyncRE algorithm with the Python implementation for HPC clusters and university grids. We have summarized the major differences and provided more implementation details in Supporting Information.</p><!><p>There is no universally accepted metric for evaluating the sampling efficiency in REMD simulations but there are several methods described in previous work such as the replica relaxation correlation time and end-to-end transit time of the replica state index, 74 variances of the estimated means of relevant observables, 78–80 slowest eigenvalue of the Markov chain from analyzing the state transition matrix,81 root-mean-square deviation (RMSD) of related observables of test simulations from the corresponding reference simulation,82,83 and our previous statistical inefficiency analysis43 by extracting the total effective relaxation time of the binding energy at λ = 1.0 (λb in Eq.1 is shortened as λ for better notation here after), which includes all relaxation effects both from MD simulations and replica exchange mixing. To include all of the relaxation times at different λ values and the replica exchange mixing time, we borrow the idea of the diffusion coefficient in coordinate space84 and define the diffusion coefficient of the replica exchange system in the λ space as
(4)Dλ=<(λi−λj)2>2<τii>,</p><p>where the average is (<···>) calculated through all λ states, τii is life time of state λi measuring the total simulation time of a replica stays in the λi state before jumping to any different state λj,|λi−λj| is the jump difference of a replica from the λi state to any different state λj.</p><!><p>The toy model complex is b-cyclodextrin-heptanoate as depicted in Fig. 2 was used to benchmark the ported IMPACT MD engine and the WCG implementation of AsyncRE framework. It is one of host-guest systems well investigated in our previous algorithm developments and benchmarking work.51,52 Besides the standard tests for the IMPACT MD engine to reproduce single point energies, we performed standard 1D BEDAM simulations at 16λ values (0.0, 0.001, 0.002, 0.004, 0.01, 0.04, 0.07, 0.1, 0.2, 0.4, 0.6, 0.7, 0.8, 0.9, 0.95, and 1.0) using the OPLS-AA force field,85,86 the AGBNP2 implicit solvent model,69 and three different implementations of RE framework: a) SyncRE using MPI targeting homogeneous HPC clusters; b) AsyncRE framework implemented by Python interface for median-size campus computing grids such as Temple Grid; c) AsyncRE framework implemented by C++ for massive-scale distributed grids such as the IBM WCG grid. The MD period for the SyncRE simulations is 100 steps (0.2 ps) and 10,000 steps (20ps) for both types of AsyncRE simulations. From Fig.2b we can see that all simulations (after 10ns) converge to the consistent binding energy distributions at λ =1.0 which supports the correctness of the different implementation of RE frameworks (SyncRE on NSF XSEDE and other HPC clusters, AsyncRE on heterogeneous campus grids. and massive AsyncRE on IBM WCG).</p><!><p>The 301 ligands (53 experimental binders and 248 nonbinders, 451 protein-ligand complexes considering the protonation states and other tautormeric extensions of ligands) of HIV-1 IN complexes were prepared previously53 for the HIV-1 integrase virtual screening SAMPL4 challenge54,55 with the aim of predicting likely binders to the LEDGF site of integrase87 using the BEDAM method for binding free energy calculations8,50–52.</p><p>Docked structures for all ligands were obtained by the flexible docking procedure using AutoDock Vina for the SAMPL4 challenge54 and initial reference crystal complex with the PDB ID 3NF8 as shown in Fig. 3. Only the best docking conformation of each complex was selected to perform the binding free energy refinement for the original simulations using the SyncRE framework on HPC resources and submitted to the SAMPL4 challenge.53 The tremendous computing resources of IBM WCG grid allow us to include the remaining top 8 predicted poses from AutoDock Vina in addition to the best one for each of 451 protein-ligand complexes and extend the number of replicas per complex from the original 20 to 100 as listed in the Supporting Information. This corresponds to 900 replicas per complex (100 replicas/pose × 9 poses). Furthermore all simulation parameters, force field parameters, and restraint methods were kept the same as the previous work53 as described in more detail in the Supporting Information. We performed different types of AsyncRE BEDAM simulations for the SAMPL4 ligands as listed in Table 1 using the WCG resources. The first set of WCG simulations for the SAMPL4 library were carried out using the independence sampling scheme with the same 20 λ values as the original SAMPL4 submission but without exchange feature and the simulations were extended to all 9 predicted poses from AutoDock. A large MD period of 100,000 steps was selected in order to minimize the number of resubmissions (hence the overhead) due to the fact that an individual job on the client side is limited to a few hours of the CPU running time (after removing any pauses and interrupts) and a maximum return time of 7 days (namely a job is considered as failed if it can not be returned in 7 days) based on the previous setting for other projects on WCG. The remaining two sets of WCG simulations of the whole SAMPL4 library use the AsyncRE framework to perform exchanges among 100 replicas ( 100 λ states) but with different MD periods of 10,000 (20 ps) and 50,000 (100 ps) steps per exchange cycle respectively. For a comparison, we also include flattening AsyncRE simulations that were carried out previously76 on the Temple University Grid.</p><!><p>In the conventional parallel synchronous replica exchange (SyncRE) framework, each replica is assigned a different thermodynamic state (a λ value in this report) and periodically, after all of the replicas complete a given number of steps (MD period), replicas attempt to exchange their current state assignments with typically one of their nearest neighbors at adjacent thermodynamic states. Exchange attempts are accepted based on well-established microscopic reversibility requirements.18–22,32,33 To improve the efficiency and convergence, conventional REMD has been extended mostly through the modification of the Hamiltonian 24,25,27–29,36,88–95 as for example our recent addition of biasing potentials for BEDAM simulations.76 On the other hand, some studies focused on expanding the exchange dimensions19,21,36,96 and optimizing the setting of simulation parameters not limited to the number and distribution of replicas (such as A values97 and temperatures18,98–100), the number of exchange attempts,74,101 and exchange frequency (MD period).82,83,102,103 Although the efficiency and convergence analysis of SyncRE MD in comparison with conventional MD without exchange has been performed in many previous studies.78–80,104–108 there still exist debates on how to best select simulation parameters such as the length of individual MD simulations (MD or exchange period) and the number of exchanges attempted after a MD period. Early results showed that the efficiency of REMD could be significantly reduced when the MD period is smaller than a certain size,104–106 while recent studies82,83 found that the efficiency increases monoton-ically as the MD period becomes smaller. Previous findings74,101,106,109 also illustrated that the number of exchange attempts should be chosen as frequently as feasible so as to reach a so-called infinite swap limit.75,102,103,110</p><p>The characteristics of the WCG distributed computing grid are unsuitable for SyncRE, we implemented an asynchronous RE protocol which differs from SyncRE in various respects43,75: 1) Only a portion of all simulated replicas are available locally on the server side to participate in the exchange process and the rest (half of the total number of replicas for our WCG simulations) reside on the client side for individual MD simulations; 2) Because the randomness of available replicas makes the nearest-neighboring exchange scheme of the SyncRE not feasible, we developed a random pairwise exchange scheme whereby many pairs of replicas (not limited to the nearest neighboring replicas) in the waiting pool are randomly selected for exchange attempts so as to reach the infinite swap limit;43,74,75 3) AsyncRE simulations generally require relatively large MD periods (> 1 picoseconds = 500 MD steps) to minimize overheads (such as preparation and transportation of input and output files between the coordination server and the computing clients, random pauses of the computing client, and the job queuing time on the server) relatively to the running time of a MD period. Due to the characteristics of AsyncRE and the dynamic and heterogeneous nature of computing grids, it is difficult to perform a direct comparison on efficiency and convergence between the AsyncRE simulations on the WCG and the SyncRE simulations. To optimize and guide our AsyncRE simulations on WCG, we consider two simulation parameters for our benchmarking tests using the AVX38789 ligand from the SAMPL4 library: a) the total number of replicas varying from 20 to 1000 (Note that Table 1 only lists the simulations using 20 and 100 replicas for the whole SAMPL4 library.) to test the effect of larger numbers of replicas which can result in more overlaps among the binding energy distributions of replicas and increase the accepted ratio of exchanges; and b) the MD period varying from 10 thousand to 50 thousand steps to reduce the communication overhead at the expense of frequency of exchanges but this might lead to slower equilibration.</p><p>Table 2 shows the accepted ratios of exchanges (the total number of successful exchanges divided by the total number of attempted exchanges) calculated from different benchmarking tests of WCG AsyncRE simulations for the AVX38789 ligand in complex with HIV-1 IN. While its value is not affected by the MD period, the accepted ratio almost doubles as the number of λ states is increased from 20 to 1000 albeit not in a linear way. For example, the accepted ratio of accepted changes to exchange attempts increases only slightly (0.2 to 0.22) when the number of λ states is increased from 100 to 1000. The increase of accepted ratio is based on the fact that increasing the number of λ states increases the overlap of binding energy distributions of two exchanged replicas. The corresponding accepted ratios calculated from the SyncRE simulations for the original SAMPL4 submission are around 0.5 where the exchanges are attempted only between nearest neighboring states (with more overlaps in binding energy distributions) and results in more successful exchanges when the total number of attempted exchanges is fixed. Instead, AsyncRE simulations require many more attempted exchanges to reach a similar swap limit.43,74,75 Due to the insensitivity to the MD period, the accepted ratio is not a good quantity to optimize the simulation parameters. Similar trends can be found on the mean square jump of λ values (<|Δλij|2>) as shown in Fig. 4a. On the contrary, the mean life time of replicas (<τii> in the unit of picoseconds) decreases as the MD period is decreased and its value approaches to a limit (the MD period) as the number of replicas is increased to 1000 as exhibited in Fig.4b. Figure 4c displays the diffusion coefficients (Dλ) in the A space (defined in Eq. 4) which include both effects of the mean square jump of λ values and the mean life time of replicas. It is pronounced that the diffusion coefficient decreases as the MD period is increased when the number of replica is fixed. The diffusion coefficient also exhibits a large increase (almost doubled) when the number of replicas is increased from 20 to 100. However, there are no significant changes (only ~10%) to the diffusion coefficient as the total number of replicas is increased from 100 to 1000 with the MD period is fixed. Based on this observation we selected 100 replicas for each AsyncRE WCG simulation of the whole SAMPL4 library. Because there are more than 4 fold increase of the diffusion coefficient, in principle we should select 10 thousand steps for the MD period instead of 50 thousand if there exists no overhead (or the same percentage of overhead) for exchanges and queuing on the server between two adjacent MD periods. However, reducing the MD period will increase the number of cycles required of replica resubmission (when the total number of simulated MD steps is fixed for each replica) and increase the percentage of the overhead that dynamically depends on the available WCG resources and the settings of clients and server. We will discuss the effects of the MD period in more detail in the following section where we analyze two sets of WCG AsyncRE simulations for the whole SAMPL4 library using MD periods of 10 and 50 thousand steps respectively.</p><!><p>We conducted two sets of AsyncRE simulations on WCG (with MD periods of 10K and 50K respectively). Each set includes 451 complexes and 9 predicted poses from AutoDock for each complex. There are 100 replicas for each BEDAM simulation started from each pose. In total, there are 2 × 451 × 9 × 100 = 811, 800 replicas and with half of them (405,900) submitted for MD simulations on the WCG client side and half of them checkpointed to disk and participating in exchanges on the WCG master server side. The total available computing units for the FightAIDS@Home Phase 2 project fluctuates around 100,000. Therefore at any one time roughly three quarters of submitted MD jobs are waiting for execution on the WCG BOINC server and only one quarter are running on client CPUs. To perform the overhead analysis we define four quantities that can be found or calculated from the returned info files of WCG simulations: 1) the "MD CPU time" which represents the wall clock time for a computing unit to complete a period of MD simulation of 10K or 50K steps without any pauses (temporally idling due to the occupation of the volunteer CPU by other tasks) on the client side; 2) the "replica client time" which denotes the total wall clock time that a replica resides on a client node including restarting if any from a crashed job and the delay in returning a completed MD job (A MD job is considered as failed if it can not be returned within one week); 3) the "replica server time" which accounts for the total wall clock time that a replica stays on the WCG server side including the time spent on the exchange process of AsyncRE and the time spent on the MD execution queue of BOINC server when the replica is selected for the next cycle of MD period after exchanges; 4) The "replica cycle time" (=replica client time + replica server time) which is calculated from the wall clock time difference between two sequentially returned MD jobs for the same replica (Its mean value multiplied by the total number of RE cycles is equal to the mean value of wall clock time to finish the total simulation per replica). There is no direct way to obtain the replica server time (from the returned info files of WCG simulations) but we can get the value by subtracting the replica client time from the corresponding replica cycle time after performing a one to one match.</p><p>The probability distributions of MD CPU times calculated from two sets of WCG AsyncRE simulations are displayed in Figure 5a. The most probable value (mode) for a client to complete a MD period of 50K steps is ~ 5 hours, almost 5 times longer than that for 10K steps (~1 hour) which is also consistent with the mean values as listed in Table 3. Figure 5b shows the corresponding distributions of replica client times, with the most probable values ( ~ 1.1 and 9.1 hours for 10K and 50K respectively) not more than twice larger in comparison with the corresponding values of MD CPU time as listed in Table 3. The mean values of the replica client times behave in a different way for these two sets of WCG AsyncRE simulations: The mean value from 10K simulations ( ~ 6.9 hours) is almost 5 times larger than the mean value ( ~1.3 hours) of MD CPU time but there is only twice larger for 50K simulations. Namely the 10K simulations have much larger overhead (5 .vs. 2) on the client side due to the fact that although the returning time of a MD job can be distributed in a wide distribution due to the heterogeneity of CPU speeds and operation settings of client computers, the maximum return time for a MD job is the same (7 days) for both sets of simulations as exhibited by the unsymmetrical long tail distributions in Figs.5b. Similar distributions of replica cycle times are shown in Fig. 5c. The mode values (see Table 3) from the simulations of 10K and 50K MD periods are around 60 and 70 hours respectively and there is an exponential tail at the longer time side due to the random selection of exchanged replicas to be submitted to a MD simulation on the client side, with the mean values of ~74.6 and 113.2 hours (for 10K and 50K respectively). The reason that those two modes are so close is that the replica server time (calculated by subtracting the replica client time from the replica cycle time) is more than 60 hours for those two sets of simulations, and their mean values (~67.7 and 99.6 hours for 10K and 50K respectively) take the major part (91% and 88% for 10K and 50K respectively) of the mean replica cycle time (~74.6 and 113.2 hours for 10K and 50K respectively) as listed in Table 3. For an ideal case when there is no overhead from the server, theoretically the minimum value of mean replica sever time should be the same as the mean value of replica client time since only half number of total simulated replicas are assigned to run MD jobs on the client side and half of them stay on the server side. Namely, on average a replica has to wait for an additional period of replica client time to be resubmitted for the next cycle of MD simulation. Hence, the minimum value of the mean replica cycle time is twice ( ~ 13.9 and 27.2 hours for 10K and 50K respectively) as the mean value of the replica client time. Assuming that we can reduce the mean value of replica server time by a factor of 4 if we do not over submit the WCG grid and the mean value of replica client time is the same as before, then the mean values of replica cycle times are ~23.8 and 38.5 hours for 10K and 50K respectively (~1.7 and 1.4 times of the minimum values of mean replica cycle time). Hence, when considering the overhead from both the client and the server sides together, the final overall overheads are ~8.5 (5 on the client side × 1.7 on the server side) and 2.8 (2 on the client side × 1.4 on the server side) times for 10K and 50K simulations respectively when the WCG grid is not over submitted. Because of the heterogeneity of WCG simulations and the different percentages of overhead as illustrated in Fig.5 and Table 3, the average value of the total wall clock time to finish a total of 600,000 MD steps (1.2ns) for each replica using 10K MD period is around 3 times larger than that when using a 50K MD period.</p><p>Although the percentage of overhead is considerable large (8.5 and 2.8 times for 10K and 50K respectively without over-submission) , the advantage of WCG simulations becomes more obvious when we consider the total number of available computing units and the total number of simulations we can submit. Table 4 lists the estimations of wall clock time for different types of BEDAM simulations performed for all 451 SAMPL4 HIV-1 IN complexes. The total MD length per replica has been normalized to 1.2ns for a consistent comparison. Our original SAMPL4 submission required ~2 millions service unit (SU, 1 core per hour) in total on XSEDE HPC clusters and 48 days to finish all jobs when using 1000 cores at the same time. One set of WCG AsyncRE simulations for the SAMPL4 library needs 90M SUs (around 9 years allocations for a big computing project) and would take ~2160 days to finish on XSEDE if the same 1000 cores are used to running all jobs without considering the queuing time. Using the current available resources (~100,000 cores) for the FightAIDS@Home Phase 2 project the required time to finish 1.2ns per replica for one set of AsyncRE simulations is only ~56 days for the 50K MD period and ~180 days for 10K MD period respectively. We should emphasize that the number of required days can be reduced greatly ( ~19 and 59 days respectively) if the WCG grid is not over submitted (with ~400,000 CPUs available) to reduce the average queuing time on the BOINC server (and the replica server time), by ~14 and 36 days respectively for the ideal case without any overhead from the server. Furthermore, a significant reduction can also be achieved on the client side by using a much smaller value of the maximum return time of MD jobs and asking the clients return the completed jobs as soon as possible.</p><!><p>To perform quantitative comparisons for the retrospective predictions of the whole library from different types of BEDAM simulations, we calculated the receiver operating characteristic (ROC) curves as shown in Fig. 6 using the experimental crystallographic data available.55 Note that there are 53 binders based on the available crystal structure data for the SAMPL4 library and they are expected to be weak binders since most of them have affinities more than 200 μM via surface plasmon resonance (SPR) and only 8 ligands were provided binding affinities by the SAMPL4 Challenge organizers.55 To aggregate the WCG simulations which include 451 complexes (301 unique ligands plus protonation and tautomeric states) and 9 poses for each complex, we have applied the following strategy: a) select the minimum value of binding free energies calculated from 9 BEDAM simulations started from 9 different poses for each complex; b) then select the minimum value from the simulations with different chemical extensions (such as the different protonation and tautormeric states) but for the same ligand as done for previous work.53,76 The ROC curves calculated from the binding free energy (ΔGbo) predictions are shown in Fig. 6a and corresponding values of the area under the curves (AUC) are listed in Table 5. It can be seen that the free energy scores from all BEDAM simulations (obtained with or without exchange) increase the AUC value in comparison with that from the AutoDock predictions. The AUC value from the BEDAM simulations with Independence sampling (represented by "WCG-Indp") have the smallest increase (~9%). The AUC values calculated from the other two sets of WCG simulations with the exchange feature have larger improvements of ~20% (50,000 MD steps, denoted as "WCG-50K") and ~22% (10,000 MD steps, denoted as "WCG-10K"). The most significant improvement (~50%) is from the BEDAM simulations with both the selected intramolecular torsional and nonbonded interactions flattened (represented by "Flat-All"). We should emphasize that it is well-known that docking methods are better suitable for libraries with very diverse ligands instead of the focused libraries as is the SAMPL4 library in this report. Furthermore, the AutoDock predictions are the raw output and there are post processing procedures such as pharmacophore model filtering and visual checking that can improve the docking predictions.5453,54 Our BEDAM binding free energy prediction provides an alternative approach to the refinement of docking results. In real applications, our BEDAM calculations should start from the best docking predictions after post analysis of filtering and prioritization using other fast models. In comparison with the AUC value from the original SAMPL4 BEDAM simulations using SyncRE framework on XSEDE,53 both WCG AsyncRE simulations with exchange show only slight improvements (~1% for WCG-50K and 2% for WCG-10K respectively) but have already achieved the same rank (2nd out of 25 submissions) in SAMPL4 Challenge. This result is mainly due to two facts: 1) The WCG AsyncRE simulations have an exchange frequency more than 100 times slower than the SyncRE simulations on XSEDE; 2) Although the AsyncRE simulations on the WCG have 45 times larger aggregated simulation time (9 poses and 5 times more lambda values), simulations from pose 6 to 9 (as we show below) have worse predictions than that of pose 1 to 5 and reduce the effective aggregated simulation time to only 25 times larger. The corresponding ROC curves using the average binding energy (ΔGIIo=<u>RL) predictions are illustrated in Fig. 6b. It is obvious that all AUC values using average binding energies from all BEDAM simulations are greatly reduced in comparison with that from the binding free energy predictions. Especially at the cutoff of top 20% of inactives the distinguishing binders from nonbinders is close to a random selection. The inferior prediction from the average binding energies implies that besides the binding energy component the reorganization free energy component is also crucial and they are correlated in a way to enhance the discrimination between binders and nonbinders.</p><p>The AUC values focus on the performance of predictions for the whole library. Instead the enrichment factor measure as defined in the Supporting Information is better suitable for evaluating the prediction power of early recognition from a library, which tells us how many more actives (binders in our analysis) can be found within an "early recognition" fraction (such as 20%) of the ordered library based on a given score relative to a random distribution. Table 5 also shows calculated enrichment factors at 20% (EF20) and 40% (EF40) cutoff. Similar observations as the AUC values can be found: All EFs using binding free energies (ΔGbo) calculated from BEDAM simulations have more than 20% improvements in comparison with the raw AutoDock predictions; The WCG simulations with independence sampling have the smallest increases (21% for EF20 and 26% for EF40); The EF values for the two sets of WCG AsyncRE simulations show improvements of 20~30% for EF20 and 40~50% for EF40; Obviously the flattening simulations have the largest increases (46% for EF20 and 53% for EF40). From Table 5 it is also clear that the EF values using average binding energies are worse than that from binding free energies indicating the importance of reorganization free energies.</p><!><p>As with most free energy calculation methods, the BEDAM method relies on docked poses or crystal structures to start a simulation. The tremendous WCG resources allow us to refine the docking predictions not only based on the best predicted pose but also alternative poses according to the docking scores. The WCG results shown in Fig. 6 and Table 5 were obtained by selecting the minimum (most favorable) binding free energy scores among the 9 scores calculated from individual BEDAM simulations started from each of 9 different poses of the same complex. In Fig. S2 of Supporting Information we display corresponding ROC curves of the binding free energies calculated from all trajectories of 9 poses by simple averaging (denoted as pos avg). It is clear that such equal treatment for all 9 poses results in a worse prediction of AUC values as listed in Table 5, implying that many simulations are still trapped in a local equilibrium state instead of reaching a global one. Such local equilibrium properties become more obvious when we check the Spearman ranking order correlations between the results calculated from 1.2ns trajectories and that from 3.0ns as shown in Fig. S3 of Supporting Information. Both the WCG AsyncRE simulations using a MD period of 50K steps and the flattening BEDAM simulations have rank order correlations larger than 0.95 (highly correlated) at two different simulation times. To optimize the computing resources for future simulations we are interested in the prediction performance from individual docking poses in addition to that aggregated from all 9 poses for the same complex. Fig. S4 exhibits the ROC curves of binding free energies calculated from the WCG AsyncRE BEDAM simulations using the MD period of 10K steps and started from each one of 9 different poses predicted by AutoDock (See Fig. S5 in Support Information for corresponding simulations using 50K MD period). The corresponding AUC values are plotted in Fig. 7 (see Table S1 for the values in detail). Generally the AUC values calculated from individual poses are smaller than that obtained by aggregating to the minimum value of binding free energies from all 9 poses of the same complex implying that aggregation by the minimum value does help to improve the accuracy of overall prediction. They are also correlated among the different WCG simulations (See Fig. S3 for the comparisons of individual binding free energies) which is consistent with the fact that different types of simulations were started from the same individual poses. Using a smaller MD period increases the AUC values, which is also in agreement with the diffusion coefficient analysis shown in 4c. A more careful comparison also shows that the AUC values have large fluctuations from pose 1 to 5 implying the best pose (top 1) predicted by AutoDock might not be the best candidate for BEDAM simulations (for example, WCG BEDAM simulations from the top 2 and top 5 poses have larger AUC values). We can also find that the AUC values decrease to lower values from pose 6 to 9 which means that it is not necessary to include all 9 poses into the WCG AsyncRE simulations and top 1 to 5 poses from the AutoDock predictions are good enough for this SAMPL4 library. We should point out that this is only an observation based on Fig. 7 and for a more rigorous comparison we should perform the accumulative analysis of AUC values by aggregating the minimum values of binding free energies from top 1 to top 9 poses. The conclusion might differ depending on which ligand and also the parameter selection for the AutoDock predictions. In fact for more advanced BEDAM simulations using flattening potentials, we can only use one pose but achieve a high AUC value as shown in Fig.6 and Table 5.</p><!><p>We implemented and optimized the AsyncRE framework of BEDAM method on the IBM WCG volunteer grid for massive-scale binding free energy calculations using the OPLS-AA force field for describing receptor-ligand complexes and the AGBNP model for implicit solvent effects. We discussed the general procedure on how to select the simulation parameters of BEDAM and refine the AutoDock prediction. We also performed the first massive-scale binding free energy calculations using distributed computing gird and AsyncRE framework with 2 × 451 × 9 = 8,118 complexes and 811, 800 replicas submitted to the WCG volunteer gird at same time. Namely for the WCG AsyncRE simulations we extended the number of lambda states to 100 and included 9 poses from the AutoDock predictions for each HIV-1 IN complex to remediate the slower diffusion in lambda space due to the use of much longer MD periods than the HPC SyncRE simulations to reduce the overhead. The free energy scores obtained from WCG AsyncRE simulations are comparable with those from the SyncRE simulations and show significant improvements over the initial AutoDock predictions although they are worse than those from more advanced REMD that includes the flattening potentials to selected degrees of freedom of ligands. Running asynchronous REMD on WCG is a tradeoff: We sacrifice rapid simultaneous RE exchanges (in depth) with little overhead in order to gain many more simulations (in breadth) on heterogeneous hardware but more overhead per replica.</p><p>The binding free energies from the flattening AsyncRE BEDAM simulations show significant improvement on discriminating binders and nonbinders in comparison with the initial AutoDock predictions and BEDAM simulations without the flattening feature. Using flattening simulations can also reduce the number of poses for each simulated complex and the requirements on computing resources since the flattening potentials can help to accelerate sampling of internal degrees of freedom of the ligand/receptor and adjust the pose during a simulation. However there is also a tradeoff between the selected number of degrees of freedom for flattening and the simulation time to sample all of the relevant space similar to the REST methods in FEP calculations.27,88 Namely improvement from the convergence of simulations in the relevant coordinate space might be weakened by including unnecessary internal coordinates into flattening. 27,88 Hence, using flattening BEDAM simulation is the best strategy to refine the AutoDock results when we have prior knowledge about ligand libraries and/or particular receptors. For example we can perform cluster analysis of the binding pocket to identify related degrees of freedom of receptor flexibility using crystallographic data, MD simulations, or induced-fit docking studies. On the other hand we can consider the degrees of freedom from ligand flexibility using our chemical knowledge and other studies such as flexible docking that can specify the rotatable dihedral angles to explore the conformational flexibility of ligands. The flattening BEDAM simulations reported here were performed in a AsyncRE framework using a local campus grid on Temple University including 5000 CPUs. The implementation of flattening BEDAM on the IBM WCG grid is underway.</p><p>The use of AsyncRE simulations without flattening on the WCG Grid is a good choice when we do not have sufficient prior knowledge concerning which degrees of freedom to apply flattening potentials to or when the number of complexes is very large and an automatic workflow to select related internal degrees of freedom is not available. Considering the smaller overhead and insignificant reduction to AUC values, choosing an MD period of 50K and 100 λ states to run on the WCG is a better choice than that of 10K. We also found that the AUC values from the BEDAM simulations started from individual poses decrease significantly at large number of pose (after 5). This means that we could perform simulations for the 5 top poses predicted by AutoDock which will reduce the requirement of computing resources close to 45%. We note that a factor of ~4 × attributable to overhead is from the waiting (queuing) times on the IBM BOINC server after jobs are submitted (this is included in the replica server time). The reason for this is that the number of jobs submitted to the grid is more than three times the number of CPUs available (~0.1 million) to our FightAIDS@Home Phase 2 project (though there are ~3.0 million CPUs in the entire WCG grid). This part of the overhead can be removed when the available CPUs are increased to a desired number or the number of jobs is reduced using fewer poses or ligands. The real overhead (~8.5 and 2.8 times for 10K and 50K respectively without over-submission) can also be reduced further by using a much smaller value of the maximum return time of MD jobs on the client side or/and increasing the frequency of checking returned MD jobs and performing exchanges on the server side. In fact, we were able to tweak the settings of client machines and the server in such a way that the overall is reduced to less than a factor of 2 × on our local Temple University Grid. We should point out that although the largest benchmark test on the SAMPL4 library is 1000 lambda values per complex, the AsyncRE framework and computing resources support running hundreds of thousands of replicas per complex if it is required such as in multi-dimensional REMD simulations. For the WCG simulations of the whole SAMPL4 library, we submitted jobs including 8,118 complexes and 811,800 replicas at the same time. In principal, we can use all 811,800 replicas for the same complex. We did not run the simulation in this way because the gain from increasing the total number of replicas becomes less significant after its value is increased to 100 based on the diffusion coefficient analysis.</p><p>In summary, combining the BEDAM method with the distributed computing network enables the massive-scale refinement directly on docking predictions and generates good candidates for further lead optimization using more high resolution FEP methods in explicit solvent. We expect that such protocols will become more routine in the near future as massive-scale distributed computing resources such as WCG gird, Amazon Cloud, and Internet of Things, become more widely available.</p>