%%%%%%%%% PROPOSAL -- 15 pages (including Prior NSF Support)
% use letters for the subsections with this command
\def\thesection{\Roman{section}.}
\def\thesubsection{\alph{subsection}.} 
\def\thesubsubsection{\alph{subsection}.\roman{subsubsection}}

%\title{\large \bf GPU-accelerated calculation of ligand binding affinities to biological macromolecules}
%\date{}
%\maketitle
{\large \bf GPU-accelerated calculation of ligand binding affinities to biological macromolecules}
\vfill
\noindent {\small
\begin{tabular}{lll}
{\bf John D. Chodera} & Distinguished Postdoctoral Fellow, QB3, University of California, Berkeley & PI\\
{\bf Michael R. Shirts} & Assistant Professor, Department of Chemical Engineering, University of Virginia & co-PI\\
{\bf David L. Mobley} & Assistant Professor, Department of Chemistry, University of New Orleans & co-PI\\
{\bf Kim M. Branson} & Research Scientist II E, R\&D, Vertex Pharmaceuticals & collaborator\\
{\bf Vijay S. Pande} & Associate Professor of Chemistry, Stanford University, Simbios & collaborator
\end{tabular}
}
\vfill

{\bf A detailed physical understanding of how biological macromolecules bind their ligands (or inhibitors) with affinity and specificity would be of tremendous value to science and medicine.}
Predictive models would facilitate the rational design of pharmaceuticals and tools for chemical biology, as well as provide insight into the biophysics of binding, catalysis, and allosteric regulation.  
Alchemical free energy calculations using molecular mechanics forcefields have made great strides in this direction, allowing binding affinities to be computed using a rigorous but efficient treatment of the relevant statistical mechanics using molecular dynamics or Monte Carlo simulations.
However, up to now, the tremendous expense of these calculations, requiring cluster-weeks of computer simulation for even simple systems, has prohibited large-scale tests of the \emph{accuracy} of current-generation force fields for this task.

{\bf The rise of general-purpose GPU computing is changing the landscape of high-performance computing} by providing immense floating point computing power: a single NVIDIA Tesla C1060 GPU has a theoretical peak of over 1 TFLOP, which the next-generation Fermi architecture promises to double.
Exploiting this level of parallelism remains challenging; the application must be efficiently mapped onto this architecture to approach peak utilization.
Simbios, the NIH Center for Physical Simulation at Stanford University, has succeeded in producing an open library (OpenMM) for accelerating molecular dynamics calculations on GPUs, obtaining over 200x speedup when compared to the single-CPU core performance of standard optimized simulation packages like AMBER using a single C1060 GPU. 
Through close collaboration with the OpenMM developers, we have built a new, freely-available open-source code (YANK) making use of the OpenMM toolkit to accelerate alchemical binding free energies, providing an open extensible platform for experimenting with enhanced sampling algorithms and aiding in the rational physically-driven design of new ligands for biological macromolecules.
{\bf The use of multiple GPUs allows a reduction in the wall clock time required for a single rigorous absolute alchemical binding free energy calculation to a few hours.}  
Preliminary calculations run using a startup allocation on NCSA Lincoln also confirm the feasibility of free energy calculations with generalized Born-surface area (GB/SA) in under a day, as well as the observation that OpenMM provides 200--250x speedup over single CPU cores.

{\bf We propose a number of experiments that will facilitate the development of methods for the quantitative computation of small-molecule binding affinities and the rational, physically-driven design of potent, selective ligands using alchemical free energy calculations in combination with GPU-accelerated computing.}  
These experiments will allow us to evaluate the accuracy and sensitivity of the current generation of implicit solvent models in computing binding affinities on an experimentally well-characterized set of protein systems for which high-quality calorimetric or binding data is available.  
Comparison with explicit solvent binding free energies for the same forcefield will provide valuable information guiding the development of a more accurate generation of models.  
Additional experiments will probe the sensitivity of computed binding affinities to various physical phenomena present in real biological systems critical for pharmaceutical design or biologically relevant interactions, such as protonation state variation at constant pH.  
Finally, we will evaluate advanced sampling methods, including replica exchange among alchemical states, expanded ensembles methods for sampling over chemical species, and mixtures of Monte Carlo and molecular dynamics moves in a Markov chain Monte Carlo framework for further enhancing sampling and thus reducing the computational expense of these simulations.
In addition to reporting the results of these studies in publications, these algorithms will be made available to the modeling and commercial pharmaceutical community through YANK as free, open-source (GPL) software.

