Steered MD allows reactions to be simulated that are so improbable on the time scale of molecular dynamics that they would effectively never happen during a simulation. This approach takes advantage of the Jarzynski equality which allows equilibrium properties to be obtained from averages over nonequilibrium processes
[25]. Following the approach of Bustamante
et al.
[19], consider a system at temperature

whose equilibrium state is determined by a control parameter

. Let the system initially be in state

with control parameter

. If the system is evolved via a nonequilibrium process by changing

along a given path

to some final value

, the Jarzynski relationship states that
where

is the free energy difference between equilibrium states

is the work done along the trajectory

is Boltzmann's constant and

denotes an average over an infinite number of such nonequilibrium experiments repeated under the protocol

. While the derivation of the Jarzynski relationship invokes fairly deep statistical physics, the effect of using Jarzynski exponential averaging for PMF determination along a path is to minimize the impact of higher energy (and, therefore, lower probability) steered trajectories in favor of the lower energy (i.e., more probable) trajectories. The definition of the Jarzynski relationship is as an average over an infinite number of trajectories. Averaging over a finite number of trajectories introduces a truncation bias in the computation. Gore,
et al.
[24] have introduced methods that correct for that bias to a large degree, and have also furnished equations for the mean square errors associated with such a bias–corrected Jarzynski average. These methods have been used in place of a more classical definition of confidence limits. Further details are given in the
file S4.
In this study, state

was a snapshot from an equilibrium MD trajectory carried out with a collective progress variable (CV) restrained to its initial value. Multiple such starting structures were obtained by periodically sampling the restrained trajectory. To avoid spurious correlations between SMD replicates started with these snapshots, they were not only obtained from different points in the restrained ensemble trajectory, but also simulated with different values of the pseudo random number generator (PRNG) seed (the system clock in microseconds was used to generate the seed value). The collective variable (a linear combination of interatomic distances) constituted

, and it was varied from initial to final values under a harmonic restraint. Since in this study covalent bonds were being broken, very stiff restraints were required

. Preliminary studies with restraints in the range of

showed that the actual collective variable values achieved sometimes lagged the targeted restrained values. This behavior was followed by discontinuous “jumps” in the CV coordinate as its values “caught up” with the targeted restrained values. These discontinuities complicated the analyses, so

was chosen. This value ensured that the energy of the restraint at all points exceeded the energy of the bonds being broken. Under these conditions, plots of restrained (targeted) CV value versus actual (achieved) CV value were straight lines with slight scatter (governed by the finite width of the underlying harmonic restraint parabola) but no large discontinuities. In all cases, the final value of the CV was attained at the end of the SMD simulations to within an interval consistent with the magnitude of the steering restraint.
illustrates the initial active site quantum region structure for the SMD simulation analyzed in . The atom names are those referred to in the collective variable definition in , Panel B. Several distances are also defined:

,

, and

, where

is the distance between atoms

and

. The starting value of the SMD CV was its value at the end of the 1 ns QM/MM relaxation run (8.92 Å). The ending value (3.52 Å) was taken to be the sum of the bond lengths:

. In this and similar simulation setups, the quantum water molecule(s) were restrained to be within a flat well potential spherical volume around

with weak

harmonic restraint walls. Further details are given in the
file S3.
shows the PMF determination for carbinolamine formation mediated by an intermediate water acting as a proton transfer reagent. Thirty SMD replicates (colored lines in Panel A) were Jarzynski averaged (solid black line, both panels.) The definition of the elements of Panel B are given in the legend to this figure. The atom names in the CV are identified in .
is a stereogram illustrating a typical initial active site quantum region for the SMD simulations further analyzed in (imine formation without involvement of an intermediate water). The atom names referred to in the CV definition in , Panel B are identified. Interatomic distances are defined similarly to those in .
shows the PMF determination for direct proton transfer carbinolamine formation without water involvement. Work functions for 30 SMD replicates are shown in Panel A, and atom names referred to in the CV defined in Panel B are identified in . The SMD CV end points were 8.47 Å for the starting value (measured at the end of the 1 ns QM/MM relaxation run) and 2.51 Å at the end:

in .
In preparation for SMD studies of imine carbinolamine dehydration, simulations were carried out with a quantum water restrained in the active site by a constant “collective variable”. The starting structure for these simulations was obtained from one of the product structures in the carbinolamine formation study. In these preliminary simulations to generate “snapshots” for SMD, carbinolamine dehydration was sometimes observed without SMD forcing. To follow up on these observations, three simulations were carried out with identical initial structures. In one of the simulations, no dehydration was observed, in a second, dehydration occurred with the participation of the restrained water as an intermediate proton transfer reagent, and in a third, dehydration occurred by direct transfer of the protonated Glu22 hydrogen to the carbinolamine oxygen without participation of the restrained water.
is a stereogram of the quantum region for these unforced carbinolamine dehydration simulations, and this diagram corresponds to the active site analyzed in and . As this trajectory was not “steered” between differing endpoints, the CV did not change over its course. The starting (and ending)

