text:
<!doctype html>

	<meta charset="utf-8">
	<meta http-equiv="X-UA-Compatible" content="IE=edge">
	<title>Manager | mini.ts</title>
	<meta name="description" content="">
	<meta name="viewport" content="width=device-width, initial-scale=1">
	<link rel="stylesheet" href="../assets/css/main.css">

						Search
						<input id="tsd-search-field" type="text" />

Preparing search index...

The search index is not available

					mini.ts

						Options

								All

Public

Public/Protected

All

							<input type="checkbox" id="tsd-filter-inherited" checked />
							Inherited
							<input type="checkbox" id="tsd-filter-externals" checked />
							Externals
							<input type="checkbox" id="tsd-filter-only-exported" />
							Only exported

					Menu

					Globals

					"mini/keyboard"

					Keyboard

					Manager

Class Manager

Hierarchy

						InputEventManager

								Manager

Index

Constructors

constructor

Properties

tracked

Methods

action

activate

de<wbr>Activate

is<wbr>Active

is<wbr>Repeated

is<wbr>Valid

keydown

keyup

no<wbr>Repeat<wbr>Action

set

set<wbr>Repeated

Constructors

constructor

new <wbr>Manager(trackedKeys: KeyCodes[]): Manager

Overrides InputEventManager.constructor

Defined in mini/keyboard.ts:110

Parameters

trackedKeys: KeyCodes[]

Returns Manager

Properties

Protected tracked

tracked: InputEventMap

Inherited from InputEventManager.tracked

Defined in mini/input_event.ts:37

Methods

action

action(trigger: InputID, action: function): void

Inherited from InputEventManager.action

Defined in mini/input_event.ts:72

Calls the action function when a trigger key is active. Will call as long as the key is active

Parameters

trigger: InputID

action: function

(): void

Returns void

Returns void

Protected activate

activate(key: InputID): void

Inherited from InputEventManager.activate

Defined in mini/input_event.ts:81

Parameters

key: InputID

Returns void

Protected de<wbr>Activate

de<wbr>Activate(key: InputID): void

Inherited from InputEventManager.deActivate

Defined in mini/input_event.ts:85

Parameters

key: InputID

Returns void

Protected is<wbr>Active

is<wbr>Active(key: InputID): boolean

Inherited from InputEventManager.isActive

Defined in mini/input_event.ts:39

Parameters

key: InputID

Returns boolean

Protected is<wbr>Repeated

is<wbr>Repeated(key: InputID): boolean

Inherited from InputEventManager.isRepeated

Defined in mini/input_event.ts:43

Parameters

key: InputID

Returns boolean

Protected is<wbr>Valid

is<wbr>Valid(key: InputID): boolean

Inherited from InputEventManager.isValid

Defined in mini/input_event.ts:51

Parameters

key: InputID

Returns boolean

Private keydown

keydown(e: KeyboardEvent): void

Defined in mini/keyboard.ts:104

Parameters

e: KeyboardEvent

Returns void

Private keyup

keyup(e: KeyboardEvent): void

Defined in mini/keyboard.ts:108

Parameters

e: KeyboardEvent

Returns void

no<wbr>Repeat<wbr>Action

no<wbr>Repeat<wbr>Action(trigger: InputID, action: function): void

Inherited from InputEventManager.noRepeatAction

Defined in mini/input_event.ts:58

Calls the action function when a trigger key is active. Must had not been called before

Parameters

trigger: InputID

action: function

(): void

Returns void

Returns void

Protected set

set(key: InputID, value: boolean): void

Inherited from InputEventManager.set

Defined in mini/input_event.ts:90

Parameters

key: InputID

value: boolean

Returns void

Protected set<wbr>Repeated

set<wbr>Repeated(key: InputID): void

Inherited from InputEventManager.setRepeated

Defined in mini/input_event.ts:47

Parameters

key: InputID

Returns void

						Globals

						"mini/keyboard"

								Keyboard

						Key<wbr>Codes

						Manager

								constructor

								tracked

								action

								activate

								de<wbr>Activate

								is<wbr>Active

								is<wbr>Repeated

								is<wbr>Valid

								keydown

								keyup

								no<wbr>Repeat<wbr>Action

								set

								set<wbr>Repeated

						COMBO_<wbr>WASD

Legend

Module

Object literal

Variable

Function

Function with type parameter

Index signature

Type alias

Enumeration

Enumeration member

Property

Method

Interface

Interface with type parameter

Constructor

Property

Method

Index signature

Class

Class with type parameter

Constructor

Property

Method

Accessor

Index signature

Inherited constructor

Inherited property

Inherited method

Inherited accessor

Protected property

Protected method

Protected accessor

Private property

Private method

Private accessor

Static property

Static method

Generated using TypeDoc


text:
\section{Introduction}

Fast and accurate predictions of binding free energies
between proteins and small organic ligands would have significant impact on designing drugs \cite{Michel2010,Chodera2011,Mobley2012a} and other modulators of biological processes.
The clear relevance of protein-ligand binding affinity prediction in chemical biology and drug discovery has inspired a vast array of physics-based methods (for a broad review, see \cite{Gilson2007}), each with a different trade-off between computational accuracy and speed.

On one extreme, molecular docking focuses on speed.
Docking algorithms are designed to quickly obtain plausible configurations of a protein-ligand complex.
Scoring functions are then used to rank one configuration versus another.
Docking programs are commonly assessed by their ability to redock ligands into crystallographic structures from which they have been removed.
Comparative studies \cite{Warren2006,Kim2008} and blinded exercises \cite{Damm-Ganamet2013} consistently show that docking methods are adept at generating the native pose but are less competent at giving it the highest rank.
In this context, it is not surprising that docking scores are poorly correlated with binding free energies \cite{Wang2004,Warren2006,Kim2008,Smith2011}.

In contrast, alchemical pathway methods are based on rigorous statistical mechanics.
The methods involve sampling from a series of possibly nonphysical thermodynamic states
in between end-states where the receptor and ligand are bound and unbound.
In the unbound state, the receptor-ligand nonbonded interaction terms may be switched off or the species may be physically separated.
In accordance to the established statistical mechanics of noncovalent binding \cite{Gilson1997}, the receptor is usually allowed full flexibility.
Unfortunately, sampling a fully flexible complex from multiple statistical distributions along an alchemical pathway generally requires substantial computing resources;
it has been suggested that most published studies are not fully converged \cite{Mobley2012b}!
Even sampling the ligand binding pose is challenging \cite{Mobley2012b}.
To bypass this difficulty, most pathway calculations pursue relative binding free energies between similar molecules and assume that the binding mode does not change.
In absolute binding free energy calculations, ligand sampling issues are usually alleviated by confining the molecule to a specific pose \cite{Boresch2003, Mobley2006}.
Within this restricted range of problems, alchemical pathway calculations are amassing a growing track record of accurate prediction (e.g. \cite{Mobley2007,Michel2008,Boyce2009,Deng2009,Ge2010,Malmstrom2011,Wang2012,Zhu2013,Wang2015}).
\emph{The success of these methods suggests that docking may be substantially improved by incorporating rigorous statistical mechanics.}

Implicit ligand theory (ILT) \cite{Minh2012}, a recently derived statistical mechanics framework for noncovalent association, has the potential to inspire new methods that combine the speed of docking and the rigor of alchemical pathway methods.
ILT formally separates receptor and ligand sampling into two distinct stages.
Because the first stage of receptor sampling does not require a ligand, receptor conformations can be sampled once and used with many different ligands.
In contrast, conventional alchemical pathway methods require thorough receptor sampling for every receptor-ligand pair.
The second stage is to calculate the binding potential of mean force (BPMF) - the binding free energy between a ligand and a \emph{rigid} receptor configuration (Eq. \ref{eq:BPMF}) - for each receptor configuration.
Finally, the standard binding free energy is an exponential average of BPMFs (Eq. \ref{eq:FE_ILT}).

While multiple BPMFs are required to estimate a standard binding free energy, the overall calculation should require less computer time than conventional methods because BPMFs are easier to estimate than binding free energies with a fully flexible receptor. The key reason for this speedup is that a BPMF calculation only requires sampling of the ligand, which usually has many fewer degrees of freedom than the complex. Furthermore, nonbonded interactions between a rigid receptor and ligand can be treated by interpolating precomputed three-dimensional grids, a strategy first developed for docking \citep{Meng1992}. Once the grid is stored, calculation time no longer depends on the size of the receptor. In contrast, conventional alchemical pathway methods require frequent force evaluation between flexible receptor atoms. As the number of pairwise interactions scales as $O(N^2)$ with N receptor atoms (neglecting cutoffs), the relative efficiency of ILT-based methods will be more pronounced as receptor size increases.

In the first paper on ILT, BPMFs were estimated for simple host-guest systems \cite{Minh2012}.
Here, the focus is on precise BPMF estimation for protein-ligand systems.
\citet{Mobley2007, Ucisik2014} previously computed binding free energies between simple ligands and a rigid protein, T4 lysozyme.
In contrast, the present work involves more diverse systems.
\citet{Mobley2007} were equally rigorous, but did not implement specialized methods
to significantly speed their calculations compared to flexible-receptor calculations.
\citet{Ucisik2014} developed a fast method based on a number of approximations.
As in other BPMF \cite{Minh2012} and standard binding free energy \cite{Gallicchio2012, Wang2013b, Wang2015} calculations,
Hamiltonian replica exchange is applied.
The main methodological differences between the current and previous work are the adaptive initialization of thermodynamic states and the use of linearly-scaled interaction grids \cite{Meng1992} with transformation-based smoothing \cite{Venkatachalam2003};
these will be discussed in detail in the section on \emph{Theory and Methods}.
The method is demonstrated on the Astex diverse set \cite{Hartshorn2007},
a curated database of 85 high-quality protein-ligand complexes of pharmaceutical or agrochemical interest.

