Each enzyme-substrate system was then subjected to six initial equilibration runs of 40 ns at 300 K with different initial velocities, heating from 0 to 300 K as described in refs. (Figures and). This was then followed by five EVB simulations per equilibration run, performed using 51 mapping frames of 200 ps in length each, resulting in a total of 30 individual EVB trajectories for each system, and 600 EVB trajectories over all systems. The starting points for the EVB simulations were generated by assigning new random seeds at the end point of the preceding equilibration run, and performing a further 110 ps of equilibration to generate the new starting structures for the EVB simulations. All equilibration and EVB simulations were performed at the transition state for the proton transfer reaction, based on the valence bond states shown in Figure , as described in our previous work, in order to enable simultaneous propagation of EVB trajectories downhill towards the reactant and the product state. All simulations were performed using a 1 fs timestep, to a total equilibration time of 240 ns per system (4.80 μs over all systems) and a total EVB simulation time of 306 ns per system (6.12 μs over all systems). Note that we used the same EVB parameters (Hij and α) as in our prior studies to describe the deprotonation of substrates DHAP and GAP. These parameters can be found in Table . For a description of these parameters, see ref. , and for a description of how the EVB fitting was performed for these reactions see ref. . All EVB parameters used in this work can be found in the Supporting Information of ref. , and full protocols and simulation details can be found in our previous work on TIM. All energy analysis in this study was performed using the QCalc module of the Q simulation package. VMD version 1.9.1 was used to perform all structural analysis. Representative structures at various stationary points for the TIM-catalyzed deprotonation of DHAP and GAP in wild-type and P166A TIM were obtained by performing clustering analysis, using the method of Daura et al. as implemented in GROMACS version 4.6.5. A cutoff distance of 0.5 Å was used for the analysis, which was based on structures sampled every 2 ps of the simulation time.