. Note that in this kind of restraint, the positions of the individual atoms may change as only the sum of interatomic distances is restrained. The atoms are shaded as in .
In AMBER, residue and atom names do not change over the course of a simulation. At the end of the run analyzed in , the protonated Glu22 carboxylate hydrogen (

) was found associated with the water oxygen (

) initially present in the active site, as was one of the original water hydrogens (

), the hydrogen in the restrained CV. In addition, the product water was composed of the other water hydrogen (

), not involved in the restrained CV, the carbinolamine

oxygen and

hydrogen. This unambiguously identified the trajectory in as that of unforced carbinolamine dehydration with the assistance of the initial active site water as an intermediate proton transfer reagent. This trajectory also illustrates that the restraint operates only on the sum of distances:

shrank to the equilibrium water

bond length, while

increased to maintain the constant sum.
The axes of , Panel A are proton transfer reaction coordinate variables:

on the abscissa, and

on the ordinate. In these types of proton transfer variables, the value changes from negative to positive as the proton transfers from the donor to the acceptor, with

indicating the hydrogen at the midpoint of the donor, acceptor interatomic distance. Progress along the trajectory is coded in color (as defined in the palette). The two proton transfers associated with water-mediated carbinolamine dehydration occur within 2.5 ps of each other, which is the trajectory sampling interval.
The abscissa of , Panel B is the same as the ordinate of Panel A, whereas the ordinate reflects the “transfer of the

oxygen” (plus its bonded

hydrogen) to the unconstrained water hydrogen (

), forming the water molecule product of the dehydration reaction. As in Panel A, both reactions occur within the sampling dead time of each other. In this panel, as in , Panel B, a proton transfer coordinate is plotted against the “transfer” of the acceptor oxygen “to a proton” to form water. In the initial part of the trajectory, i.e., in the negative regions of both coordinates, these two movements are seen to be almost completely reciprocal to each other, as reflected in the relatively low scatter in the points.
Panel C shows a trajectory course view of the “transfer of the

oxygen” (and its bonded

hydrogen) to the water hydrogen to form the water product. The ordinate in this panel is the same as that of Panel B, and measures the oxygen transfer (the dehydration reaction itself). The red curve shows the distance of the water hydrogen to the dehydrating oxygen. Note that the black curve crosses

(the oxygen is at the numerical midpoint of the coordinate) at the same time the

distance becomes equal to and maintains a typical water

bond length.
In , Panel A, the black line is the distance

in and measures proton transfer from the protonated Glu22 carboxylate (

) to the carbinolamine oxygen (direct general acid catalyzed carbinolamine dehydration). The red line is

, the donor oxygen, acceptor oxygen distance for comparison. The proton transfer, at

of the trajectory duration, is accompanied by relative stability in the donor-acceptor distance. Eventually, separation of the oxygens due to diffusion diminishes the influence on the proton of the departed carboxylate oxygen at

of the trajectory duration. As mentioned above, interrogation of the product structure atom names unambiguously identified the simulated reaction course. Panel B compares the timing of the

proton transfer from

to that of the “

transfer” from

. Both reactions occur within one sampling frame.
, Panel C, shows what would be the proton transfer reaction surface as in , Panel A, for a trajectory in which no carbinolamine dehydration was observed. Note the similarity of the dynamics in this trajectory to the initial dynamics in Panel A of (in the rectangular region bounded by abscissa values of

and

, and ordinate values of

and

).
In addition to spontaneous dehydration of the carbinolamine, unforced, direct proton transfer of the protonated Glu22 hydrogen to the unprotonated N–terminal amine nitrogen was observed. To investigate this event, shows the quantum region active site used in three simulations. The structure is that of the carbinolamine formation reactants. All were started from this same structure (which had been observed to undergo this unforced proton transfer event), and with the same restrained quantum water:

at both the beginning and end of these unforced simulations. Spontaneous proton transfer occurred in one of the three simulations, but without participation of the quantum water. The atom names referred to in the variable definitions in are given in .
In , the black curve shows the reaction coordinate for the transfer of the HE2 protonated Glu22 hydrogen to the N–terminal amine nitrogen over the course of the trajectory. The proton transfer event occurred slightly after the midpoint of the trajectory duration. The red curve measures the distance between the donor and acceptor heavy atoms, and shows the interatomic distance between them decreasing directly before proton transfer. The potential for spontaneous proton transfer occurred when diffusion of the heavy atoms brought them into close proximity. Subsequent to the proton transfer event, the two curves track each other closely. Two other replicates in this series did not display spontaneous proton transfers of this type, showing only the two proton-heavy atom distances oscillating about their initial values (data not shown).