The algorithm is implemented in a new software package, Alchemical Grid Dock (AlGDock),
a python module based on the Molecular Modeling Toolkit (MMTK) 2.7.8 \cite{Hinsen2000}.
AlGDock is available under the open-source MIT license at \url{https://github.com/ccbatiit/algdock/}.

\section{\label{sec:methods}Theory and Methods}

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

\subsection{Implicit Ligand Theory}

For the noncovalent association between a receptor $R$ and ligand $L$ to form a complex $RL$, $R + L \rightleftharpoons RL$, the standard binding free energy is,
\begin{eqnarray}
\Delta G^\circ = -\beta^{-1} \ln \left( \frac{C^\circ C_{RL}}{C_R C_L} \right),
\end{eqnarray}
where $\beta = (k_B T)^{-1}$ is the inverse of Boltzmann's constant times the temperature, $C^\circ$ is the standard concentration (1 M = 1/1660 \AA$^3$), and $C_X$ is the equilibrium concentration of species $X \in \{R, L, RL\}$. Activities have been assumed to be unity, a reasonable approximation in the limit of low concentrations.

Coordinates of the complex, $r_{RL}$, are partitioned into receptor ($r_R$) and ligand internal ($r_L$) and external ($\xi_L$) coordinates. Based on this partitioning, the interaction energy is defined as $\Psi(r_{RL}) = \mathcal U(r_{RL}) - \mathcal U(r_R) - \mathcal U(r_L)$, where $\mathcal U(r) = U(r) + W(r)$ is an effective potential energy that includes the gas-phase potential energy $U(r)$ and solvation free energy $W(r)$ \cite{Gilson1997}. A BPMF is an exponential average of interaction energies over ligand coordinates in the binding site \cite{Minh2012},
\begin{eqnarray}
B(r_R) & = & -\beta^{-1} \ln \left(
\frac{ \int  I_\xi e^{-\beta \Psi(r_{RL})} e^{-\beta \mathcal U(r_{L}) } ~ dr_L d\xi_L}
       { \int  I_\xi e^{-\beta \mathcal U(r_{L}) } ~ dr_L d\xi_L } \right),
\label{eq:BPMF}
\end{eqnarray}
where $I_\xi \equiv I(\xi_L)$ is an indicator function that takes values between 0 and 1 and specifies whether the receptor and ligand are bound or not.

According to ILT, the standard binding free energy $\Delta G^\circ$ is related to an exponential average of BPMFs over Boltzmann-distributed receptor configurations $r_R$,
\begin{eqnarray}
\Delta G^\circ &=&
-\beta^{-1} \ln \left( \frac{\int e^{-\beta [B(r_R)+  \mathcal U(r_{R})]} dr_R }
{\int e^{-\beta \mathcal U(r_{R})} dr_R } \right) - \beta^{-1} \ln \left( \frac{\Omega C^\circ}{8 \pi^2} \right),
\label{eq:FE_ILT}
\end{eqnarray}
where $\Omega = \int  I_\xi(\xi_L) d\xi_L$ is the volume of the binding site.

\subsection{Thermodynamic Cycle}

A BPMF can also be expressed as a ratio of partition functions \cite{Minh2012},
\begin{eqnarray}
B(r_R) &=& -\beta^{-1} \ln \left(
\frac{ \int  I_\xi e^{-\beta \mathcal U(r_{RL})} ~ dr_L d\xi_L}
       { \int  I_\xi e^{-\beta [\mathcal U(r_{L}) + \mathcal U(r_R)] } ~ dr_L d\xi_L } \right).
\label{eq:BPMF_ratio}
\end{eqnarray}
In AlGDock, BPMFs are estimated by completing the thermodynamic cycle shown in Figure \ref{fig:thermo_cycle}.

\onecolumngrid

\begin{figure}[h]
\begin{center}
\includegraphics[width=3.36in]{fig_thermo_cycle.png}
\end{center}
\caption{
{\bf Thermodynamic cycle for BPMFs: Schematic.} Milestone thermodynamic states are labeled with letters in parentheses. Counter-clockwise from the top left: the non-interacting receptor and ligand are desolvated (a\&b), the temperature is increased to 600 K (b\&c), receptor-ligand interaction grids are scaled from 0 to 1 and the temperature is decreased to 300 K (c\&d), and the complex is solvated and grid energies replaced with direct calculations (d\&e). Arrows with orthogonal lines indicate multiple intermediate thermodynamic states. For BPMF calculations, configurations are sampled from thermodynamic states with the rounded boxes and from their intermediates.
\label{fig:thermo_cycle}}
\end{figure}

Mathematical expressions for individual free energy differences within the cycle are given in Table \ref{tab:b_parts}.
The sum of all the free energy components in Table \ref{tab:b_parts} is $\beta B(r_R)$.

\begin{table}[h]
\caption{
\bf{Thermodynamic cycle for BPMFs: Expressions}}
\begin{tabular}{c c | c c}
Milestones & Process & Expression \\
\hline
a\&b & Receptor desolvation & $f_{ab,R} = -\ln \frac{e^{-\beta_T U(r_R)}}{e^{-\beta_T \mathcal U(r_R)}}$ \\
a\&b & Ligand desolvation & $f_{ab,L} = -\ln \frac{\int I_\xi e^{-\beta_T U(r_L)} dr_L d\xi_L}{\int I_\xi e^{-\beta_T \mathcal U(r_L)} dr_L d\xi_L}$ \\
b\&c & Receptor heating & $f_{bc,R} = -\ln \frac{e^{-\beta_H U(r_R)}}{e^{-\beta_T U(r_R)}}$ \\
b\&c & Ligand heating & $f_{bc,L} =-\ln \frac{\int I_\xi e^{-\beta_H U(r_L)} dr_L d\xi_L}{\int I_\xi e^{-\beta_T U(r_L)} dr_L d\xi_L}$ \\
c\&d & Grid scaling/cooling & $f_{cd,o} =-\ln \frac{\int I_\xi e^{-\beta_T U_g(r_{RL})} dr_L d\xi_L}{\int I_\xi e^{-\beta_H [U(r_L) + U(r_R)]} dr_L d\xi_L}$ \\
d\&e & Complex solvation & $f_{de} = -\ln \frac{\int I_\xi e^{-\beta_T \mathcal U(r_{RL})} dr_L d\xi_L}{\int I_\xi e^{-\beta_T U_g(r_{RL})} dr_L d\xi_L}$ \\
\end{tabular}
\begin{flushleft}
Expressions for free energy differences from the thermodynamic cycle in Figure \ref{fig:thermo_cycle}.
The target temperature $T_T = 300$ K is the basis of $\beta_T = \left(k_B T_T \right)^{-1}$.
Similarly, the high temperature $T_H = 600$ K is the basis of $\beta_H = \left(k_B T_H \right)^{-1}$.
\end{flushleft}
\label{tab:b_parts}
\end{table}

\twocolumngrid

Over the course of the cycle, the receptor-ligand interaction strength is scaled and the temperature is varied.
The temperature is varied because high-temperature states enhance transitions between local energetic minima.
Because the protocol involves temperature changes, the notation is simplified by using a lower-case $u$ to refer to the \emph{reduced potential energy} \cite{Shirts2008}, a log probability density that incorporates $\beta$.
For example, at milestone c, the reduced potential energy is
$u(r_{RL}) = \beta_H \left[ U(r_L) + U(r_R) \right]$.
Likewise, the reduced free energy difference between two milestones $x$ and $y$ is $f_{xy}$.
Because converting from reduced to standard potential energies and free energies involves dividing by $\beta$,
the units of these reduced quantities are specified as $k_B T$.

For thermodynamic states that require $I_\xi$, the ligand is confined to the binding site using a flat-bottom harmonic potential \cite{Gallicchio2012,Wang2013b,Minh2012},
\begin{eqnarray}
	u_I(d) = \left\{
		\begin{array}{lcc}
			0 & \textrm{if} & d \leq d_0 \\
			\frac{1}{2}\beta k(d-d_o)^2 & \textrm{if} & d > d_0 \\
		\end{array}
	\right.,
\end{eqnarray}
where $k = 10000$ kJ/(mol nm$^2$) is the spring constant, $d$ is the distance between the ligand center of mass and the center of the binding site, and $d_0 = 6.0$ \AA~is the radius of the binding site. There is no restriction on ligand rotation.

At milestone d, the reduced potential energy is,
\begin{eqnarray}
u_g(r_{RL}) = \beta_T \left[ U(r_L) + U(r_R) + \Psi_g(r_{RL}) \right],
\end{eqnarray}
where $\Psi_g(r_{RL})$ is the grid interaction energy.
The grid interaction energy,
\begin{eqnarray}
\Psi_g(r_{RL}) = \Psi_{PBSA}(r_{RL}) + \Psi_{vdW}(r_{RL}),
\end{eqnarray}
is based on one electrostatic and two van der Waals grids.

The electrostatic interaction energy $\Psi_{PBSA}(r_{RL})$ is evaluated by multiplying atomic partial charges with the electrostatic potential. The electrostatic potential at each ligand position is obtained by trilinear interpolation of a precomputed grid. The grid is produced by solving the linear Poisson-Boltzmann equation around the minimized receptor molecule using APBS 1.4 \cite{Baker2001} with sequential focusing.
Coarse grids are at least 1.5 times larger than the range of the receptor molecule in each dimension.  Fine grids have the same size as the van der Waals grids, and a spacing of 0.5 \AA.
Coarse grids use multiple Debye-Huckel boundary conditions, and fine grids use coarse-grid solutions as boundary conditions.
Both grids are solved with the following options:
a quintic B-spline charge discretization, spline window width of 0.3, protein dielectric of 2.0, solvent dielectric of 80.0, solvent density of 10.0, solvent radius of 1.4 \AA, smoothed dielectric and ion-accessibility coefficients, and temperature of 300.0 K.

The van der Waals energy $\Psi_{vdW}(r_{RL})$ is evaluated by an analogous grid-based procedure pioneered by \citet{Meng1992}.
In the AMBER \citep{Ponder2003} molecular mechanics force field, the receptor-ligand van der Waals interaction energy is represented as,
$\sum_{i=1}^{N_{lig}} \sum_{j=1}^{N_{rec}} \left[ \frac{A_{ij}}{r_{ij}}^{12} - \frac{B_{ij}}{r_{ij}}^{6} \right]$,
a double sum over ligand atoms $i \in {1, ..., N_{lig}}$ and receptor atoms $j \in {1, ..., N_{rec}}$. $A_{ij}$ and $B_{ij}$ are the repulsion and attraction parameters and $r_{ij}$ is the distance between atoms $i$ and $j$. Molecular docking programs usually model these interactions by a geometric mean approximation \cite{Meng1992}, such that $A_{ij} = \sqrt{A_{ii}} \sqrt{A_{jj}}$ and $B_{ij} = \sqrt{B_{ii}} \sqrt{B_{jj}}$.
With this approximation, the receptor-ligand van der Waals interaction energy at a point $k$ is,
\begin{eqnarray}
\Psi_{vdW}(r_{RL}) &=& \sum_{i=1}^{N_{lig}} \left[ \sqrt{A_{ii}} a_k - \sqrt{B_{ii}} b_k \right] \\
a_{k} &=& \sum_{j=1}^{N_{rec}} \frac{\sqrt{A_{jk}}}{r_{jk}^{12}} \\
b_{k} &=& \sum_{j=1}^{N_{rec}} \frac{\sqrt{B_{jk}}}{r_{jk}^{6}}.
\end{eqnarray}
The values $a_k$ and $b_k$ are precomputed on a three-dimensional grid with a spacing of 0.25 \AA.
Grids span at least $10 + 1.5 d_{max}$ \AA~ in each dimension surrounding the ligand center of mass in the crystallographic binding mode, where $d_{max}$ is the maximum distance from any ligand atom to its center of mass.
For energy and force evaluations, $a_k$ and $b_k$ are approximated at actual ligand atom positions by trilinear interpolation.

Because van der Waals potentials are highly nonlinear, straightforward trilinear interpolation of van der Waals grids provides energies that poorly match directly calculated energies \cite{Oberlin1998}. Hence, van der Waals repulsive grid energies are calculated using a transformation, trilinear interpolation, and inverse transformation procedure \cite{Venkatachalam2003}: If $f(\mathbf{x}_i)$ is the van der Waals attractive or repulsive energy at $\mathbf{x}_i$, a position vector on a grid corner, then $g(x_i) = f(x_i)^{-1/m}$ is evaluated for all grid corners. To estimate $f(\mathbf{x})$ at a desired position vector $\mathbf{x}$, $g(\mathbf{x})$ is obtained by trilinear interpolation of $g(x_i)$ for the eight $x_i$ surrounding $x$, and $f(x) = g(x)^{-m}$. Previously, when equivalent powers were used for the van der Waals repulsive and attractive grids, $m=2$ was found to dramatically reduce interpolation error \cite{Venkatachalam2003}. Here, $m=2$ is used for the repulsive grid and no transformation for the attractive grid; this combination was found to yield the most accurate van der Waals energies for docked poses.

Alchemical transformations that modulate the strength of interactions between atoms often face an ``end-point catastrophe'' where free energy changes are numerically unstable \cite{Michel2010}.
To circumvent this issue, a set of soft Lennard-Jones repulsive and electrostatic grids is introduced between milestones c\&d. In these soft grids, the original grid value, $v_{o}$, is replaced with $v_{max} \tanh \left( v_o/v_{max} \right)$. (\citet{Gallicchio2012} also used a hyperbolic tangent energy cap.) For the soft Lennard-Jones repulsive grid, $v_{max} = 10.0$ kJ$^{1/2}$ mol$^{-1/2}$. To prevent ligands from being ``pinned'' by electrostatics, the soft electrostatic grid uses as maximum value such that the electrostatic energy is less than or equal to the soft Lennard-Jones repulsive energy for every heavy atom at every grid point. This is established by setting $v_{max}$ for the electrostatic grid to 10 times the minimum ratio of Lennard-Jones and electrostatic scaling factors. The reduced potential energy is switched according to the protocol,
\begin{eqnarray}
u_\alpha(r_{RL}) &=& \frac{1}{k_B T(\alpha)} \left[U(r_L) + U(r_R) \nonumber \right. \\
&=& \left. + \alpha_{sg}(\alpha) \Psi_{sg}(r_{RL}) + \alpha_{g}(\alpha) \Psi_g(r_{RL}) \right] \\
\alpha_{sg}(\alpha) &=& -\left(2 \alpha - 1\right)^2 + 1 \nonumber \\
\alpha_{g}(\alpha) &=& \frac{\left(2 \alpha - 1\right)^2}{1 + \exp \left[-1000 (\alpha - \frac{1}{2}) \right]} \nonumber \\
T(\alpha) &=& (T_T - T_H) \alpha + T_H \nonumber
\end{eqnarray}
This protocol turns on the soft grids first, and then the unperturbed grids (Figure \ref{fig:grid_scaling}).
The potential is consistent with milestone c at $\alpha=0$ and milestone d at $\alpha=1$.

\begin{figure}[t]
\begin{center}
\includegraphics[width=3.36in]{fig_grid_scaling-eps-converted-to.pdf}
\end{center}
\caption{
{\bf Grid scaling between milestones c\&d.} $\alpha_{sg}$ (dashed line) and $\alpha_g$ (solid line) as a function of the progress variable $\alpha$.
\label{fig:grid_scaling}}
\end{figure}

\subsection{Sampling and Estimation}

The sampling and estimation strategy in AlGDock is summarized as follows:
\begin{enumerate}
\item \textbf{Ligand preparation:} the ligand is minimized with 500 conjugate gradient steps. The temperature is ramped from 20 K to 300 K over 30 geometrically spaced simulations of 500 steps each.
\item \textbf{Thermodynamic state initialization:} starting from 50 seed configurations, simulations of 1000 steps were used to initialize each thermodynamic state between milestones b\&c and c\&d.
\item \textbf{Hamiltonian replica exchange:} simulations were conducted between milestones b\&c and subsequently between c\&d.
\item \textbf{Postprocessing:} samples from milestones b and d were postprocessed using the generalized Born/surface area (GBSA) water model \cite{Massova2000}. For milestone e, grids were replaced by the rigid receptor configuration and pairwise interactions were directly computed.
\item \textbf{Estimation:} Free energy differences that sum up to the BPMF are estimated.
\end{enumerate}

During thermodynamic state initialization and Hamiltonian replica exchange, the No-U-Turn Sampler (NUTS) \cite{Hoffman2011} is used to propagate the ligand from one configuration to the next.
NUTS is a hybrid molecular dynamics (MD)/Markov chain Monte Carlo (MCMC) sampling technique \cite{Duane1987} that eliminates the need to select the MD trajectory length between Monte Carlo trials.
In other hybrid MD/MCMC techniques, the trajectory length is an adjustable parameter that must be carefully tuned.
NUTS automatically truncates MD trajectories when they make a ``U-Turn'' that brings trajectory end points closer to one another, or when energy changes are too large.
In the present implementation, the change in potential energy or total energy (potential and kinetic energy)
is required to be less than 125 kJ/mol.
The implementation also truncates trajectories so that the allotted number of MD integration steps is not exceeded.

Another adjustable parameter, the MD time step, is selected using dual averaging - an extension to NUTS \cite{Hoffman2011}. Dual averaging adjusts the time step such that an acceptance statistic $H_t$ becomes close to $\delta$. Parameters for dual averaging were $\delta=0.6$, $\gamma=0.05$, $t_0 = 10$, $\kappa = 0.75$, and $\mu = \ln(10 \Delta t)$. Because adjusting the time step disturbs the equilibrium distribution, dual averaging is only used during initialization.

To accelerate transitions between binding poses, external coordinate Markov chain Monte Carlo moves are attempted for thermodynamic states between milestones c\&d for $\alpha < 0.01$. The external coordinate move consists of:
\begin{enumerate}
\item Random rotation. A random quaternion is converted into a matrix that is used to rotate the molecule about its center of mass.
\item Random translation. The magnitude of translation in each dimension is drawn from a Gaussian distribution with a standard deviation of 0.6 \AA.
\end{enumerate}
The move is accepted or rejected according to the Metropolis criterion. Moves are not attempted for $\alpha > 0.01$ due to low acceptance probability.

\subsubsection{Thermodynamic state initialization}

For milestones b\&c,
the first thermodynamic state ($k=0$) is at 300 K and the ligand is steadily warmed to 600 K.
The first thermodynamic state between milestones c\&d depends on whether
there is a fully-bound pose available in the binding site.
If a pose is available, then the first state is fully bound at 300 K.
Otherwise, the first state is fully unbound at 600 K.
In the latter situation, 50 randomly selected configurations from milestone c are placed in the binding site at 453 random center-of-mass positions (0.5 positions per \AA$^3$) and rotated with 100 different random orientations.

After state $k$ is initialized, state $k+1$ is initialized as follows:
\begin{enumerate}
\item \textbf{Parameter selection:} Parameters for state $k+1$ are selected using a new algorithm designed to separate states at approximately even intervals in thermodynamic length.

Thermodynamic length is a metric of the distance on the manifold of thermodynamic states \cite{Weinhold1975}. For a sequence of states, the statistical error in free energy calculations is minimized and the replica exchange frequency is nearly maximized when the intermediate states are equidistant in thermodynamic length \cite{Shenfeld2009}.
Suppose that thermodynamic parameters (e.g. temperature, pressure, grid scaling strengths) are specified by a vector $\mathbf{\lambda}$ with components $\lambda^i$.
Let $\gamma \equiv \gamma(\alpha)$ describe the dependence of $\mathbf{\lambda}$ on the variable $\alpha$, such that $\gamma(0)$ is the initial and $\gamma(1)$ is the final thermodynamic state.
For microscopic systems, the thermodynamic length is defined by the path integral \cite{Crooks2007, Shenfeld2009},
\begin{eqnarray}
\mathcal L & \equiv & \int_0^1 \sqrt{\sum_{i,j} \frac{\partial \gamma^i}{\partial \alpha} g(\gamma)_{ij} \frac{\partial \gamma^j}{\partial \alpha} } d\alpha.
\end{eqnarray}
Given a parameter vector $\lambda$, the reduced potential energy is $u_\lambda(x) = \mathcal U_\lambda(x)/(k_B T_\lambda)$, where $\mathcal U_\lambda(x)$ is the effective potential energy and $T_\lambda$ is the temperature. The normalized log probability of observing a configuration $x$ is $l_\lambda(x) = -u_\lambda(x) - \ln Z_\lambda$, where $Z_\lambda = \int e^{-u_\lambda(x)} dx$ is the partition function. These quantities are used to define elements of the Fisher information matrix,
\begin{eqnarray}
g(\gamma)_{ij} \equiv \sigma^2_\lambda \left[ \partial^i l_\lambda, \partial^j l_\lambda \right],
\end{eqnarray}
where $\sigma^2_\lambda$ is the covariance in state $\lambda$ and $\partial^i$ denotes a partial derivative with respect to $\lambda^i$.
For a protocol in which only one parameter $\lambda^i$ varies with $\alpha$, the length is,
$\mathcal L = \int_0^1 \frac{\partial \gamma^i}{\partial \alpha} \sigma_\lambda \left[ \partial^i l_\lambda \right] dt$.

Numerical estimates of $\mathcal L$ are most accurate when samples are drawn from many states between $0<\alpha<1$ \cite{Crooks2007}. Such exhaustive sampling, however, is unavailable during initialization. A simple approximation for the thermodynamic length when one parameter changes is,
$\mathcal L = \Delta \lambda^i ~ \sigma_0 \left[ \partial^i l_\lambda \right]$,
where $\Delta \lambda^i$ is the total change in the value of the parameter $\lambda^i$ and $\sigma_0 \left[ \partial^i l_\lambda \right]$ is a standard deviation in the initial state. Thus, if one desires $\mathcal L$ to be approximately constant between different intermediate stages in a protocol, then the change in parameter should be inversely proportional to $\sigma_0 \left[ \partial^i l_\lambda \right]$,
\begin{eqnarray}
\Delta \lambda^i = \frac{s}{\sigma_0 \left[ \partial^i l_\lambda \right]},
\end{eqnarray}
where $s$ is an adjustable parameter, the \emph{thermodynamic speed}.

Between milestones $b\&c$, the parameter that varies is the temperature T.
As the log probability of a ligand configuration is
$l_\lambda = -\frac{U(r_L)}{k_B T} - \ln Z_\lambda$,
T is incremented by,
\begin{eqnarray}
\Delta \lambda^i &=& -\frac{s_{bc} k T^2}{\sigma_\lambda[U]},
\end{eqnarray}
where $s_{bc} = 0.1$.
Between milestones $c\&d$,
the log probability of $r_{RL}$ is
$l_\lambda = -u_\alpha(r_{RL}) - \ln Z_\lambda$.
$\alpha$ is incremented by,
\begin{eqnarray}
\Delta \lambda^i = s_{cd} \left[
\abs{\d{\alpha_{sg}}{\alpha}} \frac{\sigma_\lambda[\Psi_{sg}]}{k_B T(\alpha)}
+ \abs{\d{\alpha_{g}}{\alpha}} \frac{\sigma_\lambda[\Psi_{g}]}{k_B T(\alpha)} \right. \nonumber \\
\left. + \abs{T_T - T_H} \frac{\sigma_\lambda[u_\alpha(r_{RL})]}{T(\alpha)} \right]^{-1},
\end{eqnarray}
with $s_{cd} = 0.2$.

If the targeted value of the parameter is exceeded (e.g. temperature increased above 600 K), then the targeted value is used.

\item \textbf{Seed selection:} 50 configurations from state $k$ are resampled as starting seeds for simulations in state $k+1$.

Configurations are drawn from state $k$ with weights proportional to $\exp[u_k(x_i) - u_{k+1}(x_i)]$, where $u_k(x) = \mathcal U_k(x)/(k_B T_k)$ is the reduced potential energy in state $k$. In the limit of infinite sampling of state $k$, resampled configurations would be Boltzmann-distributed in state $k+1$. With imperfect sampling of state $k$, resampled configurations approximate the Boltzmann distribution in state $k+1$.

\item \textbf{Sampling:} Simulations of 1000 steps are run from each seed.

NUTS is run with dual averaging to adjust the time step. Each simulation of 1000 steps has its own estimate of the ideal time step. The median time step, restricted to between 0.25 and 2.5 fs, is subsequently used in Hamiltonian replica exchange.

\item \textbf{Verification:} The mean replica exchange probability, $\E{p_{acc}}$, is used to verify that thermodynamic states are not too distinct nor similar.

It is estimated by taking the sample mean of $p_{acc}$ (Eq. \ref{eq:repXacc}) for every pair of initial samples (at the same time index) from states $k$ and $k+1$. If $\E{p_{acc}}$ is estimated to be too low (below 0.3), then parameters for state $k+1$ are reselected with a smaller increment (thermodynamic speed is adjusted by a factor of 4/5) and simulations are repeated. If it is too high (above 0.99), then state $k$ is removed.
\end{enumerate}

\subsubsection{Hamiltonian replica exchange}

Initialization is followed by Hamiltonian replica exchange \cite{Jiang2010,Gallicchio2012,Wang2013b,Minh2012} between milestones b\&c (15 cycles) and c\&d (18 cycles).
Replica exchange is a Markov chain Monte Carlo move that swaps the configurations of a pair of simulations at different thermodynamic states. (Equivalently, it may be regarded as swapping the states.)
Consider the thermodynamic states $a$ and $b$ with reduced energies $u_a$ and $u_b$, respectively.
If $x$ is the original configuration in state $a$ and $y$ the original configuration in state $b$, then the acceptance probability,
\begin{eqnarray}
p_{acc} &=& \min \left[ 1, e^{-u_a(y)-u_b(x)+u_a(x)+u_b(y)} \right],
\label{eq:repXacc}
\end{eqnarray}
preserves the Boltzmann distribution in both states.

Typical replica exchange protocols attempt exchanges between pairs of neighboring thermodynamic states, but this restriction is unnecessary.
As replica exchange is a type of Gibbs sampling \cite{Chodera2011a}, an arbitrary number of attempts can be made between arbitrary pairs of states.
In the present simulations, each replica exchange cycle consists of 1000 iterations with 50 NUTS steps, 20 external coordinate MCMC moves (if appropriate), and a sweep of replica exchange. Each sweep of replica exchange includes attempts to swap configurations between pairs of states that are $1, 2, ..., \min(5, K)$ states apart, where $K$ is the total number of thermodynamic states in the direction.

The statistical inefficiency $g$ is estimated from the integrated autocorrelation time of replica exchange indices, as described \cite{Chodera2007}.
Multiple samples (up to 3 between milestones b\&c and 25 between milestones c\&d) were saved for every $g$ replica exchange sweeps.

\subsubsection{Estimation}

BPMFs were estimated according to,
\begin{eqnarray}
B(r_R) &=& f_{ab,R} + f_{ab,L} + f_{bc,L} + f_{cd} + f_{de}.
\end{eqnarray}
This expression replaces the sum $f_{bc,R} + f_{cd,o}$ with,
\begin{eqnarray}
f_{cd} &=& -\ln \frac{\int I_\xi e^{-\beta_T [U(r_L) + \Psi_g(r_{RL})]} dr_L d\xi_L}{\int I_\xi e^{-\beta_H U(r_L)} dr_L d\xi_L},
\end{eqnarray}
which can be evaluated without determining the receptor internal energy $U(r_R)$.
The receptor desolvation free energy is estimated by the difference, $f_{ab,R} = \beta_T (U(r_R) - \mathcal U(r_R))$.

Other free energy differences are estimated based on equilibrated samples from replica exchange
between milestones b\&c or c\&d.
Equilibration is determined by calculating the mean and standard deviation of the energy in milestones b or d, respectively.
Systems are considered to be equilibrated once the mean energy is within a standard deviation of the mean energy in the last cycle.
$f_{ab,R}$ and $f_{de}$ are estimated by free energy perturbation \cite{Zwanzig1954} using configurations drawn from milestones a and d, respectively.
$f_{bc,L}$ and $f_{cd}$ are estimated by the multistate Bennett acceptance ratio \cite{Shirts2008}, which uses potential energies from every replica.

\subsection{Astex diverse set BPMF calculations}

The Astex diverse set is a subset of the protein data bank carefully curated for quality and diversity  \cite{Hartshorn2007}.
The members of the set will be referred to be their protein data bank identifiers, e.g. 1s3v.
Each system in the set was set up using AmberTools 14 \cite{Case2012}.
Proteins and ions (from protein.mol2) were parameterized using the AMBER ff14SB force field. Zinc parameters were also used for iron, mercury, and magnesium. Non-standard amino acid residues were converted to their standard counterparts. Using the python wrapper to the OpenEye OEChem Toolkit, all crystallographic ligands were parameterized with the Generalized Amber Force Field \cite{Wang2004a} and AM1BCC charges \cite{Jakalian1999,Jakalian2002}. The main ligands from ligand.mol were regenerated to lose memory of the crystallographic pose. Noncrystallographic atoms in each receptor were minimized for 2500 steps in gas and 2500 steps in GBSA solvent using NAMD 2.9 \cite{Phillips2005}.

For each system, 15 independent simulations were performed between milestones b\&c.
Subsequently, simulations were continued between milestones c\&d from three different starting points:
crystallography, docking, or the lowest-energy pose (according to $u_g$) sampled in all the other simulations.
In most situations, starting from docked poses will be the most practical approach, but the others are a useful comparison.

Docked poses were obtained using the anchor and grow algorithm in UCSF DOCK 6 \cite{Lang2009} to redock ligands into their respective proteins. First, sphgen was run with default parameters: dotlim = 0.0, radmax = 4.0 \AA, and radmin = 1.4 \AA. DOCK was run with 20000 maximum orientations, using internal energy with an exponent of 12, a flexible ligand, an a minimum anchor size of 40. Pruning was performed with clustering, 1000 maximum orientations, a clustering cutoff of 1000, and conformer score cutoff of 25.0. A bump filter was used with max\_bumps\_anchor = 12 and max\_bumps\_growth = 12. Final conformations were clustered with a root mean square deviation (RMSD) threshold of 2.0 \AA. For two systems, 1t46 and 1n46, no docked poses in the binding site were obtained.

After the starting poses were minimized for 1000 conjugate gradient steps in $u_g$, the lowest-energy pose was used to initialize milestone d. After the thermodynamic states were initialized, docked poses within the binding site were also used as starting points for each state in replica exchange. The lowest-energy pose was used in milestone d. Higher-energy poses were used in intermediate replicas to fill all available thermodynamic states. If there were more states than docked poses, the lowest-energy pose was duplicated. For simulations starting from a single pose, replicas were not reinitialized at the beginning of replica exchange.

Calculations were performed on the Open Science Grid \cite{Pordes2007}, Duke shared computing resources, or the Minh group computing cluster at IIT. Each node on the Minh group cluster includes dual AMD Opteron 6376 2.30GHz 16MB cache sixteen-core processors and 64GB DDR3 1600 ECC registered memory. Benchmark calculations were run on the Minh group cluster.

To quantify the precision of the results, we consider the root second moment about a reference value $x_r$,
\begin{eqnarray}
\sigma[x,x_r] = \sqrt{ \frac{1}{N} \sum_{n=1}^N \left( x_n - x_r \right)^2 },
\label{eq:rsmr}
\end{eqnarray}
for $N$ samples indexed by $x_n$. When $x_r$ is the mean, then Equation \ref{eq:rsmr} is the standard deviation, which will simply be referred to as $\sigma[x]$.
We also consider the root second moment about the minimum, $\sigma[x,x_{min}]$, and about the sample mean at the end of a simulation, $\sigma[x,\bar{x}_{end}]$.

\section{Results}

\subsection{Thermodynamic state initialization}

\subsubsection{Protocol length}

Overall, the present method for thermodynamic state initialization robustly yields protocols that are adapted to specific systems.
The large range in the number of thermodynamic states, $N_{states}$, provides evidence of system-specific adaptation (Figure \ref{fig:n_states}).
Between milestones b\&c, the average number of states, $\bar{N}_{states}$ ranges from 35.1 to 93.6.
Between milestones c\&d for simulations initialized from the crystallographic pose,
the $\bar{N}_{states}$ ranges from 84.6 to 190.1.
For most systems, the standard deviation of the number of states, $\sigma[N_{states}]$ is small relative to the $\bar{N}_{states}$, indicating that the procedure is robust.
The robustness and apparent system-specificity of the protocols implies that a universal thermodynamic protocol may strongly deviate from optimality, e.g. $\E{p_{acc}}$ may be inconsistent across replicas.

\begin{figure}[h]
\begin{center}
\includegraphics[width=3.36in]{fig_n_states-eps-converted-to.pdf}
\end{center}
\caption{
{\bf Number of thermodynamic states} (a)
between milestones b\&c, and (b) between milestones c\&d.
The marker indicates the mean value and error bars the standard deviation of 15 independent simulations.
Between milestones c\&d,
simulations are initialized from the crystallographic (x),
lowest-energy docked (circles),
or lowest-energy observed (downward-facing triangles) pose.
They are ordered by the mean number of states
for simulations initialized from the crystallographic pose.
\label{fig:n_states}}
\end{figure}

In some systems, however, $N_{states}$ is not robust; $\bar{N}_{states}$ is dependent on the starting pose and/or $\sigma[N_{states}]$ is large.
Protocols will vary when the initialization stage samples from distinct local minima and estimates different thermodynamic lengths.
Between milestones b\&c, the $\sigma[N_{states}]$ is larger than 2 in three sets of simulations: 1tz8 (6.1), 1r1h (7.3), and 1vcj (13.6).
The larger variance in protocol lengths occurs for different reasons (Figure \ref{fig:cool_protocols}).
In 1r1h simulations, the protocols require variable number of states to reach 400 K,
after which the slope is similar.
For 1tz8 and 1vcj, there are two classes of protocols, indicative of sampling from distinct local minima. With 1tz8, the two protocols are obtained with approximately equal probability. With 1vcj, the short protocol was observed in only 2 of 15 simulations.

\begin{figure}[h]
\begin{center}
\includegraphics[width=3.36in]{fig_cool_protocols-eps-converted-to.pdf}
\end{center}
\caption{
{\bf Protocols between milestones b\&c}
for 15 independent simulations of the ligand
from (a) 1tz8, (b) 1r1h, and (c) 1vcj.
\label{fig:cool_protocols}}
\end{figure}

Protocols between milestones c\&d are more variable than those between milestones b\&c. Interaction grids introduce more possibilities for trapping in local minima. For all three starting points, 1l7f has the largest variance in $N_{states}$. The protocols are most similar for the first 40 thermodynamic states and the last 10, but take different paths in between (Figure \ref{fig:dock_protocols_1l7f}). It is worth noting, however, that high variation in protocols does not necessarily lead to inaccurate results, as will be explained in the following sections.

\begin{figure}[h]
\begin{center}
\includegraphics[width=3.36in]{fig_dock_protocols_1l7f-eps-converted-to.pdf}
\end{center}
\caption{
{\bf Protocols between milestones c\&d}
for 15 independent simulations of the complex from 1l7f,
starting from the crystallographic (a and b),
docked (c and d),
and lowest-energy observed (e and f)
poses.
The left panel shows the protocol on a linear scale for $\alpha > 0.01$,
indexed such that the first state is at $\alpha=1.0$.
The right panel shows the protocol on a semilog scale for $\alpha < 0.01$,
indexed such that the first state is at $\alpha=0$.
\label{fig:dock_protocols_1l7f}}
\end{figure}

\subsubsection{Time steps}

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

\begin{figure}[h]
\begin{center}
\includegraphics[width=3.36in]{fig_time_steps-eps-converted-to.pdf}
\end{center}
\caption{
{\bf Time step statistics.}
Statistics for time steps from dual averaging, based on all protocols from all systems.
Histograms of time steps between 1.0 and 2.0 fs
for (a) between milestones b\&c and (c) between milestones c\&d.
The largest bin count is 10220 between milestones b\&c and 63514 between milestones c\&d.
Five protocols with the (b) shortest and (d) longest
observed time steps between milestones c\&d.
The shortest time steps are observed in 2bm2 (x) and 1n46 (other markers).
The longest time steps are observed in 1jd0.
\label{fig:time_steps}}
\end{figure}

\begin{table}[h]
\caption{
\bf{Time step statistics}}
\begin{tabular}{c | c  c  c  c}
Milestones & Mean & Standard Deviation & Min & Max \\
\hline
b\&c & 1.58 & 0.11 & 1.27 & 2.09 \\
c\&d & 1.55 & 0.12 & 0.87 & 2.08 \\
\end{tabular}
\begin{flushleft}
Statistics for time steps (in fs) from dual averaging, based on all protocols from all systems.
\end{flushleft}
\label{tab:time_steps}
\end{table}

\subsubsection{Replica exchange acceptance probabilities}

A good replica exchange protocol has a reasonable $\E{p_{acc}}$ between all neighboring states.
The key benefit of replica exchange is to spread sampled configurations across a range of different thermodynamic states.
A bottleneck in the exchange of configurations across the pair of states can eliminate this benefit of replica exchange; the simulations effectively become two independent sets of simulations.
Low exchange probabilities are also indicative of poor phase space overlap,
which can limit the convergence of free energy estimates \citep{Lu2001}.

Because of its importance to efficient simulation, AlGDock includes an estimate of $\E{p_{acc}}$ to verify new thermodynamic states during the initialization step. (The notation $\bar{p}_{acc}$ is used to refer to an \emph{estimator} for $\E{p_{acc}}$). In some simulations, however, the phase space explored during replica exchange can be distinct from that explored during initialization, leading to a substantial change in the observed $p_{acc}$. Thus, thermodynamic state verification was followed up by estimating $\E{p_{acc}}$ based on equilibrated samples during replica exchange.

Estimates of $\E{p_{acc}}$ indicate that the initialization protocol does not lead to any simulations with replica exchange bottlenecks (Figure \ref{fig:mean_acc} and Table \ref{tab:mean_acc}).
Between milestones b\&c, $\bar{p}_{acc}$ estimated during replica exchange are high and have low variance. This outcome is consistent with achieving the goal of nearly equal thermodynamic length between adjacent thermodynamic states. It also implies that the configuration spaces explored during initialization and replica exchange are largely the same. Between milestones c\&d, the replica exchange rates are also high but there is a larger variance. While most $\bar{p}_{acc}$ are between 0.7 and 1.0, there are a few simulations where $\bar{p}_{acc}$ is much lower; the lowest observed values are during simulations of 2bm2 (0.21) and 1gm8 (0.28). In these simulations, the acceptance probability dips around $\alpha=0.75$, but are high for most other values of $\alpha$. These drops indicate that different conformations are explored in replica exchange compared to during thermodynamic state initialization. However, even these low $\bar{p}_{acc}$ are not low enough to be considered a bottleneck; they are expected to allow configurations to pass through the thermodynamic states several times per cycle.

\begin{figure}[t]
\begin{center}
\includegraphics[width=3.36in]{fig_mean_acc-eps-converted-to.pdf}
\end{center}
\caption{
{\bf Mean acceptance probability statistics.}
Statistics for mean acceptance probabilities from replica exchange, based on all simulations.
Histograms of $\bar{p}_{acc}$ between 0.7 and 1.0
for (a) between milestones b\&c and (c) between milestones c\&d.
The largest bin count is 35785 between milestones b\&c and 63592 between milestones c\&d.
$\bar{p}_{acc}$ for fifteen protocols with the lowest observed $\bar{p}_{acc}$ (b) between milestones b\&c and (d) and between milestones c\&d are shown with the a line connecting neighboring states.
Between milestones b\&c,
the temperature increases with $\alpha$ such that $T(\alpha) = (T_H - T_T) \alpha + T_T$.
The simulations are from 1s3v (7), 2bm2 (3), 1r55 (3), 1z95 (1), and 1v48 (1).
Between milestones c\&d, the simulations are from 1v48 (10), 1t40 (3), 2bm2 (1), and 1gm8 (1).
\label{fig:mean_acc}}
\end{figure}

\begin{table}[h]
\caption{
\bf{Mean acceptance probability statistics}}
\begin{tabular}{c | c  c  c}
Milestones & Mean & Standard Deviation & Minimum \\
\hline
b\&c & 0.94 & 0.010 & 0.84 \\
c\&d & 0.91 & 0.039 & 0.21 \\
\end{tabular}
\begin{flushleft}
Statistics for mean acceptance probabilities from replica exchange, $\bar{p}_{acc}$,
based on all simulations.
\end{flushleft}
\label{tab:mean_acc}
\end{table}

The observed $\bar{p}_{acc}$ also underscore the point that large variation in protocols is not necessarily problematic.
While simulations of 1tz8, 1r1h, and 1vcj have a relatively large variation in $N_{states}$ between milestones b\&c, all $\bar{p}_{acc}$ estimated from replica exchange are greater than 0.9 (Figure \ref{fig:cool_mean_acc}).
As is particularly clear for low temperature states of 1r1h, longer protocols have a higher $\bar{p}_{acc}$.
Another bump in $\bar{p}_{acc}$ occurs between the last pair of thermodynamic states, where an even interval in thermodynamic length would necessitate going higher than the temperature of interest.

\begin{figure}[t]
\begin{center}
\includegraphics[width=3.36in]{fig_cool_mean_acc-eps-converted-to.pdf}
\end{center}
\caption{
{\bf Mean acceptance probabilities}
for 15 independent simulations of the ligand
from (a) 1tz8, (b) 1r1h, and (c) 1vcj between milestones b\&c.
$\bar{p}_{acc}$ is shown as a line connecting neighboring temperatures.
\label{fig:cool_mean_acc}}
\end{figure}

\subsection{Simulation time}

Due to the heterogeneity of computing resources, variability of queue times, and the diverse nature of the systems, calculations took a variable length of time to complete. Some calculations completed in one day, and all were done within a week.

Total CPU time of benchmark simulations spanned a large range, from 13.6 to 71.1 hours
(Figure \ref{fig:benchmark_timings}).
In all simulations, the majority of time was spent in replica exchange between milestones c\&d, and the second largest fraction in replica exchange between milestones b\&c. Initialization and analysis between milestones c\&d consumes a small but noticeable fraction of simulation time, whereas the analogous steps between milestones b\&c consume a nearly negligible fraction of simulation time.
There is a trend relating the calculation time affiliated with milestones a to c and milestones c to e, but the correlation is imperfect.

~ \\

\begin{figure}[t]
\begin{center}
\includegraphics[width=3.36in]{fig_benchmark_timings-eps-converted-to.pdf}
\end{center}
\caption{
{\bf Benchmark simulation times}
starting from the crystallographic pose.
In ascending order, lines depict the cumulative time for initialization, replica exchange, and analysis (running GBSA energy and MBAR calculations), first for milestones a to c and then milestones c to e.
Systems are ordered along the x axis by the total simulation time.
\label{fig:benchmark_timings}}
\end{figure}

\subsection{Sampling}

At milestone c, sampled configurations are uniformly distributed in a sphere. As $\alpha$ increases from 0 to 1 between milestones c\&d, the phase space of the ligand is gradually restricted. First, the soft grids prevent ligand atoms from overlapping with receptor atoms. At intermediate values of $\alpha$, the ligand may assume several poses (e.g. Figure \ref{fig:dock_intermediates_1l7f}). Finally, at milestone d, the ligand usually but does not always sample from a single minimum.

\onecolumngrid

~ \\

\begin{figure}[h]
\begin{center}
\includegraphics[width=6.68in]{fig_dock_intermediates_1l7f_100dpi-eps-converted-to.pdf}
\end{center}
\caption{
{\bf Samples from evenly spaced thermodynamic states between milestones c\&d},
taken from a representative simulation of 1l7f starting from docked poses.
\label{fig:dock_intermediates_1l7f}}
\end{figure}

In many simulations, the phase space sampled at higher $\alpha$ is a subset of that sampled at lower $\alpha$ (e.g. Figure \ref{fig:dock_intermediates_1kzk}). In others, however, the conformational minima for $\alpha \approx 0.5$, where the soft grids are at full strength, are entirely distinct from the minima for $\alpha \approx 1$ (e.g. Figure \ref{fig:dock_intermediates_1l7f}). When there are shifts in important phase space, the thermodynamic states at $\alpha \approx 0.5$ are less beneficial to sampling from milestone d. Nonetheless, phase space overlap between adjacent states is a sufficient condition for precise free energy estimates.

\begin{figure}[h]
\begin{center}
\includegraphics[width=6.68in]{fig_dock_intermediates_1kzk_100dpi-eps-converted-to.pdf}
\end{center}
\caption{
{\bf Samples from evenly spaced thermodynamic states between milestones c\&d},
taken from a representative simulation of 1kzk starting from docked poses.
The protein structure structure is shown with ribbons and
the crystallographic ligand pose is shown with a thick licorice representation
and purple carbon atoms.
The same illustration scheme is used in Figures \ref{fig:dock_intermediates_1l7f} and \ref{fig:dock_last}.
Figures \ref{fig:dock_intermediates_1kzk}, \ref{fig:dock_intermediates_1l7f}, and \ref{fig:dock_last}
were generated with VMD \cite{HUMP96}.
\label{fig:dock_intermediates_1kzk}}
\end{figure}

\newpage

At milestone d, most simulations sample from a single minimum (Figure \ref{fig:dock_last}),
but there are several other situations. These include,
\begin{enumerate}
\item Sampling several alternate poses that are (e.g. Figure \ref{fig:dock_last}(a))
or are not (e.g. Figure \ref{fig:dock_last}(d)) clearly separated;
\item Sampling a clearly defined but sparsely populated alternative pose (e.g. Figure \ref{fig:dock_last}(b));
\item Sampling alternate poses that share a common warhead position, but have a floppier tail (e.g. Figure \ref{fig:dock_last}(g)).
\end{enumerate}

\begin{figure}[h]
\begin{center}
\includegraphics[width=6.68in]{fig_dock_last_100dpi-eps-converted-to.pdf}
\end{center}
\caption{
{\bf Samples from milestone d}, taken from representative simulations.
\label{fig:dock_last}}
\end{figure}

\twocolumngrid

\newpage

\subsection{Native pose identification}

A large fraction of simulations include the crystallographic pose among samples from milestone d (Figure \ref{fig:rmsds}).
Unsurprisingly, it is most often observed in simulations starting from the crystallographic pose.
After equilibration, all the simulations still sample a pose with a heavy-atom RMSD less than 0.75 \AA.
In contrast, a native pose (defined as having a heavy-atom RMSD less than 2.0 \AA) is observed in 85.6\% of simulations starting from the lowest energy pose and 77.7\% of simulations starting from docked poses. According to $u_g$, the lowest-energy observed pose is a native pose in 61/85 = 71.8\% of systems. This implies that Boltzmann sampling was able to start at another pose and find the native pose in 13.8\% of simulations. Starting docked poses include a native pose in 75/85 = 88.2\%, and rank the native pose with the lowest $u_g$ (and lowest UCSF DOCK 6 grid score) in 48/85 = 56.5\% of systems. This implies that the sampling procedure was able to bring the native pose into milestone d, via Monte Carlo moves or replica exchange, in 21.2\% of simulations.

\begin{figure}[h]
\begin{center}
\includegraphics[width=3.36in]{fig_rmsds-eps-converted-to.pdf}
\end{center}
\caption{
{\bf Heavy-atom RMSDs from the native state,}
for samples from milestone d (after equilibration) from all of the simulations starting from
(a) the crystallographic, (b) docked, and (c) lowest-energy observed poses.
Symbols denote the lowest observed RMSDs in the simulation (x), or
RMSDs of the lowest-energy pose according to $u_g$ (circles) or the GBSA force field (triangles).
Each data set is sorted in ascending order.
To best emphasize native-like poses, the ordinate is truncated at 3.0 \AA.
\label{fig:rmsds}}
\end{figure}

Overall, the GBSA force field is better at identifying the native pose than the grid energy $u_g$.
According to the grid energy,
the lowest energy configuration has a RMSD less than 2.0 \AA~ in 78.5\% of simulation starting from the crystallographic pose, 68.9\% of simulations starting from the lowest energy pose, and 55.6\% of simulations starting from docked poses.
According to the GBSA force field, the analogous fractions are 87.1\%, 75.1\%, and 61.5\% of simulations, respectively.

\subsection{Estimation}

\subsubsection{Standard deviation of free energy estimates}

In most systems, the standard deviation of the BPMF converges to within 5 $k_B T$ (Figure \ref{fig:std_fe} and Table \ref{tab:nconverged}).
Within the thermodynamic cycle, free energies estimated based on multiple intermediate states ($f_{bc}$ and $f_{cd}$) were more precise, and those estimated by single-step perturbation ($f_{ab}$ and $f_{de}$) were less precise.
Simulating from multiple intermediate states promotes phase space overlap between adjacent states, leading to more precise free energy estimates.
Indeed, $\sigma[f_{bc}]$ and $\sigma[f_{cd}]$ are less than 5 $k_B T$ except in two systems with no docked poses in the binding site: 1n46 and 1t46. For these simulations, randomly placed ligands end up sampling different local minima in independent simulations and result in distinct free energy estimates.
Among all free energy differences between adjacent milestones, $f_{de}$ is usually the least precise (Table \ref{tab:least_precise}).

\onecolumngrid

\begin{figure}[h]
\begin{center}
\includegraphics[width=6.68in]{fig_std_fe-eps-converted-to.pdf}
\end{center}
\caption{
{\bf Standard deviation of free energy estimates}
for (a) 15 independent simulations of the ligand or
of the complex starting from (b) the crystallographic,
(c) docked, and (d) lowest-energy observed poses.
Panel (a) is ordered by the standard deviation of the free energy difference between milestones a\&c.
The other panels are ordered by the standard deviation of the free energy difference between milestones a\&e.
\label{fig:std_fe}}
\end{figure}

\twocolumngrid

\begin{table}[h]
\caption{
\bf{Number of converged systems}}
\begin{tabular}{c c | c c | c}
Milestones & Starting point & $<1.5 ~ k_BT$ & $<5.0 ~  k_BT$ & False \\
\hline \hline
a\&b & - & 69 & 82 & - \\
b\&c & - & 82 & 85 & - \\
a\&c & - & 64 & 80 & - \\
\hline
c\&d & xtal                & 84 & 85 & 5 \\
c\&d & dock              & 82 & 83 & 9 \\
c\&d & lowest $u_g$ & 77 & 85 & 4 \\
c\&d & all                  & 48 & 71 & - \\
\hline
d\&e & xtal                & 34 & 76 & 4 \\
d\&e & dock              & 26 & 73 & 8 \\
d\&e & lowest $u_g$ & 37 & 80 & 11 \\
d\&e & all                  & 17 & 63 & - \\
\hline
a\&e & xtal                & 28 & 74 & 6 \\
a\&e & dock              & 20 & 68 & 11 \\
a\&e & lowest $u_g$ & 27 & 77 & 12 \\
a\&e & all                  & 15 & 52 & - \\
\end{tabular}
\begin{flushleft}
Number of systems (out of 85) in which the free energy difference has a standard deviation of less than 1.5 or less than 5.0 $k_BT$. False convergence means that the system $\sigma[x] < 5.0 ~k_BT$, but $\sigma[x,x_{min}] > 15.0 ~ k_B T$. For the xtal, dock, and lowest $u_g$ starting points, the standard deviation is based on 15 independent simulations. The starting point of `all' combines data from the other starting points, and is thus based on 45 independent simulations.
\end{flushleft}
\label{tab:nconverged}
\end{table}

\begin{table}[h]
\caption{
\bf{Least precise free energy components}}
\begin{tabular}{c c | c c c c }
Starting point & Subset & a\&b & b\&c & c\&d & d\&e \\
\hline \hline
xtal                & all & 11 & 0 & 5 & 69 \\
dock              & all & 10 & 1 & 5 & 69 \\
lowest $u_g$ & all & 12 & 1 & 9 & 63 \\
all                  & all & 6 & 0 & 20 & 59 \\
\hline
xtal                & imprecise & 5 & 0 & 0 & 10 \\
dock              & imprecise & 3 & 0 & 2 & 18 \\
lowest $u_g$ & imprecise & 4 & 0 & 1 & 13 \\
all                  & imprecise & 3 & 0 & 11 & 25 \\
\end{tabular}
\begin{flushleft}
The least precise free energy estimates between adjacent milestones. Counts are based on all simulations or the imprecise subset, when the standard deviation of the total BPMF is greater than 5.0.
\end{flushleft}
\label{tab:least_precise}
\end{table}

For systems with poor convergence of BPMF estimates, longer sampling is not necessarily an effective way to improve the precision of free energy estimates.
While $\sigma[f,f_{end}]$ usually decreases monotonically with sample size, this trend is not universally true (Figure \ref{fig:std_fe_cycle}).
In some systems there are temporary jumps in $\sigma[f,f_{end}]$ due to instability in determining the equilibrated cycle.
More importantly, in other systems, $\sigma[f,f_{end}]$ appears flat, with little or no change with increased sampling.

\onecolumngrid

\begin{figure}[h]
\begin{center}
\includegraphics[width=6.68in]{fig_std_fe_cycle-eps-converted-to.pdf}
\end{center}
\caption{
{\bf Standard deviation of free energy estimates as a function of cycle}:
$\sigma[f_{ac}, \bar{f}_{ac, end}]$
based on 15 independent simulations of (a) the ligand or
$\sigma[f_{ae}, \bar{f}_{ae, end}]$
based on the ligand simulations and 15 independent simulations
of the complex starting from the (b) crystallographic,
(c) docked, and (d) lowest-energy observed poses.
\label{fig:std_fe_cycle}}
\end{figure}

~ \\

\twocolumngrid

In some systems, the variance is low but free energy estimates are dependent on the starting poses.
When results from all simulations between milestones c\&d are grouped together, the standard deviation in many systems no longer converges to within 5 $k_B T$.
Because converged thermodynamic quantities should not depend on initial conditions, this inconsistency implies that some simulations are not fully converged.

When free energy estimates contradict, it is not completely clear which result is most consistent with the force field.
False convergence is quantified based on the minimum value, using $\sigma[x,x_{min}]$. A system is considered to be falsely converged when $\sigma[f]$ is less than 5 $k_B T$ but $\sigma[f,f_{min}]$ is greater than 15 $k_B T$.
Based on this definition, false convergence occurs in a nontrivial fraction of systems (Table \ref{tab:nconverged}).

Ultimately, poor convergence and false convergence are caused by incomplete sampling.
The most straightforward explanation for incomplete sampling is that a simulation starts and remains trapped in a local minimum distinct from the global minimum energy.
Indeed, between milestones c\&d, most falsely converged calculations starting from the crystal structure (4/5) and docked poses (7/11) have no starting pose within 1.0 \AA~ of the lowest-energy observed pose.
However, sampling the lowest-energy structure is not a sufficient condition for complete sampling,
as other poses may have a lower \emph{free energy}.
For simulations starting from the lowest-energy observed pose, four systems appear falsely converged between milestones c\&d (formatted by PDB identifier ($\sigma[f], \sigma[f,f_{min}]$) with $\sigma$ in units of $k_B T$):
1hww (0.353, 16.1), 1r1h (3.95, 15.3), 1hp0 (0.496, 16.7), 1gm8 (0.422, 15.6).
For 1hww (Figure \ref{fig:dock_last}(c) and (d)) and 1gm8 (Figure \ref{fig:dock_last}(g) and (h)), starting from the crystallographic pose leads to a lower estimate of $f_{cd}$ than the lowest-energy pose.
For 1hp0 (Figure \ref{fig:dock_last}(e) and (f)), a docked pose has the lowest estimated $f_{cd}$.
1r1h does not appear to sample different poses;
the large $\sigma[f,f_{min}]$ is the result of an outlier with a significantly lower $f_{cd}$ estimate than average.

\section{Discussion}

The selection of intermediate states between two thermodynamic milestones of interest is a ubiquitous problem in molecular simulation. The usual approach is an iterative trial-and-error process starting from a naive protocol, checking for issues such as replica exchange bottlenecks, and inserting and removing states as necessary. I have developed a simple and robust approach to initialize a series of thermodynamic states based on only a single adjustable parameter, the thermodynamic speed. Due to this automated procedure, I was able to run a large number of simulations on a diverse array of protein-ligand complexes. In the vast majority of cases, the protocols yielded consistent replica exchange rates across neighboring thermodynamic states without further fine-tuning. The described approach to trailblazing thermodynamic state space may find use in other classes of simulations.

Replica exchange calculations in this present study include more thermodynamic states than most published molecular simulations. Conventional wisdom about replica exchange is that an optimal number of replicas will maximize efficiency. With too few replicas, there is limited phase space overlap between neighbors and exchange rates are vanishingly small. With too many replicas, metrics of replica exchange efficiency, such as a mean round-trip time, diminish. However, in a recent study involving extensive simulation of several distinct processes, it was found that if there are no bottlenecks, the number of states has little impact on the convergence of free energy estimates (the manuscript is under preparation). Hence, I chose to include a large number of states to allow for the possibility that later sampling will explore different regions of phase space and reduce the exchange rate between neighboring states.

While useful, the described thermodynamic state initialization process remains imperfect. Replica exchange is particularly beneficial when the important phase space of a thermodynamic states is a subset of the important phase space of another. The present procedure is limited to the variation of a single thermodynamic parameter between milestones, and can do little to promote the subset relationship. Future improvements could accomodate varying multiple parameters (e.g. separate parameters for the temperture, van der Waals grids, and electrostatic grids) in between thermodynamic states of interest, while maintaining a phase space subset relationship.

While a direct comparison is difficult, the described procedure appears less successful than the best docking programs at identifying the native pose. \citet{Hartshorn2007} found that with a larger site (10 \AA), the top-ranked GOLD solution is within 2 \AA ~ of the experimental binding mode in 80.4\% of calculations. In contrast, only 56.5\% of the present UCSF DOCK 6 calculations were successful, suggesting that setup and parameterization, and perhaps the docking algorithm, could be fine-tuned. Native pose identification based on AlGDock performs similarly to UCSF DOCK 6, suggesting that replica exchange sampling from the Boltzmann distribution does not significantly increase the diversity of poses. As evidenced by the variable performance as a function of starting pose(s), improving Boltzmann sampling schemes should improve native pose identification.

Nonetheless, the sampling and estimation approach is able to attain precise BPMF estimates in the majority of systems. As for the others, it is unsurprising that the same approach does not work equally well in all cases. With a diverse set, ligands and receptors will exhibit distinct levels of complexity and may be driven to bind by different factors. Indeed, the inadequate sampling that led to convergence issues for several systems is probably not unique to BPMF calculations \cite{Mobley2012b}. Because of the relative simplicity of the calculations, I have been able to conduct more replicates from more distinct starting points than other alchemical standard binding free energy calculations. The sheer number of simulations has allowed us to expose problems that may not be clear from fewer calculations. Our results suggest that Hamiltonian replica exchange with molecular dynamics is not always an adequate approach for sampling distinct ligand binding poses. Specialized sampling methods will likely be required in a more broadly effective protein-ligand BPMF estimation procedure. It is also clear that for some systems, there may be limited phase space overlap between different implicit solvent models. More precise BPMF estimation can be attained by introducing intermediate states or implementing more similar force fields for milestones a\&b and d\&e.

\section{Conclusions}

I have developed a method to initialize thermodynamic states and estimate BPMFs for protein-ligand systems. The largest sources of imprecision are found to be ligand pose sampling and phase space overlap between implicit solvent models. Intriguingly, I found that starting from the lowest-energy pose is not a sufficient condition for converged BPMF estimates.

Future studies may build on the current work to compare BPMFs with molecular docking scores as classifiers of binding. Another useful direction will be to assess the convergence of standard binding free energy estimates upon increased receptor conformational sampling.

\section{Acknowledgments}

I thank Michael Shirts (University of Virginia), John Chodera (MSKCC), and Trung Hai Nguyen (IIT) for helpful discussions.
Peter Eastman (Simbios) and John Chodera assisted with implementing an early version of the code in OpenMM.
Michael Sherman (Simbios) made an invaluable suggestion of using precomputed grids.
Rob Gardner, Lincoln Bryant, and Balamurugan Desinghu (Open Science Grid) and
Tom Milledge (DSCR) provided assistance with computing resources.
David Beratan (Duke) was a supportive postdoctoral advisor.
OpenEye Scientific Software, Inc. and UCSF provided academic licenses to their software.
At the beginning of this project, I was a postdoctoral scholar at Duke, supported by NIH 2P50 GM-067082-06-10.
I also spent a month as an OpenMM visiting scholar at Simbios.

\section{References}

redpajama_set_name: RedPajamaArXiv


I offer complete marketing coordination, from creating a strategic marketing plan and budget, to hands-on execution of all tactical work and vendor negotiations, and finally, to measurement of campaign effectiveness. My extreme attention to detail and organizational skills are well suited to small teams that need to accomplish promotional goals with limited resources. When efficiency is paramount and teams are challenged to do a lot with a little, I provide an all-in-one marketing solution.


Impact Mill Size And Shape Of Gold Grains - … Impact Mill Size And Shape Of Gold Grains. We have sold many of crushing and grinding units to more than 130 countries!
SAG mills are primarily used at gold, impact mill size and shape of gold grains. Impact mill size and shape of gold grains Manganese Crusher... Learn More.


text:
Daily Ratings and News for Exxon Mobil
Complete the form below to receive the latest headlines and analysts' recommendations for Exxon Mobil with our free daily email newsletter:
Short Interest in Acushnet Holdings Corp (NYSE:GOLF) Drops By 36.2%
Legal & General Group (LON:LGEN) Stock Rating Reaffirmed by Credit Suisse Group
Weekly Analysts' Ratings Updates for Aravive (ARAV)
Repay (NASDAQ:RPAY) Price Target Raised to $20.00
Short Interest in First Industrial Realty Trust, Inc. (NYSE:FR) Declines By 26.6%
HC2 Holdings Inc (NYSE:HCHC) Short Interest Down 11.0% in December
Fortis Inc (NYSE:FTS) Short Interest Down 49.3% in December
Huami Corp (NYSE:HMI) Short Interest Down 19.0% in December
Short Interest in FGL Holdings (NYSE:FG) Decreases By 36.7%
Insider Buying: Osisko Metals Inc (CVE:OM) Director Buys 50,000 Shares of Stock
J.Jill (NYSE:JILL) Releases Q4 2019 IntraDay Earnings Guidance
Foster Victor Wealth Advisors LLC Purchases 4,283 Shares of Verizon Communications Inc. (NYSE:VZ)
Philip Morris International (NYSE:PM) Downgraded by ValuEngine
Canacol Energy Ltd (TSE:CNE) Director Sells C$67,425.00 in Stock
Tilly's (NYSE:TLYS) Releases Q4 Earnings Guidance
Best Buy Co Inc (NYSE:BBY) Sees Large Drop in Short Interest
Brookfield Business Partners LP (NYSE:BBU) Short Interest Update
Brown-Forman Co. (NYSE:BF.A) Short Interest Update
Short Interest in Brandywine Realty Trust (NYSE:BDN) Increases By 27.1%
Biglari Holdings Inc (NYSE:BH.A) Short Interest Update
Calamos Wealth Management LLC Sells 2,045 Shares of Exxon Mobil Co. (NYSE:XOM)
Posted by Amber Abramson on Jan 16th, 2020
Calamos Wealth Management LLC lessened its stake in Exxon Mobil Co. (NYSE:XOM) by 4.4% in the 4th quarter, according to the company in its most recent Form 13F filing with the SEC. The firm owned 44,538 shares of the oil and gas company's stock after selling 2,045 shares during the period. Calamos Wealth Management LLC's holdings in Exxon Mobil were worth $3,108,000 at the end of the most recent reporting period.
Several other hedge funds and other institutional investors also recently added to or reduced their stakes in XOM. BlackRock Inc. lifted its stake in Exxon Mobil by 0.4% during the second quarter. BlackRock Inc. now owns 280,882,107 shares of the oil and gas company's stock worth $21,523,998,000 after purchasing an additional 1,069,038 shares during the period. Nuveen Asset Management LLC lifted its stake in Exxon Mobil by 5,514.1% during the second quarter. Nuveen Asset Management LLC now owns 22,672,032 shares of the oil and gas company's stock worth $1,737,358,000 after purchasing an additional 22,268,191 shares during the period. California Public Employees Retirement System lifted its stake in Exxon Mobil by 22.2% during the third quarter. California Public Employees Retirement System now owns 10,458,189 shares of the oil and gas company's stock worth $738,453,000 after purchasing an additional 1,902,091 shares during the period. Eaton Vance Management lifted its stake in Exxon Mobil by 0.4% during the second quarter. Eaton Vance Management now owns 8,034,970 shares of the oil and gas company's stock worth $615,684,000 after purchasing an additional 28,217 shares during the period. Finally, Rhumbline Advisers lifted its stake in Exxon Mobil by 1.2% during the third quarter. Rhumbline Advisers now owns 7,678,548 shares of the oil and gas company's stock worth $542,182,000 after purchasing an additional 93,918 shares during the period. Institutional investors and hedge funds own 54.51% of the company's stock.
Get Exxon Mobil alerts:
NYSE XOM opened at $69.09 on Thursday. The firm has a market cap of $292.54 billion, a price-to-earnings ratio of 14.01, a PEG ratio of 2.19 and a beta of 1.00. The company has a 50 day moving average of $69.50 and a 200-day moving average of $70.74. The company has a current ratio of 0.78, a quick ratio of 0.50 and a debt-to-equity ratio of 0.13. Exxon Mobil Co. has a 1-year low of $66.31 and a 1-year high of $83.49.
Exxon Mobil (NYSE:XOM) last posted its quarterly earnings data on Friday, November 1st. The oil and gas company reported $0.68 earnings per share for the quarter, topping the Thomson Reuters' consensus estimate of $0.64 by $0.04. The company had revenue of $65.05 billion during the quarter, compared to analyst estimates of $65.86 billion. Exxon Mobil had a net margin of 5.43% and a return on equity of 7.45%. The company's revenue was down 15.1% on a year-over-year basis. During the same quarter last year, the business posted $0.82 earnings per share. On average, analysts forecast that Exxon Mobil Co. will post 2.42 earnings per share for the current year.
A number of research analysts have commented on XOM shares. Goldman Sachs Group cut their price target on shares of Exxon Mobil from $74.00 to $72.00 and set a "neutral" rating for the company in a research note on Friday, January 3rd. Citigroup dropped their price objective on shares of Exxon Mobil from $75.00 to $68.00 and set a "neutral" rating for the company in a report on Tuesday. Bank of America raised shares of Exxon Mobil to a "buy" rating in a report on Tuesday, January 7th. Piper Jaffray Companies initiated coverage on shares of Exxon Mobil in a report on Tuesday, November 26th. They set a "hold" rating and a $83.00 price objective for the company. Finally, Berenberg Bank reaffirmed a "neutral" rating and set a $67.00 price objective on shares of Exxon Mobil in a report on Tuesday. One analyst has rated the stock with a sell rating, nineteen have assigned a hold rating and three have issued a buy rating to the company. The company currently has an average rating of "Hold" and a consensus price target of $78.05.
In related news, VP David S. Rosenthal sold 9,615 shares of the firm's stock in a transaction that occurred on Wednesday, December 4th. The shares were sold at an average price of $69.51, for a total value of $668,338.65. Also, VP Liam M. Mallon sold 32,287 shares of the firm's stock in a transaction that occurred on Wednesday, November 20th. The shares were sold at an average price of $68.26, for a total transaction of $2,203,910.62. Following the completion of the sale, the vice president now directly owns 211,289 shares of the company's stock, valued at approximately $14,422,587.14. The disclosure for this sale can be found here. Company insiders own 0.05% of the company's stock.
Exxon Mobil Profile
Exxon Mobil Corporation explores for and produces crude oil and natural gas in the United States, Canada/Other Americas, Europe, Africa, Asia, and Australia/Oceania. It operates through Upstream, Downstream, and Chemical segments. The company is also involved in the manufacture, trade, transport, and sale of crude oil, petroleum products, and other specialty products; and manufactures and markets petrochemicals, including olefins, polyolefins, aromatics, and various other petrochemicals.
Featured Story: What sectors are represented in the Nikkei Index?
Want to see what other hedge funds are holding XOM? Visit HoldingsChannel.com to get the latest 13F filings and insider trades for Exxon Mobil Co. (NYSE:XOM).
Receive News & Ratings for Exxon Mobil Daily - Enter your email address below to receive a concise daily summary of the latest news and analysts' ratings for Exxon Mobil and related companies with MarketBeat.com's FREE daily email newsletter.

redpajama_set_name: RedPajamaCommonCrawl


Paper is the regular reward for very first year anniversary and even though it can just mean a easy card, there are a good deal of one year anniversary reward concepts that are paper but a little bit kiss japan tour Setlist much more enjoyable and imaginative.
Along with all of the sightseeing and wonderful food, you'll have the prospect to fulfill some interesting individuals along the way. Pros in the Town is a great corporation that presents chances to meet up with new individuals. You'll go household with so a lot of new pals you will want to appear again up coming calendar year to celebrate with them all over once again.
Toby's most recent album, Bullets in the Gun has currently produced three prime twenty hits, which includes the title observe. His most up-to-date monitor, Designed in The us, is currently at #22 on the Billboard Nation Chart. Toby fans will no doubt walk away from this live performance with a few new preferred songs u kiss concert 2016 to insert to the huge collection of hits that Toby has presented them above the a long time.
Okay enable's say that you and this identical guy have various events going 1 night-allow's say he has a sporting celebration he desires to go to and you experienced kiss concert new zealand kiss tickets burgettstown to your beloved band. A sacrifice would be offering up those kiss tickets sydney to go observe the sport with him, being aware of it is crucial to him, even even though you had been absolutely psyched to see the kiss concert 10/29/14.
Enjoy Live Music furnished by neighborhood New Zealand functions, a picnic styled lunch served at the Wynyard backyard with gourmand foodstuff and wines created domestically. There's every little thing you want to take pleasure in a working day in the sun at the Culverden Christmas Country Fete.
On Sunday, May possibly nine, 2010 at 2:00 p.m., http://meftahelkheir.com/Index.php/profile/77 the Alabama Shakespeare Festival will present Hamlet by William Shakespeare. Ticket selling prices are $42 for Adult Orchestra seating, $32 for Adult Balcony seating, $35 for Associates, and $28 for these beneath 25 a long time aged.
Give them some thing that could enable them modify to their lifestyle immediately after retirement. Permit your favourite retiree know that there is certainly existence following retirement. And most likely a single of the finest ways to do it is by offering them a properly viewed as guide concerning the matter. There are numerous excellent publications that permit the reader enjoy how he or she can dwell lifestyle to the fullest, even following indicating goodbye to the workforce. Offering books that could assist them learn the potentials of other worthwhile cash flow-making functions is also an great idea.


text:
Vaux's Happening Home
Monroe: Where it all began
Washington State boasts only two well-known Vaux's Swift migratory communal roost sites, and barely known, is probably more accurate. One of the first American structures that greeted Canadians entering at Sumas was the Old Customs House. Large numbers of swifts had been seen entering its brick chimney as early as the 1950s, especially in the spring. In 1991 a successful effort was made to save this building from the wrecking ball. Success was assured by the National Historic Preservation Act, but the birds were also a factor. Enjoy the riveting story as told by the Sumas Astonisher.
The second site belongs to the Monroe School District. In the middle of town sits a snazzy 1930s brick schoolhouse. Its name honors Frank Wagner, whoever he was. It also has a 4x4 foot slightly tapered brick chimney that extends 25 feet through the school's flat roof. It's not the largest chimney around, but for the swifts, it's a comfy five star accommodation. One local citizen claims knowing about the swifts for 40 years, but the birding community wasn't aware on this avian spectacle until the fall of 2004, when a few posts surfaced on Tweeters. As of July 2008, the majority of Monroe town folk have no idea that they have a twice-yearly bird invasion. Audubon is working hard to change that. The problem is that the building is old and functional, but the chimney isn't, and a brick storm at recess could be a fatal event. Better tear that thing down.
Audubon quickly launched operation S.O.S. = SAVE OUR SWIFTS. The Monroe School District's response to these mostly outside bird lovers was, well, actually, extremely positive. In March the Everett Herald came out with an article about the unfolding drama.
But if the Wagner chimney is torn down, surely the birds will just go down the street to some other chimney. Won't they?
Vaux's Happening would document just how significant and important Monroe was to migrating swifts statewide. Long before the last bird was counted, Monroe's significance was obvious. Of the 183,463 documented Washington State Spring 08 Vaux's Swift chimney encounters, a stunning 74% occurred in Monroe. The May 7 high count of over 21,000 had to be a standing room only evening for the little birds, and probably doomed the shorten-the-chimney option.
What is it about Monroe's Wagner Elementary School Chimney that makes it so attractive? A nice clear approach to the chimney entrance may be a part of it, but as every real estate agent is fond of saying, it's location, location, and location.
Frank Wagner Elementary School
Monroe Swifts website
Inside the Chimney
The security cameras installed in late August 2010 in the Frank Wagner Elementary chimney in Monroe captured these clips on September 4, 2010.


This inquiry-based unit presents students with a range of visual primary sources to spark curiosity about life in the 1800s. Each activity introduces a new concept related to the Australian Gold Rushes.
In this unit of work students are asked to consider how the National Museum of Australia's First Australians gallery presents Indigenous cultures and history. The unit includes a number of activities that help students explore the stories and experiences of Australia's Indigenous peoples.
This resource is a webpage with information, study guide and resources on the process of analysing and evaluating historical sources to support Stage 3, 4 and 5 HSIE and the Australian Curriculum: History.
This resource is about exploring Governor Macquarie's preferred residence while learning more about the early colony through this virtual tour of Old Government House at Parramatta.
This resource is a YouTube playlist containing a series of videos taken as a group of senior high school students are given a guided tour of NSW Government House in 2010. The tour covers primary sources such as architecture, furniture and images significant to the history of Australia and NSW.


Paul Simon's career began with Art Garunkel as the folk duo Tom and Jerry in Forest Hills, New York. Their partnership flourished in the mid-60's with numerous hits. Simon went solo in the early 70's, and continued to top the charts on his own, including the ground- breaking 1986 album "Graceland."


text:
Q: Algorithm for Client-Server Games For stand alone games, the basic game loop is (source: wikipedia)
while( user doesn't exit )
  check for user input
  run AI
  move enemies
  resolve collisions
  draw graphics
  play sounds
end while

But what if I develop client-server-like games, like Quake, Ragnarock, Trackmania, etc,
What the Loop/Algorithm for client and the server parts of the game?

A: It really isn't a simple problem. At a most basic level you could say that the network provides the same data that the MoveEnemies part of the original loop did. So you could simply replace your loop with:
while( user doesn't exit )
  check for user input
  run AI
  send location to server
  get locations from server
  resolve collisions
  draw graphics
  play sounds
end while

However you need to take into account latency so you don't really want to pause your main loop with calls to the network. To overcome this it is not unusual to see the networking engine sitting on a second thread, polling for data from the server as quickly as it can and placing the new locations of objects into a shared memory space:
while(connectedToNetwork)
    Read player location
    Post player location to server
    Read enemy locations from server
    Post enemy locations into shared memory

Then your main loop would look like:
while( user doesn't exit )
  check for user input
  run AI
  read/write shared memory
  resolve collisions
  draw graphics
  play sounds
end while

The advantage of this method is that your game loop will run as fast as it can, but the information from the server will only be updated when a full post to and from the server has been completed. Of course, you now have issues with sharing objects across threads and the fun with locks etc that comes with it.
On the server side the loop is much the same, there is one connection per player (quite often each player is also on a separate thread so that the latency of one won't affect the others) for each connection it will run a loop like
while (PlayerConnected)
    Wait for player to post location
    Place new location in shared memory

When the client machine requests the locations of enemies the server reads all the other players locations from the shared block of memory and sends it back.
This is a hugely simplified overview and there are many more tweaks that will improve performance (for instance it may be worth the server sending the enemy positions to the client rather than the client requesting them) and you need to decide where certain logically decisions are made (does the client decide whether he has been shot because he has the most up to date position for himself, or the server to stop cheating)

A: The client part is basically the same, except replace
run AI
move enemies
resolve collisions

with
upload client data to server
download server updates

And the server just does:
while (game is running)
{
    get all clients data
    run AI
    resolve collisions
    udpate all clients
}

A: You can use almost the same thing, but the most of you logic would be on the server, you can put timers, sounds, grafics, and other UI components on the client app.
Any business rule (AI,Movements) goes in the server side.

A: It would be something like
Client:
while( user does not exit )
    check for user input
    send commands to the server
    receive updates about the game from the server
    draw graphics
    play sounds
end

Server:
while( true )
    check for client commands
    run AI
    move all entities
    resolve collisions
    send updates about the game to the clients
end

A: A very useful and I would argue pertinent paper to read is this one: Client-Server Architectures
I gave it a read and learned a lot from it, a lot of sense was made.  By separating out your game into strategically defined components or layers, you can create a more maintainable architecture.  The program is easier to code, and more robust than a conventional linear program model like the one you've described.
That thought process came out in a previous post here about using a "Shared Memory" to talk between different parts of the program, and so overcoming the limitations of having a single thread and step-followed-step game logic.
You can spend months working on the perfect architecture and program flow, read a single paper and realise you've been barking up the wrong tree.
tldr; read it.

A: Client:
connect to server
while( user does not exit && connection live)
    check for user input
    send commands to the server
    estimate outcome and update world data with 'best guess'
    draw graphics
    play sounds
    receive updates about the game from the server
    correct any errors in world data
    draw graphics
    play sounds
end

Server:
while( true )
    check for and handle new player connections
    check for client commands
    sanity check client commands
    run AI
    move all entities
    resolve collisions
    sanity check world data
    send updates about the game to the clients
    handle client disconnects
end

The sanity checks on the client commands and world data are to remove any 'impossible' situations caused either by deliberate cheating (moving too fast, through walls etc) or lag (going through a door that the client thinks is open, but the server knows is closed, etc).
In order to handle lag between the client and server the client has to make a best guess about what will happen next (using it's current world data and the client commands) - the client will then have to handle any discrepancies between what it predicted would happen and what the server later tells it actually happened. Normally this will be close enough that the player doesn't notice the difference - but if lag is significant, or the client and server are out of synch (for example due to cheating), then the client will need to make an abrupt correction when it receives data back from the server.
There are also lots of issues regarding splitting sections of these processes out into separate threads to optimise response times.
One of the best ways to start is to grab an SDK from one of the games that has an active modding community - delving into how that works will provide a good overview of how it should be done.

redpajama_set_name: RedPajamaStackExchange


text:
Humanscale to Exhibit Designs that Create an Ergonomic Workstation at Bureaux Expo
Ergonomic workstations promote comfort, health and productivity while minimising the risk of pain and reducing the risk of musculoskeletal disorders
PARIS — Ergonomics leader Humanscale will exhibit products that reflect their vision of an ergonomic workstation at Bureaux Expo from April 16th to 18th in Paris. Bureaux Expo will showcase 200 exhibitors and host more than 10,000 visitors this year. Visitors to the Humanscale booth will get to experience the Float™ height-adjustable table, the Element Disc™ LED task light and the M/Flex™ monitor support system. Each of these office tools represents a key component of Humanscale's ergonomic workstation. When used together, these products increase comfort, health and productivity among office workers.
Humanscale has created the Float height-adjustable table to promote movement and a healthy lifestyle. With recent research suggesting that sitting all day is harmful to one's health, height-adjustable work surfaces such as Float are more essential than ever, as they encourage a combination of sitting and standing. The Float height-adjustable table meets the demands of healthier work environments with an unprecedented combination of adjustability for a vast range of heights (51 cm of height adjustability) and weights (up to 59 kg of desktop weight, including the tabletop), along with streamlined good looks. With a unique counterbalance mechanism that provides constant counter-weighting force, it enables simple, one-hand weight adjustment, and fluid height adjustment for ergonomic sit-stand usage—all without the use of an obstructive crossbeam or a cumbersome crank. The Float height-adjustable table is also built with measures to prevent the possibility of injury―its innovative safety feature only allows adjustment when the desktop and contents are properly balanced.
Additionally, the Float height-adjustable table accommodates a range of tabletop widths from 120 cm to 180 cm, and can be used with any manufacturer's tabletops to ensure aesthetic consistency with any installed or specified systems furniture. Made for a long life, it is composed predominantly of recyclable steel and is designed with minimal parts and manufacturing processes. Unlike electric-powered height-adjustable tables, it is able to balance the table's load without using energy, which makes it exceptionally sustainable.
In addition to showing a height-adjustable table, Humanscale will exhibit a task light that offers customisable illumination and fingertip adjustability to accommodate all of the user's needs. Providing complete user control, the Element Disc LED task light features intuitive, smart dimming functionality, which affords the user seven levels of brightness and uses memory to save the last preferred light setting selected by the user. Creating an ultra-wide footprint of warm, glare-free light, it uses the innovative Thin Film LED Technology that was first manifested in the award-winning Horizon™ LED task light. For exceptional sustainability, the Element Disc LED task light is equipped with a passive infrared (PIR) occupancy sensor, which automatically turns the light off or on when the user leaves or enters the area.
In addition to these components of an ergonomic workstation, Humanscale will show the M/Flex monitor support system, which supports six monitors from a single post and up to 12 monitors on two posts united via crossbars. Monitor arms such as M/Flex bring the work to the user, clear valuable desktop space and encourage easy positioning to support the user's every task. With its innovative design, users can choose a scalable platform to accommodate many monitor display applications (and their various specific hardware, desk configurations or multiple monitor requirements), including notebooks, touchscreens and monitors in a wide range of weights. All M/Flex monitor support systems use the award-winning M2™ and M8™ monitor arm components, which can accommodate many configurations to meet the user's exact requirements. The scalable M/Flex monitor support system is an eco-conscious, long-lasting solution. With fewer parts, it is tremendously durable, since there are fewer components that can wear out or break over time. By utilising fewer parts, the M/Flex monitor support system is also lighter and smaller, which means that less material and energy are used in the manufacturing and shipping of the product.
These and other innovations can be experienced in person at booth K29. For more information about Humanscale's ergonomic solutions, visit www.humanscale.com or call 800-400-0625.