8488
J. Phys. Chem. B 2006, 110, 8488-8498
Sampling of Rare Events Using Hidden Restraints Markus Christen, Anna-Pitschna E. Kunz, and Wilfred F. van Gunsteren* Laboratory of Physical Chemistry, Swiss Federal Institute of Technology Zu¨rich, ETH, CH-8093 Zu¨rich, Switzerland ReceiVed: January 24, 2006; In Final Form: March 2, 2006
A method to enhance sampling of rare events is presented. It makes use of distance or dihedral-angle restraints to overcome an energy barrier separating two metastable states or to stabilize a transition state between the two metastable states. In order not to perturb these metastable end states themselves, a prefactor is introduced into the restraining energy function, which smoothly increases the weight of this function from zero to one at the transition state or on top of the separating energy barrier and then decreases the weight again to zero at the final state. The method is combined with multi-configurational thermodynamic integration and applied to two biomolecular systems, which were difficult to treat using standard thermodynamic integration. As first example the free energy difference of a cyclic R-aminoxy-hexapeptide-ion complex upon changing the ion from Cl- to Na+ was calculated. A large conformational rearrangement of the peptide was necessary to accommodate this change. Stabilizing the transition state by (hidden) restraints facilitates that. As a second example, the free energy difference between the 4C1 and the 1C4 conformation of β-D-glucopyranoside was calculated. In unrestrained simulations the change from the 4C1 into the 1C4 conformation was never observed because of the high energy barrier separating the two states. Using (hidden) restraints, the transition from the 4 C1 into the 1C4 state and back could be enforced without perturbing the end states. As comparison, for the same transitions the potential of mean force as obtained by using dihedral-angle constraints is provided.
1. Introduction Computer simulation is increasingly used to investigate, in atomic detail, dynamic molecular processes. Atomistic molecular dynamics (MD) simulations of systems comprising tens to even hundreds of thousands of atoms are now possible covering nanoseconds to microseconds. Yet, conformational changes in biomolecules, ligand binding or structural organization occur roughly in the millisecond range, while protein folding may take up to minutes. This means that these interesting events in biomolecular systems are often not or only rarely observed in MD simulations. Therefore, methods to enhance the sampling of these rare events during an MD simulation have been developed over the years.1-3 Often, one is only interested in a few selected degrees of freedom of the system under the average influence of the many other degrees of freedom. Then, statistical mechanics provides a means of expressing the average dependence of the degree of freedom of interest on the residual degrees of freedom in the form of a potential of mean force (PMF).4 The potential of mean force can be defined as minus kBT times the logarithm of the probability that the system is found at a specific position along the degree of freedom under investigation (at temperature T and with kB being the Boltzmann constant), and it can be interpreted as a projection of the free energy on one (or more) coordinates or degrees of freedom of interest. Using MD simulation, the potential of mean force may be calculated by constraining the system to a specified point along the investigated degree of freedom and calculating the average value of the constraint force magnitude.5-7 Another way to calculate a potential of mean force is by using umbrella sampling.8 Sampling is biased toward otherwise unfavorable * Corresponding author e-mail:
[email protected].
regions of the configurational space by performing simulations with an additional artificial biasing potential energy term, the so-called umbrella potential. The resulting probability distribution from the simulation of the unphysical system can be corrected afterward to yield the probability distribution of the corresponding physical, unbiased, system. The biasing potential may be used to focus sampling in a specific region, a so-called window, along the degree of freedom for which the potential of mean force is to be determined. By performing multiple simulations with shifted focuses, the complete degree of freedom can be sampled segment by segment. Each segment of the potential of mean force calculated from a different window has an arbitrary offset. To obtain the resulting continuous potential of mean force, the segments must be matched. A third method to determine a potential of mean force is a combination of the two mentioned above: the biasing force is measured from which the PMF can be constructed.9 All three ways to calculate the potential of mean force, which may then be integrated to get the free energy difference between two states, the so-called end states, on the degree of freedom considered, suffer from the disadvantage that the biasing or restraining (or constraining) potential energy term is also applied in the two end states. This may lead to artifacts in the relative free energies because only the difference between two restrained states, not the difference between two unrestrained states is calculated. In other words, the free energy difference will depend on the degree of freedom or pathway considered, which requires this pathway to be a good approximation of true pathways. In this work we propose a combination of multi-configurational thermodynamic integration with the use of a biasing potential energy term, but without restraining the end states. This makes the end states pathway independent, which allows one to choose an unphysical pathway that enhances the
10.1021/jp0604948 CCC: $33.50 © 2006 American Chemical Society Published on Web 04/06/2006
Sampling of Rare Events Using Hidden Restraints
J. Phys. Chem. B, Vol. 110, No. 16, 2006 8489 The λ derivative of the restraining term V restr is
∂ restr V (r) ) 2n+m (nλn-1(1 - λ)m - mλn(1 ∂λ n m∂ λ)m-1)V AB V AB (r; λ) , (3) restr(r; λ) + λ (1 - λ) ∂λ restr
[
]
and the force is
∂ ∂ (r; λ), fi ) - V restr(r; λ) ) -2n+mλn(1 - λ)m V AB ∂ri ∂ri restr
Figure 1. Weight factor 2n+mλn(1 - λ)m as function of λ for different values of the exponents n and m ) n. The solid line corresponds to n ) 1, the dashed line to n ) 2, the dashed-dotted line to n ) 3 and the dotted line to n ) 4.
sampling. In the next section, the method is explained, followed by an implementation using distance restraints and one using dihedral-angle restraints. Two examples will be given: first, in a cyclic R-aminoxy-hexapeptide-ion complex the ion is changed from an anion into a cation, resulting in a large conformational change of the peptide; second, the free energy change of a transition of β-D-glucopyranoside from the 4C1 to the 1C4 conformation is calculated. The work is concluded by a short discussion. In the Appendix, the formulas to obtain a potential of mean force using distance or dihedral-angle constraints are presented. 2. Method If a molecular system exhibits two (meta) stable states which are connected by a transition path that is only very rarely sampled during an MD simulation,10 it may be efficient to forcefully propagate the system from one end state to the other along this transition path. To be able to calculate an unperturbed relative free energy of the one state with respect to the other it is advantageous not to put any external forces or restraints on the system when being in these two end states. Therefore, the following general restraint formulation is proposed to enforce a transition from state A to B without influencing the end states:
V restr(r) ) 2n+mλn(1 - λ)mV AB restr(r; λ),
(1)
AB where V restr (r; λ) is a λ dependent (restraining) potential energy term which enforces the (smooth) transition of the system from state A into state B along the pathway λ ) 0 to λ ) 1. The weight factor (2n+mλn(1 - λ)m) applied to the restraining AB is shown in Figure 1. potential energy term V restr The relative free energy ∆GBA of state B with respect to state A may be calculated using the (multi-configurational) thermodynamic integration method11 (reviewed elsewhere12,13)
∆GBA )
B ∂H ∫λ)A 〈 ∂λ 〉λdλ,
(2)
with H being the Hamiltonian of the system, which also includes the restraining potential energy term V restr. The ensemble (or time) average of the derivative of the Hamiltonian with respect to λ, (∂H /∂λ), is calculated for a given number Nλ of fixed λ values using Nλ MD simulations, while the integration is carried out numerically afterward.
(4)
where r denotes a configuration of the system and ri the coordinates of particle i. The method may be used with any kind of restraining potential energy function Vrestr(r; λ). The following two sections show distance and dihedral-angle restraints as examples. 2.1. Distance Restraints. A distance restraint k between atoms k1 and k2 may be formulated using a harmonic potential energy function. In practice, it is often necessary to smoothly linearize the potential energy function beyond a specified distance rlin0 to avoid very high energies and extremely steep slopes, V disres (r) ) k
{
0 0 0 0 0 -Kdist k (|rk| - rk + 1/2rlin )rlin |rk| < rk - rlin 0 2 1/2Kdist k (|rk| - rk )
rk0 - rlin0 e |rk| e rk0 + rlin0
0 0 0 Kdist k (|rk| - rk - 1/2rlin )rlin
|rk| > rk0 + rlin0 (5)
where Kdist k is the force constant of the restraint, rk ) rk1 - rk2 is the vector connecting particle k2 to particle k1 and the restraint is linearized for distances larger than rlin0. If the distance restraint should drive the system (in a controlled way) from state A to state B, the restraint parameters have to change along this path. This change of restraint parameters along the path described by λ can be formulated as disres,AB V harm,k (r; λ) ) 0,B 2 1/2((1 - λ)KAk + λKBk )(|rk| - (1 - λ)r0,A k - λrk ) (6)
for the harmonic part and as
(
disres,AB (r; λ) ) ζ((1 - λ)KAk + λKBk ) |rk| - (1 - λ)r0,A V lin,k k -
1 0 0 λr0,B k - ζrlin rlin (7) 2
)
for the linearized part (with ζ ) -1 if |rk| < rk0 - rlin0 and ζ ) 1 if |rk| > rk0 + rlin0) using a restraint length in state A of 0,B A r0,A k and in state B of rk and corresponding force constants Kk B and Kk . The force in the harmonic part is
fi ) -
∂ disres,AB V (r; λ) ∂ri harm,k
) - ((1 - λ)KAk + λKBk )(|rk| - (1 - λ)r0,A k rk λr0,B (δ - δik2) (8) k ) |rk| ik1 and in the linearized part
8490 J. Phys. Chem. B, Vol. 110, No. 16, 2006
fi ) -
Christen et al.
∂ disres,AB V (r; λ) ) -ζ((1 - λ)KAk + ∂ri lin,k rk λKBk )rlin0 (δik1 - δik2), (9) |rk|
where δ is the Kronecker delta. Finally, to use the thermodynamic integration formula (2), the derivative of the distance restraint with respect to integration variable λ needs to be known:
disres,AB V lin,k (φ; λ) )
((1 - λ)KAk + λKBk )(ζ∆φkλ - 1/2φlin0)φlin0 (15) with ζ ) -1 if ∆φkλ < -φlin0 and ζ ) 1 if ∆φkλ > φlin0 for the is kept linearized part of the restraint. The parameter φmax k constant along the path, assuming that common maximum values of the restraining potential energy terms in states A and B can be found. The force in the harmonic part of the restraining potential energy function is
∂∆φkλ ∂φk ∂ dihres,AB Vharm,k fi ) ∂∆φkλ ∂φk ∂ri
∂ disres,AB (r; λ) ) 1/2((KBk - KAk )(|rk| - (1 - λ)r0,A V k ∂λ harm,k A B 2 0,A 0,B λr0,B k ) + 2((1 - λ)Kk + λKk )(|rk| - (1 - λ)rk - λrk ) (r0,A k
-
r0,B k )),
) -((1 - λ)KAk + λKBk )∆φkλ
(10)
(16)
and in the linear part
and
∂∆φkλ ∂φk ∂ dihres,AB Vlin,k fi ) ∂∆φkλ ∂φk ∂ri
∂ disres,AB V (r; λ) ) ∂λ lin,k 1 0,B 0 0 ζ (KBk - KAk ) |rk| - (1 - λ)r0,A k - λrk - ζrlin rlin + 2
(
(
)
((1 -
λ)KAk
+
λKBk )rlin0(r0,A k
-
r0,B k )
{
(φ) Vdihres k
)
0 0 0 -Kdih k (∆φk + 1/2φlin )φlin ∆φk < -φlin 2 1/2Kdih k (∆φk)
-φlin0 e ∆φk e φlin0 (12)
0 0 Kdih k (∆φk - 1/2φlin )φlin
∆φk > φlin0
∂φk ∂ri
(17)
with (∂φk/∂ri) equivalent to the expression used for the (physical) dihedral-angle potential energy term.15-17 Finally, the λ derivative of the restraint is given by
∂ dihres,AB 1 B V ) ((Kk - KAk )(∆φkλ)2 + 2((1 - λ)KAk + ∂λ harm,k 2 0,B λKBk )∆φkλ(φ0,A k - φk )), (18) and
∂ dihres,AB V ) φlin0((KBk - KAk )(ζ∆φkλ - 1/2φlin0) + ∂λ lin,k 0,B ((1 - λ)KAk + λKBk )ζ(φ0,A k - φk )). (19) 3. Applications
with Kdih the restraining force constant, φk0 the angle to k restrain to and ∆φk ) φk - φk0 + 2nπ, where n is chosen such max + 2π] and assuming that that φ is within the range [φmax k , φk φk0 is chosen within the same range. Using this dihedral-angle restraint formulation, φmax determines at which position the k direction of the rotation around the dihedral angle, caused by the restraint potential energy function, inverses. By aligning to the highest potential energy barrier for the rotation φmax k around the dihedral angle, it is possible to avoid pushing against this barrier, by rotating the other way instead. Similar to distance restraining a λ dependence can be introduced to enforce conformational sampling along a pathway from state A to state B disres,AB (φ; λ) ) 1/2((1 - λ)KAk + λKBk )(∆φkλ)2, (13) V harm,k
using 0,B ∆φkλ ) φk - (1 - λ)φ0,A k - λφk + 2nπ,
) -((1 - λ)KAk + λKBk )ζφlin0
). (11)
An example of using distance restraints in a relative free energy calculation is given in Section 3. 2.2. Dihedral-Angle Restraints. Many structural differences of biomolecules are related to different dihedral angles, an example being the ψ and φ angles in peptides.14 Similar to distances, dihedral angles may be restrained using a harmonic potential (for a dihedral k specified by the atoms k1 - k2 - k3 - k 4)
and
∂φk , ∂ri
(14)
In this section, results from studies on two types of problems are shown. First, a cyclic aminoxy-hexapeptide showing the interesting possibility to bind anions and (less tightly) also cations was investigated. When changing from binding an anion to binding a cation, the structure of the peptide has to change dramatically. Standard simulation techniques fail to reproduce this conformational change on a time-scale accessible to MD simulations. Second, the relative stability of chair conformations (1C4 and 4C1) of a hexopyranose was investigated. During a standard MD simulation, not enough transitions between the two conformations occur to accurately calculate their free energy difference. In both systems,the standard simulation does not appropriately sample the part of the phase space important to calculate the free energy differences between different metastable conformations of the systems. But, while in the former system, the complex between ion and peptide becomes instable upon (alchemically) changing from a cation to an anion and, therefore, needs stabilization of the transition state, the latter system needs to be forced over the barrier separating conformation A from B, i.e., needs destabilization of the transition state. Both requirements may be fulfilled by adequately choosing (hidden) restraints.
Sampling of Rare Events Using Hidden Restraints
Figure 2. Configurations of the cyclic aminoxy-hexapeptide complexed with Cl- (A) and with Na+ (B) from MD simulations in (explicit) chloroform.
3.1. Conformations of a Cyclic Aminoxy-Hexapeptide Upon Binding Cations and Anions. Recently, a cyclic D,LR-aminoxy-hexapeptide with the sequence (D)-Leu-(L)-Phe(D)-Leu-(L)-Phe-(D)-Leu-(L)-Phe was synthesized and its properties investigated18 as part of a study of the functions of foldamers, that is, of oligomers of unnatural peptide moieties that fold into well-defined secondary structures.19-22 Cyclic D,L-R amino acid peptides and β3-amino acid peptides have been found to self-assemble into nanotubes and function as transmembrane ion channels,23-28 while R-aminoxy acids, as backbone analogue of β-amino acids, can form an eightmembered intramolecular hydrogen bonded ring,29,30 and, as oligomers 1.88 helices. The cyclic peptide with its small pore size was expected to bind ions. The carbonyl groups may coordinate with cations, whereas the amide hydrogen atoms may form hydrogen bonds with anions (see Figure 2). Experiments found only halide ions, but not alkali metal ions, to bind to the cyclic hexapeptide.18 Here, conformational changes of the hexapeptide upon changing the complexated ion from Cl- to Na+ were investigated, using the multi-configurational thermodynamic integration method. The complex was simulated by integrating Newton’s equation of motion based on the leapfrog scheme31 in explicit chloroform solvent,32 at a temperature of 300 K (maintained by separate weak coupling33 of solute and solvent to a temperature bath with τT ) 0.1 ps), and constant pressure of 1 atm (τP ) 0.5 ps, kI ) 4.575 10-4 (kJmol-1nm-3)-1, using isotropic coordinate scaling.33) Bond lengths were constrained using the SHAKE algorithm34 with a relative geometric tolerance of 10-4. A triple range cutoff-scheme35 was used with a short cutoff of 0.8 nm and a long cutoff of 1.4 nm, and a reaction-field approximation36 (rf ) 4.81) was applied. Center of mass translation was removed. For interaction parameters the GROMOS 45A317,37 parameter set of the GROMOS force field38 was used together with some special aminoxy parameters.39 In Figure 2, configurations of the hexapeptide complexing a chloride and a sodium ion are shown. In the simulations using hidden restraints, six distance restraints from the amine hydrogens to the ion with a length of 0.35 nm and another six from the carbonyl oxygens to the ion with a length of 0.22 nm were used. During the thermodynamic integration, the force constants of the hidden distance restraints corresponding to the complexated anion (N-H‚‚‚Cl-) were linearly changed from 1.0 to 0.0 kJ mol-1 nm-2 and the ones corresponding to the complexated cation (CdO‚‚‚Na+) from 0.0 to 1.0 kJ mol-1 nm-2 for an alchemical change of Cl- to Na+ and the other way round for the change of Na+ to Cl-. Note that, through the prefactor from eq 1, the restraint potential energy
J. Phys. Chem. B, Vol. 110, No. 16, 2006 8491
Figure 3. Average derivative of the Hamiltonian with respect to the (alchemical) pathway parameter λ, obtained from 250 ps simulations (after 50 ps equilibration) at 11 discrete λ values for the cyclic aminoxyhexapeptide ion complex in chloroform. In the upper half (panels A and B) results from unrestrained simulations are shown, the lower half (panels C and D) represents simulations including hidden restraints. On the left side (panels A and C) thermodynamic integration along λ representing a change from Cl- (λ ) 0) to Na+ (λ ) 1) was carried out, whereas the right side (panels B and D) corresponds to the reverse transformation from Na+ (λ ) 0) to Cl- (λ ) 1).
TABLE 1: Numerically Integrated Free Energy Difference ∆G, in kJ/mol, Obtained by Multi-Configurational Thermodynamic Integration of the Cyclic Aminoxy-Hexapeptide Ion Complex in Chloroform, Where the Ion Was Changed from Cl- to Na+ and Back. Averages 〈DH/Dλ〉λ Were Obtained from 250 ps of Simulations (after 50 ps of equilibration) at 11 Discrete λ Valuesa process Cl-
Na+
f Na+ f Cl-
hidden restraints
unrestrained
-149 ( 16 150 ( 13
134 ( 10 245 ( 11
a Error estimates result from block averaging and extrapolating the block length to infinity.47
function is still scaled to zero in the end states, independent of the force constants. The ensemble average of the partial derivative of the Hamiltonian with respect to the (alchemical) λ-dependent pathway from Cl- to Na+ and backward is depicted in Figure 3, and the (numerically) integrated free energy differences ∆G are given in Table 1. In Figure 4, showing the averages of the distances between the backbone carbonyl oxygen and amid hydrogen atoms and the ion during the thermodynamic integration, the problems arising in the absence of distance restraints (panels A and B) are evident. Instead of undergoing the change necessary to switch between the cation and anion complexing conformation, the ion is just expelled from the peptide. On typical MD simulation time-scales, the ion will not find its way back to form a stable complex. Therefore, without distance restraints, no reliable free energy difference can be obtained. The hidden restraints could be successfully applied to stabilize an otherwise unusable transition path and to calculate a free energy difference between the two states (panels C and D). 3.2. Relative Stabilities of Hexopyranose in 4C1 vs 1C4 Conformation. A new GROMOS force-field parameter set, 45A4, has become available for explicit solvent (water) simulations of hexopyranose-based carbohydrates.40 It was obtained according to the following procedure: (1) reassigning the atomic partial charges based on a fit to the quantum-mechanical electrostatic potential around a trisaccharide; (2) refining the torsional potential energy function parameters associated with the rotations of the hydroxymethyl, hydroxyl, and anomeric
8492 J. Phys. Chem. B, Vol. 110, No. 16, 2006
Christen et al.
Figure 4. Average distances between the ion and the six carbonyl oxygen atoms (in red) and the six amide hydrogen atoms (in black) of the cyclic aminoxy hexapeptide ion complex are shown, obtained from 250 ps simulation (after 50 ps equilibration) at 11 discrete λ values during a multiconfigurational thermodynamic integration from Cl- (λ ) 0) to Na+ (λ ) 1) (left side, panels A and C) and backward from Na+ (λ ) 0) to Cl(λ ) 1) (right side, panels B and D), without restraints (panels A and B) and with (hidden) restraints (panels C and D).
TABLE 2: Values (in degrees) of the Six Dihedral Angles of the Sugar Ring in β-D-glucopyranoside for the 4C1 and the 1C Conformationsa 4
a
dihedral angle
4C
O5 - C1 - C2 - C3 C1 - C2 - C3 - C4 C2 - C3 - C4 - C5 C3 - C4 - C5 - O5 C4 - C5 - O5 - C1 C5 - O5 - C1 - C2
52 -57 49 -40 43 -47
1
1C
4
-50 55 -57 50 -50 52
The atom names are specified in Figure 6.
Figure 5. 4C1 (A) and 1C4 (B) conformations of β-D-glucopyranoside.
alkoxy groups by fitting to corresponding quantum-mechanical profiles for hexopyranosides; (3) adapting the rotational potential energy function parameters determining the ring conformation so as to stabilize the (experimentally predominant) 4C1 chair conformation (Figure 5). In unrestrained simulations, starting from the 1C4 conformation, isomerization to the dominant 4C1 conformation occurred within, at most, 2 ns. Starting from the 4C1 conformation, no isomerization occurred during 5 ns.40 To calculate the free energy difference between the two conformations hidden (dihedral-angle) restraints were applied to force a transition from 4C to 1C . The restraints were applied to the ring dihedral 1 4 angles, with reference values φ0,A and φ0,B in the two end states 4C and 1C according to Table 2. The atom numbering is 1 4 specified in Figure 6. During the isomerization, the dihedral angles will go through the syn conformation (torsional dihedral angles of 0°). Therefore, the maximum value of the restraining potential energy function was set to be reached at φmax ) k -180°. 200 ps of MD simulations (including first 50 ps of equilibration) were performed in explicit SPC water,41 at 11
Figure 6. Atom numbering of the β-D-glucopyranoside.
discrete λ values (λ ) 0.0, 0.1, 0.2, ..., 1.0). The exponents m and n of the hidden restraints in eq 1 were chosen as m ) 2 and n ) 2, and the force constants Kdih ) 10 kJ/rad2 with all six individual weight factors equal to 1.0. A temperature of 300 K (separate coupling of solute and solvent, τT ) 0.1 ps) and pressure of 1 atm (τP ) 0.5 ps, kI ) 4.575 10-4 (kJmol-1nm-3)-1 using isotropic scaling of coordinates) were maintained by weak coupling,33 bond lengths were constrained using the SHAKE algorithm34 (with a relative geometric tolerance of 10-4) and nonbonded interactions were handled using a triple-range cutoff scheme.35 Within a short-range cutoff radius of 0.8 nm, the interactions were evaluated every time step based on a pairlist
Sampling of Rare Events Using Hidden Restraints
J. Phys. Chem. B, Vol. 110, No. 16, 2006 8493
Figure 7. Average values of the ring dihedral angles (see Table 2) during a multi-configurational thermodynamic integration of β-D-glucopyranoside in explicit SPC water with 150 ps simulation (and 50 ps equilibration) at each discrete λ point. The dihedral-angle averages for a change from the 4C (λ ) 0) to the 1C (λ ) 1) conformation are indicated by circles, connected with dotted lines, the ones for the backward change from the 1C 1 4 4 to the 4C1 conformation are indicated by crosses. The order of the dihedral-angles as given in Table 2 corresponds to the colors black, red, green, blue, yellow, and indigo.
Figure 8. Average derivative of the Hamiltonian with respect to the pathway coordinate λ during a multi-configurational thermodynamic integration of β-D-glucopyranoside in explicit SPC water with 150 ps simulation (after 50 ps equilibration) at each discrete λ point. The dashed line with bold error bars corresponds to a conformational change from 4C (λ ) 0) to 1C (λ ) 1), the dotted line and thin error bars to the 1 4 reverse isomerization.
recalculated every five time steps. The intermediate-range interactions up to a long-range cutoff radius of 1.4 nm were evaluated simultaneously with each pairlist update, and assumed constant in between. To account for electrostatic interactions beyond the long-range cutoff radius, a reaction-field approximation36 was applied, using a relative dielectric permittivity of 66 for the solvent.42 Figure 7 shows the average values of the dihedral angles in the sugar ring at the 11 discrete λ values and the corresponding derivative of the Hamiltonian is depicted in Figure 8. Using thermodynamic integration (eq 2), the free energy difference between the two conformations is calculated to be 4.6 ( 3.4 kJ/mol for the change from 4C1 to 1C4 and -5.6 (
Figure 9. Potential of mean force with respect to the pathway coordinate λ of β-D-glucopyranoside (λ ) 0: 4C1 and λ ) 1: 1C4) in explicit SPC water with constrained dihedral angles (see Table 2), with (at least) 50 ps simulation (after 5 ps equilibration) at each discrete λ point (in total 48 points). No mean force could be calculated for λ-values of 0.48, 0.50 and 0.52 because the Lagrange multipliers lk of the dihedral-angle constraints approach zero and the iterative solution of eq 58 does not converge.
2.5 kJ/mol backward. The (internal) potential energy difference between the two conformations had been calculated to be 35 kJ/mol,40 so interactions with the solvent and entropy are lowering the difference significantly. For comparison, the potential of mean force of the transition from the 4C1 to the 1C4 state of the β-D-glucopyranoside was also calculated using dihedral-angle constraints (Figure 9). The ring dihedral angles (see Table 2) were changed linearly with the pathway parameter λ and 51 discrete λ-points were used. At each λ point a short equilibration (5 ps) was followed by (at least) 50 ps of simulation over which the contributions of the constraint forces to the potential of mean force (see Appendix,
8494 J. Phys. Chem. B, Vol. 110, No. 16, 2006
Christen et al.
eq 59) were averaged. The iterative procedure applied to solve the coupled equations with terms up to eighth power in λ generally converges well, even in case the dihedral constraints form a cycle (in the sugar ring). But for specific values of the dihedral angle (90°, 180°, 270°, and 360°), the linearized equations for the Lagrange multipliers (eq 58) yield lk’s that are close to zero. In those cases, the iterative procedure to solve for the constraint forces does not converge. Therefore, no mean force could be calculated for λ-values 0.48, 0.50, and 0.52. Integrating the potential of mean force leads to a relative free energy difference between the 4C1 and the 1C4 state of 8.6 ( 10.2 kJ/mol. Almost twice the amount of simulation time was spent to get this potential of mean force compared to the simulations using dihedral-angle restraints. Yet, the curve is quite rough and the uncertainty is quite high. Therefore, it is not possible to estimate the influence the constraining of the end-state has on the final result.
state 4C1 to 1C4 was compared to the λ-derivative of the Hamiltonian in multi-configurational thermodynamic integration using hidden restraints. Using hidden restraints was straightforward, whereas dihedral-angle constraining turned out to be difficult because the iterative procedure to obtain the Lagrange multipliers that determine the contributions of the constraint forces to the potential of mean force does not converge for all dihedral-angle values. Also, even though more simulation time was spent, the error estimates47 are considerably larger. We conclude that use of hidden restraints constitutes an efficient means to obtain free energy differences between states that are rarely sampled in unrestrained simulations.
4. Discussion Calculating relative free energies between pairs of different (meta)stable states may be quite demanding. On one hand, the two states may be separated by high potential energy barriers which have to be overcome. On the other hand, very extensive sampling may be necessary to have sufficient data on both states. Both problems can be mediated by applying restraints. But at the same time, the end states of the thermodynamic integration may be influenced by the restraints, making the free energy difference between the end states dependent on the path chosen. This pathway dependence can be avoided by the use of hidden restraints which are characterized by zero restraining energy and forces in the end states. Using hidden distance restraints, the relative free energy difference of changing from an anion to a cation in a cyclic D,L-R-aminoxy-hexapeptide-ion complex, accompanied by large structural changes in the peptide, could be calculated by guiding the simulation along a transition path. Application of hidden dihedral-angle restraints enforced the isomerization of a β-D-glucopyranoside from a 4C1 to a 1C4 conformation and backward. The former conformational change (from the more to the less stable structure), was not observed during 5 ns of unrestrained MD simulations. Additionally, harmonic dihedral-angle restraints have been enhanced by an additional parameter φmax to define the angle k with maximum restraint potential energy. Using this parameter, the direction of rotation around the dihedral angle can be controlled. In contrast to standard potential of mean-force calculations,4,43-46 hidden restraints allow relative free energy calculations of unrestrained end states, which makes the results independent of the pathway chosen. As the restraints are only used to guide the simulation and most of the time only with weak influence, the simulated system retains a certain flexibility around the specified transition path. On the other hand, if a steep barrier must be crossed, the system will first lag behind the pathway parameter λ, then suddenly cross the barrier and catch up. During the leap over the barrier, no equilibrium (average) derivative of the Hamiltonian with respect to λ can be obtained. This can be seen by peaks or jumps in the free energy derivative and by large root-mean-square deviations (or error estimates) of its average values at discrete λ points. This problem may be alleviated by lowering the high barriers that separate the two states during the thermodynamic integration, if the force-field terms contributing to these energy barriers can be identified. The potential of mean force obtained using dihedral-angle constraints of the transition of β-D-glucopyranoside from the
5. Appendix
Acknowledgment. Financial support by the National Center of Competence in Research (NCCR) Structural Biology of the Swiss National Science Foundation (SNSF) is gratefully acknowledged.
A potential of mean force can alternatively be obtained by constraining the system to a particular λ-value, calculating the derivative of the free energy with respect to λ (dF/dλ) from the constraint forces (f c) and repeat this for the range of λ-values connecting states A and B of the system. Below we present expressions for dF/dλ and f c when applying distance constraints5 or dihedral-angle constraints. 5.1. Distance Constraints. A set of k ) 1, 2, ..., Nc distance constraints between atoms with positions rk1 and rk2 can be written as
σk(r; rk1k20(λ)) ≡ rk1k22 - (rk1k20(λ))2 ) 0, k ) 1, 2, ..., Nc (20) where rk1k2 ≡ rk1 - rk2, k ≡ (k1, k2), and the distance
rk1k2 ≡ |rk1k2| ≡ ((xk1 - xk2)2 + (yk1 - yk2)2 + (zk1 - zk2)2)1/2, (21) is constrained to the λ-dependent value 0,B rk1k20(λ) ) (1 - λ)r0,A k1k2 + λrk1k2,
(22)
in which rk0,A is the distance constraint value in state A and rk0,B 1k 2 1k2 that in state B. Furthermore, we use the notation r ) (r1, r2, ..., rN) for a configuration of N atoms. The use of rk1k22 to define the constraint σk in eq 20 instead of rk1k2 leads to simpler equations to obtain the constraint forces and their contribution to the potential of mean force. Newton’s equations of motion for N atoms with masses mi including a potential energy function V(r) and the constraints σk multiplied with the Lagrange multipliers lk(t) are
d2ri(t) mi dt
2
)-
∂
Nc
∑
(V (r) + lk(t)σk(r; rk1k20(λ))), ∂ri k)1 i ) 1, 2, ..., N. (23)
The Lagrange multipliers lk(t) are to be determined such that the condition given in eq 20 is satisfied. The first term on the right in eq 23 represents the unconstrained force f uc i (t) derived from the interaction function V (r) and the second term represents the (yet unknown) constraint force f ci (t),
Sampling of Rare Events Using Hidden Restraints
f
c i (t)
)-
Nc
∂σk(r; rk1k20(λ))
k)1
∂ri
∑lk(t)
J. Phys. Chem. B, Vol. 110, No. 16, 2006 8495 constrained simulation at the value λ. The first term on the right contains the possible contribution of the kinetic energy K(p; r), the second the contribution of the unconstrained interaction terms and the third the contribution of the constraint forces, which can be expressed for the k-th constraint as
Nc
) -2
∑(δik
k)1
1
- δik2)lk(t)rk1k2(t).
(24)
The leapfrog scheme31 to integrate Newton’s equations of motion using a time step ∆t yields for the unconstrained positions at time tn + ∆t,
dFck(λ) 0,A ) - 2〈lk〉λrk1k20(λ)(r0,B k1k2 - rk1k2). dλ
The total contribution of the Nc constraints to the potential of mean force is then
-1 uc 2 ruc i (tn + ∆t) ) ri(tn) + Wi(tn - ∆t/2)∆t + mi f i (tn)(∆t) , (25)
where the atomic velocities are indicated by Wi. The constrained positions at time tn + ∆t are related to the constraint forces through 2 -1 c ri(tn + ∆t) ) ruc i (tn + ∆t) + mi f i (tn)(∆t) ,
(26)
and should satisfy the constraint eq 20,
σk(r(tn + ∆t); rk1k20(λ)) ) 0, k ) 1, 2, ..., Nc
(27)
which yields the following equations for the Lagrange multipliers lk(tn),
dFc(λ) dλ
-1
+ ∆t) - 2mk1 (∆t)
2
1 1
1 2
1 2
-1
+ ∆t) + 2mk2 (∆t)
δk2k2′)lk'(tn)rk1′k2′(tn)]2 - [(1 -
2
∑ (δk k ′ -
k′)1 λ)r0,A k1k2 +
rk6k3 ≡ rk3k2 × rk3k4,
(36)
sign(φk) ) sign(rk1k2‚rk6k3),
(37)
After linearization one finds 0,B 2 uc 2 [(1 - λ)r0,A k1k2 + λrk1k2] - [rk1k2(tn + ∆t)]
the Lagrange multipliers at time step tn. The derivative of the free energy F(λ) with respect to λ for a system including distance constraints is5
dλ
〈〉 〈 〉 〈∑ ∂λ
λ
∂V ∂λ
+
λ
∂
Nc
∂λk)1
following the IUPAC-IUB convention.48 Since 0 e arccos e π, we have
-π e φk e π.
(38)
Because of the occurrence of the arccos function in the definition of the dihedral angle φk, the constraints are to be formulated in terms of cos(φk). The occurrence of the square root functions in the distances |rk5k2| and |rk6k3| in the denominator of eq 34 suggests that the use of cos2(φk) will simplify the expressions. Thus, we consider a set of Nc dihedral-angle constraints
σk(φk(r); φk0(λ)) ≡ cos2(φk(r)) - cos2(φk0(λ)) ) 0, k ) 1, 2, ..., Nc (39) where the angle φk(r) is constrained to the λ-dependent value
, (30)
- 4(∆t)2(mk1-1 + mk2-1)rk1k2(tn)·ruc k1k2(tn + ∆t)
+
(34)
2 λr0,B k1k2] ) 0
A 0,B 2 [(1 - λ)r0, k1k2 + λrk1k2] ) 0. (29)
∂K
,
(35)
2 -1 + mk2-1)rk1k2(tn)]2 [ruc k1k2(tn + ∆t) - lk(tn)2(∆t) (mk1
)
)
|rk5k2||rk6k3|
rk5k2 ≡ rk1k2 × rk3k2,
This is a set of Nc quadratic equations in the Nc unknowns lk(tn). It can be solved by linearization (neglect of terms quadratic in lk(tn)) followed by matrix inversion or by sequentially solving the linearized equations for each constraint omitting the coupling between the different constraints (equations), and iterating through all the equations until the lk(tn) converge to a consistent value. The latter method is used in the procedure SHAKE.34 The quadratic, decoupled equation for the Lagrange multiplier lk(tn) is
dF
(
rk5k2‚rk6k3
2 1
k ) 1, 2, ..., Nc. (28)
lk(tn) )
(33)
where and
Nc
ruc k2 (tn
∑
φk ) sign(φk)arccos
∑ (δk k ′ - δk k ′)lk'(tn)rk ′k ′(tn) -
k′)1
dFck(λ) ) . k)1 dλ Nc
5.2. Dihedral-Angle Constraints. For dihedral-angle constraints, the derivation of the expressions for the constraint forces f c and their contribution dF c/dλ to the free energy F(λ) follows the same lines as that for the distance constraints. However, due to the not very simple dependence of a dihedral angle φk(r) upon the positions rk1, rk2, rk3, and rk4 of its four constituting atoms k1, k2, k3, and k4 (i.e., k1-k2-k3-k4), the formulas become much more complicated. Expressions for φk(rk1, rk2, rk3, rk4) can be found in the literature,15-17
Nc
[ruc k1 (tn
(32)
〉
lkσk(r; rk1k20(λ)) . (31) λ
The symbol 〈 ... 〉λ denotes an ensemble average over the
0,B φk0(λ) ) (1 - λ)φ0,A k + λφk ,
(40)
0,B in which φ0,A k is the φk-value in state A and φk that in state B. Newton’s equations of motion for N atoms become
Nc d2ri(t) ∂ ) - (V(r) + lk(t)σk(φk(r); φk0(λ))), mi ∂ri k)1 dt2 i ) 1, 2, ..., N (41)
∑
where the Lagrange multipliers lk(t) are to be determined such
8496 J. Phys. Chem. B, Vol. 110, No. 16, 2006
Christen et al.
that the condition given in eq 39 is satisfied. The second term on the right in eq 41 represents the (yet unknown) constraint forces, 0
Nc
f ci (t) ) -
∂σk(φk(r); φk (λ))
∑lk(t) k)1
∂ri
-1 -1 [ruc k1k2(tn + ∆t) × (mk3 ak3(tn) - mk2 ak2(tn)) -
∂φk(r) ) + lk(t)2cos(φk)sin(φk) . ∂ri k)1
∑
[( [(
δik3
|rk3k2|2
rk3k4‚rk3k2 |rk3k2|2
-1
) )
-1
|rk3k2|
|rk5k2|
or using a shorter notation bk1k2k3(tn + ∆t) for the last factor uc rk5k2(tn + ∆t) ) ruc k1k2(tn + ∆t) × rk3k2(tn + ∆t) +
|rk3k2| ∂φk(r) ) δik1 rk k + ∂ri |r |2 5 2 rk1k2‚rk3k2
-1 -1 ruc k3k2(tn + ∆t) × (mk1 ak1(tn) - mk2 ak2(tn))] (48)
(42)
Expressions for (∂φk/∂ri) can be found in the literature too,15-17
δik2
uc rk5k2(tn + ∆t) ) ruc k1k2(tn + ∆t) × rk3k2(tn + ∆t) +
lk(tn)sin(2φk(tn))(∆t)2
Nc
k5k2
cross products in eqs 35 and 36 and linearizing the resulting expressions yields
r + 2 k5k2
|rk3k2|
rk6k3 +
|rk6k3|2
lk(tn)sin(2φk(tn))(∆t)2bk1k2k3(tn, tn + ∆t), (49)
] ]
rk3k4‚rk3k2 |rk3k2|
rk k |rk3k2|2 |rk6k3|2 6 3
rk1k2·rk3k2 |rk3k2|
rk k |rk3k2|2 |rk5k2|2 5 2 |rk3k2| δik4 rk k . (43) |rk6k3|2 6 3
To shorten the expressions, we denote the four terms in eq 43, apart from the Kronecker delta’s, by ak1, ak2, ak3, and ak4, respectively. Then we have
and uc rk6k3(tn + ∆t) ) ruc k3k2(tn + ∆t) × rk3k4(tn + ∆t) +
lk(tn)sin(2φk(tn))(∆t)2 -1 -1 [ruc k3k2(tn + ∆t) × (mk3 ak3(tn) - mk4 ak4(tn)) -1 -1 ruc k3k4(tn + ∆t) × (mk3 ak3(tn) - mk2 ak2(tn))] (50)
or using the shorter notation uc rk6k3(tn + ∆t) ) ruc k3k2(tn + ∆t) × rk3k4(tn + ∆t) +
lk(tn)sin(2φk(tn))(∆t)2bk2k3k4(tn, tn + ∆t). (51)
Nc
fci (t)
)
∑lk(t)sin(2φk(t)[δik ak (t) + δik ak (t) + δik ak (t) +
k)1
1
1
2
2
3
3
δik4ak4(t)]. (44) The leapfrog scheme yields the unconstrained positions ruc i (tn + ∆t) from eq 25. The constrained positions ri(tn + ∆t) are related to the constraint forces (eq 44) through eq 26 and should satisfy the constraint eqs 39,
cos2(φk(r(tn + ∆t)) - cos2(φk0(λ)) ) 0, k ) 1, 2, ..., Nc, (45) or using eq 34,
[
rk5k2(tn + ∆t)‚rk6k3(tn + ∆t)
|rk5k2(tn + ∆t)||rk6k3(tn + ∆t)|
]
The scalar product in the numerator of the first term in eq 46 becomes after linearization uc rk5k2(tn + ∆t)‚rk6k3(tn + ∆t) ) (ruc k1k2(tn + ∆t) × rk3k2(tn + uc 2 ∆t))‚(ruc k3k2(tn + ∆t) × rk3k4(tn + ∆t)) + lk(tn)sin(2φk(tn))(∆t) uc [(ruc k1k2(tn + ∆t) × rk3k2(tn + ∆t))·bk2k3k4(tn, tn + ∆t) + uc (ruc k3k2(tn + ∆t) × rk3k4(tn + ∆t))‚bk1k2k3(tn, tn + ∆t)] (52)
or in a shorter notation
2
- cos2(φk0(λ)) ) 0.
(46)
Since rk5k2(tn + ∆t) and rk6k3(tn + ∆t) are each quadratic in the Lagrange multipliers lk(tn), both the numerator and the denominator of the left term in eq 46 contain powers of up to eight of the lk(tn). Thus a set of Nc equations consisting of terms containing up to powers of eight of the unknowns lk(tn) is to be solved. As for the case of distance constraints, this is achieved by linearizing the equations for each constraint, omitting the coupling between the different constraints (equations), and iterating through all Nc equations until the lk(tn) converge to a consistent value. Using eqs 26 and 44 we find for the k-th constraint
rk5k2(tn + ∆t)‚rk6k3(tn + ∆t) ) ck1k2k3k4(tn + ∆t) + lk(tn)sin(2φk(tn))(∆t)2dk1k2k3k4(tn, tn + ∆t). (53) The square becomes after linearization
(rk5k2(tn + ∆t)·rk6k3(tn + ∆t))2 ) (ck1k2k3k4(tn + ∆t))2 + 2lk(tn)sin(2φk(tn))(∆t)2ck1k2k3k4(tn + ∆t)dk1k2k3k4(tn, tn + ∆t). (54) The factors in the denominator of the first term in eq 46 becomes uc 2 |rk5k2(tn + ∆t)|2 ) (ruc k1k2(tn + ∆t) × rk3k2(tn + ∆t)) +
2lk(tn)sin(2φk(tn))(∆t)2
rk1k2(tn + ∆t) ) ruc k1k2(tn + ∆t) + lk(tn)sin(2φk(tn))(∆t)2(mk1-1ak1(tn) - mk2-1ak2(tn)) (47) and likewise for rk3k2(tn + ∆t) and rk3k4(tn + ∆t). Building the
uc (ruc k1k2(tn + ∆t) × rk3k2(tn + ∆t))‚bk1k2k3(tn + ∆t)
and
(55)
Sampling of Rare Events Using Hidden Restraints
J. Phys. Chem. B, Vol. 110, No. 16, 2006 8497
uc 2 |rk6k3(tn + ∆t)|2 ) (ruc k3k2(tn + ∆t) × rk3k4(tn + ∆t)) +
2lk(tn)sin(2φk(tn))(∆t)2 uc (ruc k3k2(tn + ∆t) × rk3k4(tn + ∆t))‚bk2k3k4(tn + ∆t).
(56)
The linearized denominator of the first term in eq 46 is then uc |rk5k2(tn + ∆t)|2|rk6k3(tn + ∆t)|2 ) (ruc k1k2(tn + ∆t) × rk3k2(tn + uc 2 ∆t))2(ruc k3k2(tn + ∆t) × rk3k4(tn + ∆t)) +
2lk(tn)sin(2φk(tn))(∆t)2 uc 2 [(ruc k1k2(tn + ∆t) × rk3k2(tn + ∆t)) uc (ruc k3k2(tn + ∆t) × rk3k4(tn + ∆t))‚bk2k3k4(tn, tn + ∆t) + uc 2 (ruc k3k2(tn + ∆t) × rk3k4(tn + ∆t)) uc (ruc k1k2(tn + ∆t) × rk3k2(tn + ∆t))‚bk1k2k3(tn, tn + ∆t)] (57)
Finally, the equation for the Lagrange multiplier of the k-th constraint becomes uc 2 lk(tn) ) [cos2(φk0(λ))(ruc k1k2(tn + ∆t) × rk3k2(tn + ∆t)) uc 2 2 (ruc k3k2(tn + ∆t) × rk3k4(tn + ∆t)) - (ck1k2k3k4(tn + ∆t)) ]
[2sin(2φk(tn))(∆t)2[ck1k2k3k4(tn + ∆t)dk1k2k3k4(tn + ∆t) -cos2(φk0(λ)) uc 2 ((ruc k1k2(tn + ∆t) × rk3k2(tn + ∆t)) uc (ruc k3k2(tn + ∆t) × rk3k4(tn + ∆t))‚bk2k3k4(tn, tn + ∆t) uc 2 + (ruc k3k2(tn + ∆t) × rk3k4(tn + ∆t)) uc -1 (ruc (58) k1k2(tn + ∆t) × rk3k2(tn + ∆t))·bk1k2k3(tn, tn + ∆t))]]
The derivative of the contribution of the constraint forces to the free energy for the k-th constraint becomes
dF ck(λ) 0,A ) 〈lk〉λ sin(2φk0(λ))(φ0,B k - φk ). dλ
(59)
We note that the expressions given in this Appendix for the application of dihedral-angle constraints are different from the formalism presented in ref 49, which is based on matrix inversion. References and Notes (1) van Gunsteren, W. F.; Huber, T.; Torda, A. E. In Conf. Proc. European Conference on Computational Chemistry (E.C.C.C 1); American Institute of Physics (A.I.P.): College Park, MD, 1995; Vol. 330, pp 253268. (2) Schlitter, J.; Swegat, W.; Mu¨lders, T. J. Mol. Modell. 2001, 7(6), 171-177. (3) Bolhuis, P. G.; Chandler, D.; Dellago, C.; Geissler, P. L. Annu. ReV. Phys. Chem. 2002, 53, 291-318. (4) Hill, T. L. Statistical Mechanics; McGraw-Hill: New York, 1956. (5) van Gunsteren, W. F.; Beutler, T. C.; Fraternali, F.; King, P. M.; Mark, A. E.; Smith, P. E. In Computer simulation of biomolecular systems,
theoretical and experimental applications; van Gunsteren, W. F., Weiner, P., Wilkinson, A. J., Eds.; ESCOM Science Publishers: B. V.: Leiden, The Netherlands, 1993; Vol. 2, pp 315-367. (6) den Otter, W. K.; Briels, W. J. J. Chem. Phys. 1998, 109(11), 41394146. (7) Sprik, M.; Ciccotti, G. J. Chem. Phys. 1998, 109(18), 7737-7744. (8) Torrie, G. M.; Valleau, J. P. J. Comput. Phys. 1977, 23, 187-199. (9) Billeter, S. R.; van Gunsteren, W. F. J. Phys. Chem. A 2000, 104, 3276-3286. (10) Alder, B. J.; Wainwright, T. E. J. Chem. Phys. 1957, 27, 12081209. (11) Kirkwood, J. G. J. Chem. Phys. 1935, 3, 300-313. (12) Mark, A. E.; van Gunsteren, W. F. In New PerspectiVes in Drug Design; Dean, P. M., Jolles, G., Newton, C. G., Eds.; Academic Press Ltd.: Turnberry, Scotland, 1995; pages 185-200. (13) van Gunsteren, W. F.; Daura, X.; Mark, A. E. HelV. Chim. Acta 2002, 85, 3113-3129. (14) Ramachandran, G.; Sasisekaran, V. AdV. Prot. Chem. 1968, 23, 283-437. (15) van Schaik, R. C.; Berendsen, H. J. C.; Torda, A. E.; van Gunsteren, W. F. J. Mol. Biol. 1993, 234, 751-762. (16) Bekker, H.; Berendsen, H. J. C.; van Gunsteren, W. F. J. Comput. Chem. 1995, 16, 527-533. (17) van Gunsteren, W. F.; Billeter, S. R.; Eising, A. A.; Hu¨nenberger, P. H.; Kru¨ger, P.; Mark, A. E.; Scott, W. R. P.; Tironi, I. G. Biomolecular simulation: The GROMOS96 manual and user guide; Verlag der Fachvereine: Zu¨rich, 1996. (18) Yang, D.; Qu, J.; Li, W.; Zhang, Y.-H.; Ren, Y.; Wang, D.-P.; Wu, Y.-D. J. Am. Chem. Soc. 2002, 124(42), 12410-12411. (19) Gellman, S. H. Acc. Chem. Res. 1998, 31, 173-180. (20) Seebach, D.; Matthews, J. L. Chem. Commun. 1997, 79, 20152022. (21) Hill, D. J.; Mio, M. J.; Prince, R. B.; Hughes, T. S.; Moore, J. S. Chem. ReV. 2001, 101, 3893-4011. (22) Cubberley, M. S.; Iverson, B. L. Curr. Op. Chem. Bio. 2001, 5(6), 650-653. (23) Ghadiri, M. R.; Granja, J. R.; Milligan, R. A.; McRee, D. E.; Khazanovich, N. Nature 1993, 366(6453), 324-327. (24) Ghadiri, M. R.; Granja, J. R.; Bu¨hler, L. K. Nature 1994, 369(6478), 301-304. (25) Fernandez-Lopez, S.; Kim, H.-S.; Choi, E. C.; Delgado, M.; Granja, J. R.; Khasanov, A.; Kra¨henbu¨hl, K.; Long, G.; Weinberger, D. A.; Wilcoxen, K. M.; Ghadiri, M. R. Nature 2001, 412(6845), 452-455. (26) Clark, T. D.; Bu¨hler, L. K.; Ghadiri, M. R. J. Am. Chem. Soc. 1998, 120(4), 651-656. (27) Seebach, D.; Matthews, J. L.; Meden, A.; Wessels, T.; Baerlocher, C.; McCusker, L. B. HelV. Chim. Acta 1997, 80(1), 173-182. (28) Bong, D. T.; Clark, T. D.; Granja, J. R.; Ghadiri, M. R. Angew. Chem., Int. Ed. 2001, 40(6), 988-1011. (29) Yang, D.; Ng, F.-F.; Li, Z.-J.; Wu, Y. D.; Chan, K. W. K.; Wang, D.-P. J. Am. Chem. Soc. 1996, 118(40), 9794-9795. (30) Yang, D.; Li, B.; Ng, F.-F.; Yan, Y.-L.; Qu, J.; Wu, Y.-D. J. Org. Chem. 2001, 66(22), 7303-7312. (31) Hockney, R. W. Methods Comput. Phys. 1970, 9, 136-211. (32) Tironi, I. G.; van Gunsteren, W. F. Mol. Phys. 1994, 83, 381403. (33) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684-3690. (34) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327-341. (35) van Gunsteren, W. F.; Berendsen, H. J. C. Angew. Chem., Int. Ed. Engl. 1990, 29, 9, 992-1023. (36) Tironi, I. G.; Sperb, R.; Smith, P. E.; van Gunsteren, W. F. J. Chem. Phys. 1995, 102(13), 5451-5459. (37) Schuler, L. D.; Daura, X.; van Gunsteren, W. F. J. Comput. Chem. 2001, 22, 1205-1218. (38) Oostenbrink, C.; Villa, A.; Mark, A. E.; van Gunsteren, W. F. J. Comput. Chem. 2004, 25(13), 1656-1676. (39) Peter, C.; Daura, X.; van Gunsteren, W. F. J. Am. Chem. Soc. 2000, 122, 7461-7466. (40) Lins, R. D.; Hu¨nenberger, P. H. J. Comput. Chem. 2005, 26, 14001412. (41) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, The Netherlands, 1981; pp 331-342. (42) Gla¨ttli, A.; Daura, X.; van Gunsteren, W. F. J. Chem. Phys. 2002, 116, 9811-9828. (43) Tobias, D. J.; Sneddon, S. F.; Brooks, C. L., III J. Mol. Biol. 1990, 216(3), 783-796. (44) Lazaridis, T.; Tobias, D. J.; Brooks, C. L., III; Paulaitis, M. E. J. Chem. Phys. 1991, 95(10), 7612-7625. (45) Beutler, T. C.; van Gunsteren, W. F. J. Chem. Phys. 1994, 100, 1492-1497.
8498 J. Phys. Chem. B, Vol. 110, No. 16, 2006 (46) Berendsen, H. J. C. Simu Newsletter 2001, 3, 33-50. (47) Allen, M. P.; Tildesley, D. J. Computer simulation of liquids; Oxford University Press: New York, 1987.
Christen et al. (48) IUPAC. Biochemistry 1970, 9, 3471-3479. (49) Tobias, D. J.; Brooks, C. L., III J. Chem. Phys. 1988, 89(8), 51155127.