{\bf The collaboratory effort we have assembled for this project is uniquely suited to carry out these experiments, and will provide numerous benefits to the field of computational chemistry and biology both in the discoveries that will be made and the tools produced as a result of this study.}  
JDC, MRS, and DLM have extensive experience with alchemical free energy calculations and large-scale / high-performance computing, JDC and MRS with enhanced sampling algorithms and statistically efficient analysis of datasets, MRS with the design of efficient alchemical schemes, and KMB with extensive experience in industrial computational drug design and access to medicinally relevant datasets.
Both the datasets and the software will be released into the public domain.  
VSP is the primary PI for OpenMM, the GPU-enabled physics based simulation toolkit.  
Between the primary participants (JDC and the groups of MRS and DLM), we have sufficient manpower to exploit the resources provided by NCSA Lincoln for these aims.  
{\bf A Teragrid allocation on Lincoln is particularly important for this project because GPU computing resources of this magnitude are unavailable at a university level.}

\required{Relevant peer review}

The OpenMM library (on which YANK is built) is produced by Simbios, the National NIH Center for Biomedical Computing focusing on Physics-based Simulation of Biological Structures, a National Center for Biomedical Computing (NCBC), part of the NIH Roadmap for Medical Research (established in response to RFA-RM-04-003).  
OpenMM implementation details and early benchmarks have been described in a peer-reviewed publication~\citep{friedrichs:j-comput-chem:2009:openmm}.  
The alchemical free energy protocols discussed here have been applied to various systems in explicit solvent by the PI and co-PIs in peer-reviewed publications addressing forcefield accuracy~\citep{shirts_extremely_2003,shirts_alchemical_2007,shirts_solvation_2005,jayachandran_parallelized-over-parts_2006,mobley_confine-and-release_2007,mobley_predicting_2007,mobley_use_2006}.
Components of this proposal have been submitted in proposals to the NSF for funding by MRS and DLM.

\required{Background}

%A detailed understanding of the fundamental physical interactions and statistical mechanics governing the binding of small molecules to biological macromolecules would be of great utility in numerous areas of biophysics, biochemistry, and pharmaceutical science.
%While much is known about molecular interactions, a general tool for reconstituting this understanding to predict ligand binding affinities demonstrates both the shortcomings in this understanding and the gap between understanding and practice.

Alchemical binding free energy calculations, in which the free energy difference between bound and unbound states is estimated by way of sampling a number of nonphysical (\emph{alchemical}) intermediate states, are an efficient way to compute binding free energies rigorously~\citep{shirts_alchemical_2007}.
Quantitative free energy calculations using \emph{explicit} solvent models have been possible for the last decade, but they are expensive and frustrated by slow conformational relaxation timescales.  
Such simulations are still too slow for rational ligand design in practical workflows; but they do show predictive promise~\citep{mobley_predicting_2007,mobley:jmb:2009:lysozyme-m102q}.  
While explicit solvent simulations treat the water atomistically, which is critical for reproducing some phenomena~\citep{mobley_charge_2008}, they are still limited in how realistically they can model a biological system: protonation state and counterion number are fixed, despite their biological counterparts existing in solutions of constant pH and ion concentration.

