sing Selectively Applied Accelerated Molecular Dynamics J. Chem. Theory Comput., Vol. 6, No. 11, 2010 3287 nF exp( ∆ V(r)i) ∑ i)0 1 + exp[- (M + Wi - ∆ F)] ) nR exp( ∆ V(r)j) ∑ j)0 1 + exp[- (M + Wj - ∆ F)] (6) is the optimal free energy estimate between adjacent windows. It has recently been shown that reweighting of states in BAR to account for non-Boltzmann sampling may have practical advantages outside of aMD simulations.34 3. Computational Details Molecular Dynamics Details. MD simulations were performed with the MD package Desmond (version 2.2) developed by D. E. Shaw Research.35 Both systems were built, solvated, and ionized with Schro¨dinger’s Maestro modeling suite such that there was a minimum of 12 Å of TIP3P36 water buffer between the macromolecule and the periodic boundary and an ionic concentration of ∼150 mM NaCl was present. The CHARMM22 force ﬁeld with the CMAP correction was utilized (except where noted below in the neuraminidase calculations).37 Following 10 000 steps of minimization, systems were continuously heated to 300 K over 1.5 ps. All simulations used the Martyna-Tobias- Klein constant pressure and temperature algorithms (a combination of Nose`-Hoover constant temperature and piston constant pressure algorithms)38,39 with a reference temperature and pressure of 300 K and 1.01425 bar, respectively. Short-range nonbonded interactions were trun- cated at 12 Å, while long-range electrostatics were calculated with a particle-mesh Ewald algorithm using a sixth-order B-spline for interpolation and a grid spacing of <1 Å in each dimension.40 A time step of 2 fs was employed, and the M-SHAKE algorithm was used for constraining all hydrogen- containing bonds.41 A plugin was written for Desmond to perform aMD calculations on speciﬁed dihedrals. Alanine Dipeptide. An alanine dipeptide molecule based on a model compound was solvated in a (27 Å)3 box using Maestro and equilibrated for 5 ns following heating. Two sets of 50 ns aMD trajectories were run, one in which all dihedrals were accelerated with aMD and one in which only the two dihedrals deﬁned as φ and ψ were accelerated (selective aMD, see Figure 1). For each setup, 16 simulations spanning the parameter space of E and R were run to optimize these parameters. For the all-dihedral simulations, parameters of E ) 8 and R ) 4 had an optimal ﬁt to our metric (as discussed below), whereas for the selective aMD E ) 1 and R ) 0.75 were optimal. Additionally, a 250 ns classical MD simulation was performed. To determine the accuracy of the aMD results, well- tempered metadynamics was performed to calculate the underlying two-dimensional free energy landscape in φ/ψ space.42 In well-tempered metadynamics, the height of a Gaussian centered at position x is proportional to the Boltzmann weight of the metadynamics potential already present at x, that is the added Gaussian has a maximum value of ω0 · exp (- Vt(x)/k∆ T), with ω0 being the initial Gaussian height, -Vt(x) the metadynamics potential at x, and k∆ T a Figure 1. Conformation of the model system alanine dipep- tide can be expressed by the angles of the two torsions φ and ψ. user-deﬁned energy which limits the explored energy range. We performed a 50 ns simulation in which gaussians were added every 0.2 ps with a width of 0.1 radians and a height determined by ω0 ) 0.02 kcal/mol and k∆ T ) 2.4 kcal/ mol. To quantitate convergence we deﬁne a metric as follows ) 1 N 3.0 ∑ ∑ N · M i)1 U).6,.7,... {PPiaie,,MxUUaDct//PPiei,,xUUact if Pie,xUact > Pi,U aMD aMD (7) aMD if Pie,xUact < Pi,U For a given well (deﬁned as i in Figure 2f) this metric calculates the ratio of the population of states below the energy U (Pai,MUD) with that expected from the metadynamics results (Pie,xUact) and averages this over M energy values of U from 0.6 to 3 kcal/mol (in increments of 0.1 kcal/mol) for each of the N ) 3 wells. If a well has a population greater than that expected from the metadynamics result, the inverse of the ratio is taken to equally account for over and under sampling of the well. This metric has several advantageous features: by averaging over multiple energy levels it selects for smooth population densities (as are observed in cMD and metadynamics but may not result when the trajectory is reweighted with aMD), it treats over and under sampling a well as equally poor, and it equally weights all three of the main energy wells. ) 1 is considered an ideal reproduction of the exact population. N1 Neuraminidase. A monomer of neuraminidase bound to the inhibitor oseltamivir taken from the 2HU0 crystal structure was solvated in an approximately (70 Å)3 box. Following the minimization and heating protocol outlined above, the system was equilibrated for 5 ns with the AMBER99SB force ﬁeld43 before protein and water param- eters were changed to the CHARMM22 force ﬁeld and further equilibrated for 1 ns. Parameters for oseltamivir were previously developed for use with the AMBER99SB force ﬁeld44 and maintained throughout the CHARMM simula- tions. The switch to the CHARMM force ﬁeld was performed after testing of aMD free energy calculations revealed that increased acceleration levels in the AMBER force ﬁeld tended to disturb the electrostatic components of the free energy calculation; therefore, this hybrid AMBER (for the ligand) and CHARMM (for the remainder of the system) force ﬁeld was utilized. While the authors concede this may produce incorrect absolute binding energies, the goal of this paper is to study convergence of accelerated free energy calculations to results obtained from unaccelerated (and longer) calculations. Alchemical free energy calculations for the decoupling of the ligand in the protein’s active site were performed using 21 windows in which the electrostatics were decoupled over