\begin{wrapfigure}{L}{0.4\textwidth}
\includegraphics[width=0.38\textwidth]{figures/force-rms-vs-cutoff}
\caption{\small \label{figure:error-vs-cutoff} RMS error in the force as a function of cutoff used in Onufriev-Bashford-Case (OBC) GBSA model I (for a static structure of T4 lysozyme L99A.}
\end{wrapfigure}

There has been significant recent development in implicit solvent models; measured by small molecule hydration free energies, they do almost as well as explicit solvent simulations compared to experiment while being considerably less expensive computationally~\citep{mobley_treating_2008}.  
Extensive computational work has demonstrating that protein structure prediction and folding can also be modeled with implicit solvent simulations.  
Implicit solvent models are particularly appealing for \emph{free energy} calculations for numerous reasons:
Water degrees of freedom need not be explicitly sampled; fewer alchemical states are required; enhanced sampling schemes, such as replica exchange and the mixing of Monte Carlo moves (e.g.~sidechain flips) with dynamics, give larger sampling boosts due to the smaller particle count compared to explicit solvent; constant-pH simulation is straightforward because the solvent responds infinitely fast~\citep{mongan_constant_2004}; and the presence of counterions can be modeled implicitly using Debye-H\"{u}ckel-like screening ~\citep{case:theor-chem-acc:1999:gb-salt-screening}, facilitating binding calculations for charged ligands.
Despite these obvious benefits, their accuracy for binding affinities of ligands to proteins is largely unknown.
Few examples of true free binding free energy calculations in implicit solvent have been reported in the literature~\citep{lee_calculation_2006,michel_protein-ligand_2006,gan_binding_2008}.

One difficulty with the current generation of the most popular class of implicit solvent models---the generalized Born models---is that their CPU implementations scale poorly to proteins of pharmacologically relevant size, such as those considered here (Table~\ref{table:systems}).
This is especially true when they are employed \emph{without} cutoffs, which is likely necessary to provide consistent and accurate treatment of solvation (see Figure~\ref{figure:error-vs-cutoff}).  
As a result, there has been a proliferation of fast post-docking scoring methods, such as MM-(PB/GB)SA
~\citep{kollman.accchemres.2000.calculating,brown.jcim.2006.high-throughput} in which differences in \emph{enthalpies} are estimated from simulations of chemical endpoints, and the entropy contribution is inexpensively approximated, but these methods have been of questionable accuracy~\citep{shirts-mobley-brown}.  
%From a purely statistical standpoint, free energy differences are fundamentally easier to converge to a target statistical precision than enthalpy differences.
%Although important details of solvent-solute interactions are averaged out in implicit solvent calculations, such increases in speed might be worth the tradeoff if this allows for more converged conformational sampling, and if the results are understood as a tool for exploring and proposing new chemical intermediates.  
However, due to lack of data, we do not yet know if direct implicit solvent free energy methods with improved sampling will offer significantly improved accuracy over more approximate methods.

Facilitated by advances in GPU computing and the OpenMM library, we are now in a position where implicit free energy calculations are fast enough to be eminently practical---consuming less than 1 wall clock day on a relatively small number of GPUs.
The remaining question is whether they are \emph{sufficiently accurate to be useful}.
Even if they are not, the experiments outlined here will provide a great deal of information about how to improve these models, as well as demonstrating the sensitivity to model parameters and to changes in the model due to physical effects like protonation states~\citep{mobley-dill:structure:2009:review}.

\section*{Computational methodology}

\subsection*{OpenMM}

Preparing GPU code from scratch is an arduous task, requiring great effort on the part of the programmer that may not always result in substantial speedup over CPU performance.  
We are therefore taking advantage of OpenMM~[\url{http://www.simtk.org/home/openmm}], a library intended to simplify and accelerate codes employing molecular mechanics, currently under active development by Simbios, the National NIH Center for Biomedical Computing focusing on Physics-based Simulation of Biological Structures~[\url{http://simbios.stanford.edu/summary.htm}].  
OpenMM is licensed as LGPL, and is hosted on the open-source collaboratory repository SimTK~[\url{http://simtk.org}].  
While a first official release of OpenMM is planned for Q1 2010, we have been working with the developers through Vijay Pande, the principal investigator in the Simbios effort behind OpenMM.

\begin{wraptable}{L}{0.48\textwidth}
\small
{\bf OpenMM benchmarks on NCSA Lincoln} \\
\begin{tabular}{|l|l|l|l|l|} \hline
& & \multicolumn{2}{c|}{ns/day} & \\
system & residues & CPU & GPU  & speedup \\ \hline
FKBP & 107 & 0.60 & 130 & 217x \\
T4 lysozyme & 162 & 0.26 & 60.1 & 231x \\
CDK2 & 297 & 0.077 & 19.6 & 255x \\ 
p38  & 349 & 0.057 & 14.6 & 261x \\ \hline
\end{tabular}
\caption{\label{table:openmm-benchmarks}\small Comparison of OpenMM simulation speed on one NCSA Lincoln Tesla GPU and Amber {\tt sander} on one Lincoln CPU core for some systems listed in Table~\ref{table:systems}.  The benchmark used Amber ff96 with OBC GB/SA ({\tt sander} igb=2) and no cutoffs.}
\end{wraptable}

OpenMM provides a simplified API with plugin architecture for common molecular mechanics operations, including energy and force calculations as well as common molecular dynamics integrators and thermostats.  
OpenMM's API is forcefield-agnostic, as the interface for defining molecular mechanics systems support multiple popular forcefields (e.g.~AMBER, CHARMM, OPLS, and GROMOS).
A CPU-only reference implementation, which is extensively validated against other molecular dynamics codes and simple test systems, emphasizes clear, straightforward implementation of the algorithms to provide a "gold standard" against which optimized implementations can be checked.
Currently, a high-performance CUDA implementation provided, with a OpenCL implementation is planned for a subsequent Q3 2010 release (see the OpenMM Roadmap at \url{http://wiki.simtk.org/openmm/RoadmapTimeline}).

The CUDA implementation, which we make use of here, is able to achieve high throughput on CUDA-capable NVIDIA GPUs, especially the GT200-class GPUs present in the Tesla S1070 accelerator units in NCSA Lincoln.  
For systems of the size relevant to the ligand binding studies here (see Table \ref{table:systems}), dynamics can be simulated on a Tesla GPU up to 261x faster than AMBER's {\tt sander} on a single CPU core for equivalent run configurations (Table \ref{table:openmm-benchmarks}).
A detailed discussion of the way in which OpenMM efficiently implements force computations as CUDA kernels is beyond the scope of this proposal, but details are provided in~\citep{friedrichs:j-comput-chem:2009:openmm}.  
One factor contributing to high throughput is ensuring the dynamical integration scheme is handled on the GPU without the need to transfer coordinates back and forth to the host CPU every timestep.

In addition to the native C++ interface, a Python interface ({\tt pyopenmm}) has recently been developed [\url{http://www.simtk.org/home/pyopenmm}].
{\tt pyopenmm} interoperates with the numpy library for efficient array and matrix computations, and uses a powerful physical unit-handling system to guard against programming errors.  
{\tt pyopenmm} allows developers to take advantage of the flexibility, and expressiveness, and large package library of Python to assemble scientific computing applications, while still retaining the benefits of OpenMM's high-performance GPU-based molecular mechanics operations.  
%Provided the majority of execution time is spent in OpenMM calls of substantial duration (as in the applications here), almost no additional overhead is introduced by constructing applications in Python instead of C++, with the benefit of considerable reductions in development time.

\subsection*{YANK}

\begin{wrapfigure}[27]{L}{0.45\textwidth}
\includegraphics[width=0.45\textwidth]{figures/yank-lysozyme-obcgbsa} 
%\includegraphics[width=0.48\textwidth]{figures/mobley-lysozyme-explicit}
\caption{\small \label{figure:yank-preliminary-data} Estimated free energies from 5 ns/replica Hamiltonian-exchange free energy calculations using the C++ version of YANK. replica exchange diagnostics showed the benzene calculation failed to converge; all other systems mixed well. 
Published explicit-solvent data~\citep{mobley_predicting_2007-1} for this system using the same forcefield had a coefficient of determination $R^2 = 0.61$ and rank correlation coefficient $\rho = 0.75$.} 
\end{wrapfigure}

The code used in the experiments proposed here is called YANK~\footnote{Why the name `YANK'?  In 2008, as a postdoc in the Pande lab, Kim Branson (an Australian) was completing the first version of his fast docking code based on empirical scoring functions, which he called `Ocker' (a dated slang term for 'Australian').   At the same time, John Chodera (an American), also a postdoc in the Pande lab, started work on a complementary code for rigorous free energy calculations, which he naturally called `Yank'.  The name stuck.}, and was built using the OpenMM toolkit both for efficiency and to allow us to focus on higher-level aspects of efficient binding free energy calculations instead of low-level implementation details.  Because the simplified API provided by OpenMM drastically reduces the complexity of codes using molecular mechanics, we have been able to rapidly
assemble a first prototype in C++ that utilizes the Generalized Born solvent models implemented in OpenMM, which we have tested extensively on NCSA Lincoln.
This prototype is only 30\% slower (in terms of ns/replica/day) than the peak OpenMM benchmarks listed in Table~\ref{table:openmm-benchmarks}, despite the modifications that have been made for free energy calculations.  These include introducing soft-core Lennard-Jones and generalized Born interactions, the overhead of transferring coordinates to and from the GPU every iteration and for the energy evaluations following an iteration. We expect to eliminate some of this difference as further optimizations are made.

We are currently working on the second generation of YANK, which is written in Python and employs the {\tt pyopenmm} Python interface to OpenMM.  Because nearly all of the execution time is spent in highly-optimized OpenMM library calls, moving to Python sacrifices little performance while gaining substantial flexibility and easy of implementation of new algorithms. 
In addition, simplifies the use of GPUs attached to separate host nodes in a single calculation through the use of preexisting high-level interfaces such as the 'multiprocessing' and 'parallel python' modules.

A calculation with YANK proceeds as follows:
\begin{itemize}
\item The protein and ligand are read from disk and parameters from the appropriate forcefield(s) are assigned.
\item Alchemical intermediates are created, in which interactions between the ligand and its environment are perturbed, spanning the fully interacting complex to a state in which the ligand had been `annihilated'~\citep{shirts_alchemical_2007}.
The the same is done for ligand alone in solvent in a separate calculation to complete a thermodynamic cycle equivalent to the binding process.  
A soft restraint in the vicinity of the binding pocket is introduced before the protein-ligand interactions are modulated to prevent the ligand from diffusing too far from the protein when in a noninteracting state and to allow the computed binding affinity to be corrected for the standard state.
\item A replica exchange simulation (one replica per alchemical intermediate) is conducted among the different alchemical states in order to allow enhanced exploration of configuration space.
The ligand undergoes stochastic transitions between strongly-interacting and weakly-interacting alchemical states, allowing it to unbind and reorient, and facilitating relaxation in the receptor, critical to reproducing accurate binding affinities~\citep{mobley_use_2006,mobley_confine-and-release_2007}.
The time between exchange attempts (generally $\sim$ 0.5--5 ps) is tuned to balance the desire to obtain maximum throughput (in terms of ns/replica) with enhanced mixing via exchanges.
Between exchange attempts, the replicas are propagated by Langevin dynamics simulations independently, allowing the calculation to be trivially parallelized up to the number of replicas; stalling is avoided if the number of GPUs is an integral divisor of the number of replicas.
\item After each iteration, energies for all replicas are evaluated at all alchemical intermediates, and coordinates and energies are written to a portable NetCDF file that can be analyzed by idle CPU threads while execution proceeds.  
\item The multistate Bennett acceptance ratio (MBAR) method (a provably optimal method for extracting information from the simulation data to estimate free energies)~\citep{shirts_statistically_2008} is used to compute the free energy of binding and error bars from the energies in the NetCDF file, and can provide an online estimate as the simulation proceeds.
\end{itemize}

The new Python version of YANK adopts a modular framework, allowing for additional enhanced sampling algorithms to be easily "plugged in." Planned extensions in 2009 which will be utilized in the experiments described below include the constant-pH capability for implicit
solvent calculations generalized to all forcefields~\citep{mongan_constant_2004}, the addition of Monte Carlo moves for protein degrees
of freedom \citep{smart-darting-mc}, and Monte Carlo moves among chemical compounds through the use of expanded ensembles simulation~\citep{lyubartsev_new_1992,park_comparison_2008} instead of replica exchange.  
Expanded ensembles simulations, in particular, will allow combinatorially large chemical spaces to be sampled over to discover the strongest binders.  
In addition, with optimized PME support (estimated to be twice as slow as simulations in reaction field) available in the Q1 2010 OpenMM release, we will extend some of these techniques for use in explicit solvent calculations, when sensible.

%Primary development of YANK was carried out on a Mac Pro with a NVIDIA GeForce GTX 285 card (which contains a GT200 GPU with 1 GB memory, nearly equivalent to the GPUs in the Tesla S1070 accelerator units).
%Although code development requires access to only a limited number of GPUs, a significantly larger resource is required to test the applicability to the ligand binding, especially the sensitivity of results to choice of implicit solvent method and pH.

A startup allocation on NCSA Lincoln has allowed us to run a number of preliminary tests on the first C++ version of YANK in order to assess the overall throughput, and the computational requirements to perform the test described in this proposal.  
%, and to identify whether inefficiencies would
%be introduced by moving to a Python codebase for components which are
%not time critical.
% MRS: I don't see this data anywhere 
In addition, we have been able to run a first set of tests to assess the accuracy of one of the GB models available in the current version of OpenMM [OBC GBSA~\citep{onufriev_modification_2000,onufriev_exploring_2004}] for computing binding affinities to T4 lysozyme L99A model system (See Fig.~\ref{figure:yank-preliminary-data}). 
Despite an offset (which may be related to a surface area treatment of the nonpolar effect, rather than a dispersion integral treatment), the coefficient of determination $R^2$ and rank correlation coefficient $\rho$ (0.62 and 0.72) were almost identical to previous results obtained with explicit solvent (0.61 and 0.75, respectively), suggesting that the current approach in many cases may be as accurate as explicit water solvation calculations.

%* Preliminary results / report on what we've already done
%  * Used startup allocation for testing on multiple GPUs / development / eliminating bugs / benchmarking / preliminary tests
%     * Figure: Accuracy for T4 lysozyme L99A system in GB, explicit solvent, and experiment, 
%     * Table: Scaling of YANK over multiple GPUs
%     ? Figure: Free energies are easier to compute and more accurate for predicting affinity than MM-GB/PBSA rescoring methods
%     * Figure: Why do we need to run without cutoffs? Accuracy of OBC GBSA vs cutoff distance.

\section*{Proposed Experiments and Justification of Resources}
%
%* Table/Box: Forcefields and GB models to be evaluated
%  * AMBER forcefields ff99sb and ff03, except when ff96 necessary to compare to explicit solvent data; GAFF for ligand with AM1-BCC charges
%    * OBC GBSA
%    * GB/VI
%    * Other GB models that fit CustomGBForce: What models fit the general two-pass structure?  At least all of the models in AMBER.
%  * OPLS; Maestro-based parameter assignment
%    * SGB
%    * AGBNP

\subsection*{Evaluation of generalized Born models in computing ligand binding free energies}

\begin{wraptable}[15]{L}{0.45\textwidth}
\small
\begin{tabular}{|l|c|c|} \hline
Receptor & Residues & Ligands \\ \hline
FKBP & 107 & 40 \\
IL-2 & 132 & 11 \\
T4 lysozyme L99A & 162 & 18 \\
T4 lysozyme L99A/M102Q & 162 & 21 \\
trypsin & 223 & 14 \\
CDK2 & 297 & 30 \\ 
p38 & 349 & 19 \\ \hline
$\beta$-secretase & 298 & 30 \\
thrombin & 392 & 12 \\ \hline
\end{tabular}
\caption{\label{table:systems} {\small {\bf Systems used for evaluation in this study.} This curated set will be made available online at \url{https://simtk.org/home/bindtest/}}}
\end{wraptable}

A number of second- and third-generation generalized Born (GB) models are now in common use for both biomolecular simulation and
post-docking rescoring evaluation. 
These models vary in one or more ways: the philosophy in how the effective Born radii are computed; the functional form used in computing pairwise interactions; their treatments of nonpolar contributions to the solvation energetics; and the parameters they use to approximate these contributions~\citep{onufriev_modification_2000,onufriev_exploring_2004,GBVI,mongan_gb_2006,agbnp2}. 
While the current generation of GB models perform similarly in in computing small-molecule hydration free energies, data for which the models are typically parameterized, their accuracy in modeling the energetics of larger, more complex systems is largely unknown.

We will test a number of GB models and corresponding forcefields to evaluate for their accuracy in describing ligand binding free energies
using rigorous alchemical free energy calculations.  
We divide the tests into a seven protein sensitivity test (listed in Table~\ref{table:GBSA_request} where we test two models, and a four protein set with smaller proteins, where we will test six or models.
Each target contains a number of ligands spanning a range of binding affinities, including non- and weak binders.  
For the full protein set, we will test the Amber {\em ff99sb} forcefield with both the OBC GBSA~\citep{onufriev_modification_2000,onufriev_exploring_2004} model I (equivalent to the {\tt sander} option igb=2) and the GB/VI~\citep{GBVI} model which features a London dispersion integral instead of a surface area term.
For the smaller set, we will test the combinations of Amber {\em ff96} with OBC, {\em ff03} with both OBC and GB/VI, as well as {\em ff99sb}
with a OBC model II (equivalent to {\tt sander}'s igb=5), and the GB$n$ model (equivalent to igb=7)~\citep{mongan_gb_2006}.  
Poisson-Boltzmann (PB) treatments are not currently available on the GPU architecture, but the prospect of postprocessing samples collected from GBSA simulations to estimate free energies with PB models will be explored.  
OBC and GBVI are currently implemented in OpenMM; the other GB models described here will be implemented \emph{via} the CustomGBForce two-pass framework to be implemented following the Q1 2010 public OpenMM release.

Since YANK can support a number of force field parameters, we will also test the sensitivity of the results to both AMBER ff99sb with parameters generated from GAFF and OPLS parameters with ligand parameters. 
Both parameter sets are widely used in protein and small molecule modeling with relatively good accuracy, but a large scale test of sensitivities to ligand binding will illustrate any differences that might arise for this orthogonal test of thermodynamics.  
We will test OPLS-AA with the AGBNP-2 implicit solvent model~\citep{agbnp2}.

We will compare these results with experimental binding affinities
from high-quality isothermal titration experiments to evaluate the
absolute accuracy of the current generation of solvent models, and the
utility of rapid binding affinity computation for rational ligand
design.  Comparison with explicit solvent results for the same
forcefield will allow us to isolate which contribution originate from
the difference in solvent model.  Examination of related implicit
solvent models that differ primarily by reparameterization will
provide insight as to the sensitivity of the computed free energies to
parameters when the functional form is held fixed.  Finally, comparing
different classes of implicit solvent models (in particular, those
that employ surface area-based nonpolar components and those that
employ dispersion integral-based nonpolar components) will inform us as
to which model best captures the physics of polar and nonpolar
solvation in ligand binding.

\begin{wraptable}[14]{r}{0.65\textwidth}
\small
\begin{tabular}{|lcccc|} 
\hline
         & GPU days/ & &full set & small proteins \\
Receptor & ligand & ligands  & GPU days & GPU days\\ 
\hline
FKBP & 5.5 & 40 & 440 & 1320 \\
IL-2 & 10$^{\dagger}$ & 11 & 220 & 660\\
T4L L99A & 12 & 18 & 216 &\\
T4L L99A/M102Q & 12 & 21 & 252 & 1512\\
trypsin & 25$^{\dagger}$ & 14 & 700& 2100\\
CDK2 & 37 & 80 & 740 &\\ 
p38 & 45 & 19 & 855 &\\ 
\hline 
Total & & & 430K & 540K\\
\hline
\end{tabular}
\caption{\label{table:GBSA_request} {\small {\bf Computational requirements for GBSA test sets}. $\dagger$ indicates an estimated timing based on the other benchmarks.}}
\end{wraptable}


Pilot calculations of 2--10 ns/replica on NCSA Lincoln using T4 lysozyme L99A showed strong relaxation behavior over the first 1 ns/replica before a (quasi)equilibrium state is reached, so we plan to run at least twenty times this relaxation time, 20 ns/replica, to obtain reliable statistics.  
To determine whether there are remaining issues with broken ergodicity and slow conformational exchange not sampled in this timeframe, where co-crystal structures exist of the ligand bound to the receptor, we will initiate separate simulations from the co-crystal structure and a variety of ligand poses docked into the \emph{apo} receptor.  
Agreement (within statistical error) would suggest the relevant conformational space has been sampled, while discrepancies will indicate slow relaxation behavior with timescale greater than 20 ns/replica.

Assuming the same 18 Hamiltonian replicas used in the initial data set, with 20 ns/replica, and 2 initial structures per ligand, this would require 720 ns/system.  
As there are 2 GPU's per 8 core node, this translates into 96 SU per GPU/day, {\bf so we request 430K SU for the full test set and 540K SU for the smaller protein test sets, according to the benchmarks shown in Table~\ref{table:GBSA_request}.}


\subsection*{Comparison of free energies with MM-PB/GBSA and other rescoring methods}

MM-PB/GBSA methods, in which the {\em absolute} free energy of a protein ligand complex and its individual components are estimated
from one or a few simulations have enjoyed widespread popularity due to their ease of application, despite mixed performance results \citep{shirts-mobley-brown}, and observations that these simulations must generally be extremely long to provide statistically reproducible results \citep{genheden-ryde:j-comput-chem:2009:mm-gbsa}.  
In statistical mechanical terms, these methods directly estimate an enthalpy difference between
the bound and unbound states, adding an entropy correction based on
one of several crude schemes.  We propose to conduct a thorough study of MM-GBSA methods by comparing
statistically well-converged, reproducible MM-GBSA estimates with
rigorous free energies computed by alchemical means using the {\em same}
GB model and forcefield as a control, to examine the magnitude and
type of errors these methods introduce and judge whether they are
suitable for actual use.  Our hypothesis, which will be borne out or
refuted by this study, is that (1) MM-GBSA methods inherently
introduce substantial errors into the calculation that cannot be
eliminated by longer simulation time, independent of the choice of
solvent model, and (2) estimating of binding free energies through
alchemical means is actually less time-consuming than obtaining a
reliable estimate via MM-PBSA for statistical reasons.

MM-PBSA requires multiple long simulations for each of the targets and
compounds: to collect a typical 3-fold replicate of each ligand
with 100 ns per compound, the timing estimates in
Table~\ref{table:GBSA_request}, {\bf would require 1550 GPU-days or 150K SUs}.

\subsection*{Examination of the importance of constant-pH in computing binding affinities} 

Real macromolecules in biological or experimental contexts are invariably immersed in a buffered solution of nearly fixed pH.  
As a result, a protein does not simply exist in a single protonation state; rather, a distribution of protonation states are sampled, and can vary
with conformation of in the presence of a ligand.  
In addition, ligands often have titratable groups, allowing the protonation state of the ligand to change with binding.  
In some cases, these effects will clearly be relevant (if not dominant) in the observed binding affinity~\citep{klebe:jmb:2001:thrombin}, but in many other cases, it is not clear how critically the binding affinity depends on protonation state.

We will implement a generalized form of the constant-pH simulation procedure for implicit solvent described in Mongan~\citep{mongan_constant_2004}, adding the capacity automatically recompute the required reference free energies for different forcefields.  
We propose to study in detail two systems (thromin and $\beta$-secretase) in which protonation state changes are known to play a critical role in the binding of some ligands though careful ITC experiments~\citep{klebe:jmb:2001:thrombin,merck:bio-med-chem-lett:2009:beta-secretase}.  
In addition, we will examine a third system (trypsin) in which protonation state changes are not known to play a role in the majority of ligands, but
with an eye toward discerning how critical variable protonation states are in the computed binding affinities.  
Conducting these simulations within the replica exchange framework presents no additional challenges: A step in which a number of Monte Carlo titration steps are carried out prior to exchange is added, and free energies are easily computed afterwards with the multistate Bennett acceptance ratio (MBAR) framework which provides tools for analysis of simulations at constant chemical potential~\citep{shirts_statistically_2008}.

We propose to evaluate all the ligands from these three systems listed in Table~\ref{table:systems} with one pair of 20 ns/replica replica exchange binding free energy calculations starting from both \emph{apo} and \emph{holo} structures (as a check of convergence) with constant-pH simulations, to compare to the fixed-protonation-state simulations described in a previous experiment.
One implicit solvent model will be used for this experiment, selected on the basis of its performance in the larger test set.
{\bf This experiment is estimated to require a total of 2180 GPU-days, requiring 210K SUs.}


\subsection*{Calculation of explicit-solvent binding free energies}

While the T4 lysozyme mutant systems and FKBP have been extensively studied by explicit solvent free energy calculations, IL-2, a protein with known allosteric conformational rearrangements, has not.
We propose to compute reference binding free energies for 11 ligands of IL-2 in explicit solvent with YANK using the particle-mesh Ewald (PME) facility of OpenMM using the Amber {\em ff99sb} forcefield, for comparison with the implicit solvent results.
While PME is still undergoing optimization (notably automatic selection of an efficient yet accurate choice of grid spacing), the target performance is only twice as slow as a reaction-field treatment of electrostatics, placing anticipated performance in the range of 9 ns/day for a Tesla GPU.
We therefore estimate that a 20 ns/replica binding free energy calculation for 11 ligands in explicit solvent for IL-2 will require 880 GPU-days, {\bf and request 85K SUs for this purpose}.

\subsection*{ Prospective calculation of binding affinities of 300 fragments screened against trypsin}

While calculations that attempt to reproduce known free energies can be highly informative, it is also critical that calculations are actually {\em predictive}.  
The best way to verify this is to predict binding affinities in advance of their measurement.  
This provides a wholly unbiased assessment of accuracy, and can often reveal unexpected physical effects that lead to further improvement in models \citep{mobley_predicting_2007,mobley:jmb:2009:lysozyme-m102q}.

KMB (collaborator) is part of a consortium in which high-throughput crystallography, surface plasmon resonance, and isothermal titration calorimetry (on a subset of compounds) was used to evaluate 300 soluble compounds from a fragment library for binding to bovine pancreatic trypsin as part of a dataset for the SAMPL2 challenge, a blind prediction effort to be held in 2010 before the data is revealed \citep{sampl2}.  
Later in 2010, in advance of the release of crystallographic or affinity data, we propose to use YANK to compute binding affinities of these 300 compounds in a single implicit solvent model from docked structures as a realistic predictive test for the accuracy with which binding affinities can be predicted.
The implicit solvent model and forcefield will be selected based on the best-performing combination from our earlier studies.

Running replica exchange free energy calculations for a duration 10 ns/replica (which we presume will suffice after efficiency improvements in sampling), these binding calculations are expected to consume 6.25 GPU-days/compound, or 1875 GPU-days for 300 compounds.
{\bf We therefore request a total of 180K SUs for this experiment.}

\subsection*{Testing of methodologies for improved efficiency}

It is likely that even with the improvement in sampling efficiency, we
will require better sampling methods in both conformational and
chemical space.  We will compare sampling efficiency on GPU cores with
direct MD combined with replica exchange in alchemical space. We
propose to include both torsional sampling and "smart darting MC"
steps from Straub along with a sidechain rotamer library
(e.g. Dunbrack) to allow enhanced sampling of sidechain conformational
changes, which present a sampling problem in even the simple T4
lysozyme case~\citep{mobley_confine-and-release_2007}.

\begin{wraptable}[11]{r}{0.35\textwidth}
\small {\bf Total request summary} \\
\begin{tabular}{|lc|} 
\hline
Experiment          & SU \\
\hline
All protein implicit $G_{\mathrm{bind}}$      & 430K \\
Small protein implicit $G_{\mathrm{bind}}$      & 540K \\
MM/GBSA             & 150K \\
Constant pH         & 260K \\
Explicit water IL-2 & 85K \\
Trypsin predictive  & 180K \\
Expanded Ensemble   & 270K \\
\hline 
Total              & 1.9M\\
\hline
\end{tabular}
\end{wraptable}

To truly be useful in drug lead design, even with a throughput of
several ligands per day may be insufficient. To break this barrier, we
will also implement and test methods for multiple end state
simulations, where 5-100 ligands similar ligands with moderate
similarity are run in the same simulation.  With properly optimized
simulations, many intermediate states can be shared, resulting in
significantly fewer total chemical states being required.  If the
number of ligands studied scales beyond the number of replicas, this
can be overcome using expanded ensemble methods, where Monte Carlo
steps between alchemical states are carried out in the same
simulation.~\citep{park_comparison_2008,escobedo_optimized_2007,escobedo_use_2006}. The
number of states required to provide phase space overlap will scale
more slowly than the total number of ligand simulated, making
exploration of this large chemical space feasible.

To effectively test the efficiency of the a number of sampling
improvements will require running several of the smaller comparison
tests multiple times.  Three the smaller of the implicit solvent test
sets is approximately 270K SU.

\subsection*{Request Summary}

To complete the tasks laid out here, we are requesting at total of
1.9M SU on NCSA Lincoln.  This allocation will allow the most
extensive testing of rigorous statistical mechanical ligand
calculations yet, an important test of the applicability of implicit
solvent methods.  In addition, it will help foster the continued
development of computational tools allowing atomistic molecular
simulations to effectively utilize future national computational
resources such as Oak Ridge's recently announced Fermi GPU processor
based supercomputer.

% bibliography
\bibliography{NSFlib}{}




