\({\mathrm{\Delta }}C_P^\ddagger\) determined by experiment:

As shown previously6,7, \({\mathrm{\Delta }}C_P^\ddagger\) can be determined by fitting the ln(rate)-versus-temperature plot using MMRT (Fig. 1). For MalL, curvature in this plot is very significant (Fig. 1d), and unrelated to unfolding5. This leads to a negative \({\mathrm{\Delta }}C_P^\ddagger\) value of –11.6 ± 0.4 kJ mol−1 K−1. For KSI, the curvature is less extreme, but still obvious, leading to a small negative \({\mathrm{\Delta }}C_P^\ddagger\) of –0.86 ± 0.1 kJ mol−1 K−1 (Fig. 1c). An important consequence of the \({\mathrm{\Delta }}C_P^\ddagger\) values for each enzyme is the position of the optimum temperature (T opt ) for activity as these parameters are correlated. For example, the large negative \({\mathrm{\Delta }}C_P^\ddagger\) value for MalL dictates the position of T opt at 320 K, whereas the much smaller \({\mathrm{\Delta }}C_P^\ddagger\) value for KSI places the T opt well above 320 K (in absence of protein unfolding). Similarly, the significant curvature of the MalL temperature dependence means that at lower temperatures, the rate approaches zero much faster for MalL than for KSI. Implicit in this approach is the assumption that \({\mathrm{\Delta }}C_P^\ddagger\) is independent of temperature in the temperature range studied, and is justified based on the good fit of the MMRT model to the experimental data.

Heat capacity differences from simulation

Heat capacity differences for enzyme-catalysed reactions can be calculated from \(\Delta \langle\partial H^{2}\rangle^{\ddagger }\) (Eq. (1)). To measure \(\Delta \langle\partial H^{2}\rangle^{\ddagger }\) from simulation, there are two main challenges: (a) the amount of sampling required for the system to define the enthalpy variance, and (b) an accurate and consistent representation of the reactant state (Michaelis complex) and the TS. A statistical thermodynamic analysis of a 1 ms MD simulation of the bovine pancreatic trypsin inhibitor indicated that 10s of microseconds (μs) of simulation may be needed to converge the heat capacity difference between two conformational states18. Sampling on the order of (at least) microseconds is thus expected to be required for reliable identification of heat capacity differences. Such sampling is now routinely feasible with a ‘molecular mechanics’ description of the atoms and their interactions. Molecular mechanics force fields have been developed and optimised over many years19, and can provide a generally good description of the structure and dynamics of proteins and protein–ligand binding20. They are, however, empirical potential functions and typically (for reasons of computational efficiency) lack physically important effects, such as variations in electronic polarisation. This may be related to limitations in their ability to capture details of fast dynamics21. Here, to calculate \({\mathrm{\Delta }}C_P^\ddagger\), we use extensive MD simulation to quantify the change in fluctuation between two states, A and B (Fig. 1a). The difference in heat capacity between these states can be determined by:

$${\mathrm{\Delta }}C_{{\mathrm {A,B}}} = \frac{\langle{\delta H_{\mathrm B}^2\rangle - \langle\delta H_{\mathrm A}^2\rangle}}{{k_{\mathrm B}T^2}}.$$ (2)

To sample the conformational dynamics of the reactant (E–S or RS) and ‘TS’ (E–TS) enzyme complexes consistently, electronically unstable states (e.g. with half-formed bonds involving enzyme residues) should be avoided for the ‘TS’ representation. We, thus, use molecular species that are representative for the TS (i.e. TS analogues), and predict that these will show a similar heat capacity change from the reactant state. This prediction has been demonstrated experimentally for human 5′-methylthioadenosine phosphorylase11. For KSI, a charged enediolate intermediate is formed after the first proton transfer, and this is the key species stabilised by the enzyme for catalysis of the reaction22,23. We use this intermediate state as a proxy for the two enzyme–TS complexes (one for each proton transfer) as the intermediate lies between the two TSs at similar energy. The substrate (5-androstene-3,7-dione) and intermediate complexes (Fig. 2a, Supplementary Fig. 1; Supplementary Fig. 3) were built based on KSI in complex with the inhibitor 5α-estran-3,17-dione (PDB 1OHP). For MalL, we obtained an experimental X-ray structure co-crystallised with a stable transition-state analogue (Supplementary Table 1; Supplementary Fig. 4 and 5) and use this to simulate the thermodynamics of the substrate isomaltose and a close analogue of the transition-state species at the rate determining step24 (Fig. 2b; Supplementary Fig. 2).

Fig. 2 Sampling and \({\mathrm{\Delta }}C_P^\ddagger\) calculation in simulations. a-b Histograms of energies from 50 to 500 ns MD simulations for KSI (a) and MalL (b). Thin lines are individual runs, thick lines are the average of ten runs. Insets show overlay of histograms for both states, and the structures indicate the species simulated (RS reactant state, IS intermediate state, TSA transition state analogue). c Representative structures for the two distinct conformational clusters present in the KSI simulations of both states (reactant state in blue and green, intermediate state in pale blue and green, starting structure in light grey). Box highlights the region with structural differences. d Representative structures for the six main conformational clusters in MalL reactant state simulations and their occupancies (starting structure in light grey). e Variance in energies for the two clusters identified in the KSI simulations, with cluster occupancies (in %) and weighted average variance for both states. f Convergence with moving average-window size of \({\mathrm{\Delta }}C_P^\ddagger\) values calculated for MalL, with value determined from experiment indicated by dotted line (with grey area indicating standard deviation). Error bars indicate the standard deviation of the calculated \({\mathrm{\Delta }}C_P^\ddagger\) values based on the cumulative standard deviation for each state, from 10 independent simulations Full size image

A total of 5 μs of explicit solvent MD simulation was run for KSI and MalL in both the substrate-bound and proxy TS representations over ten replicate simulations for each state. The force-field potential energy was used as an approximation for the system enthalpy, and was recalculated for the protein–ligand system (i.e. all atoms) without explicit water. Considering the variance of the enthalpy is the quantity required for C P calculation (Eq. (1)) and a difference in variance between two states is used to determine \({\mathrm{\Delta }}C_P^\ddagger\) (Eq. (2)) these approximations should be reasonable. Calculating the variance with explicit solvent is problematic because there is no clear criterion for selecting the water molecules that should be included in such a calculation (see Supplementary Note 6). We note that \({\mathrm{\Delta }}C_P^\ddagger\) values calculated with an implicit solvent are qualitatively similar to those reported below for both enzymes (Supplementary Table 5). Note that there is, in principle, an alternative approach to calculating \({\mathrm{\Delta }}C_P^\ddagger\), via the variance in entropy12. Calculating entropy from simulation accurately is much more challenging; however, this may be feasible in the future.

For KSI, the conformational space sampled is limited, with only two distinct structural clusters discernible (see Supplementary Note 5, Supplementary Fig. 6 and Supplementary Table 3). The difference between these clusters is in a small region in the unoccupied monomer (Fig. 2c). The H variance is significantly different between the clusters, however (Fig. 2e). For \({\mathrm{\Delta }}C_P^\ddagger\) calculation, we thus calculate the variance of the clusters separately, with the total variance for each state being the average variance weighted by the cluster occupation (Fig. 2e; Supplementary Note 5).

MalL samples a larger conformational space than KSI, occupying and regularly switching between a number of structural clusters along the simulation trajectories, related primarily to changes in loops surrounding the active site (Fig. 2d; Supplementary Notes 4 and 5, Supplementary Table 2, Supplementary Fig. 8, 10 and 11). Due to the presence of multiple conformational clusters, consideration of the system over the full simulation time overinflates the enthalpy variance. However, calculating variances for each cluster (as for KSI) does not take into account that frequent switches between the distinct conformational states will also contribute to the variance. In addition, several clusters are dominated by one state only (Supplementary Fig. 11). To be independent of clustering and account for switching between conformational substates, enthalpy variance was calculated using a moving window along the simulation trajectory for each simulation, and subsequent averaging. The ‘window’ for the moving average was varied between 5 and 80 ns and calculated \({\mathrm{\Delta }}C_P^\ddagger\) values converge when the window size is between 40 and 80 ns (Fig. 2f, Supplementary Table 4). Thus, the calculated \({\mathrm{\Delta }}C_P^\ddagger\) values for MalL converge on a value of −10.0 ± 1.7 kJ mol−1 K−1 (using a window of 70 ns), which is within the error range of the experimentally determined \({\mathrm{\Delta }}C_P^\ddagger\) value of –11.6 ± 0.4 kJ mol−1 K−1.

Local and global contributions to \({\mathrm{\Delta }}C_P^\ddagger\)

The observation that \({\mathrm{\Delta }}C_P^\ddagger\) values calculated from extensive conformational sampling agree with those determined experimentally allows meaningful analysis of the differences between the two ensembles. A significant part of \({\mathrm{\Delta }}C_P^\ddagger\) resides in the protein backbone (in agreement with experiments that suggest side-chains contribute only a fraction to the total protein heat capacity17; Supplementary Table 5), although the contribution of side-chains cannot be ignored. Striking results emerge from analysing contributions from different regions of the enzymes, by calculating \({\mathrm{\Delta }}C_P^\ddagger\) values for parts of the structures (by recalculating energies and their variances for specific regions only; Fig. 3). Energy contributions from interactions with neighbouring regions are not included, and therefore one should not expect these ‘partial’ \({\mathrm{\Delta }}C_P^\ddagger\) values to add-up to the total value. They do, however, offer new quantitative insights. Conceptually, one may expect differences in partial \({\mathrm{\Delta }}C_P^\ddagger\) values to align with regions that differ in flexibility. This is largely true for some small regions with clear differences in flexibility (e.g. residues 46–70 for KSI; residues 250–321 and 374–459 for MalL; see Fig. 3), but is not obvious throughout the structure, especially for larger regions. Crucially, differences in \({\mathrm{\Delta }}C_P^\ddagger\) are distributed across the full protein structure, whereas significant differences in flexibility are limited to regions that interact with the ligand bound in the active site. This observation bears similarity with the findings of Homans and others regarding entropy differences upon protein–ligand binding: unfavourable entropic contributions (restricted protein dynamics) around the binding site were observed to be (partially) offset by increases in the amplitude of motions in adjacent protein regions15,16.

Fig. 3 Structural fluctuations and partial heat capacity differences between reactant state and transition-state analogue complexes. Top: root-mean square fluctuations from 50–500 ns MD simulations for KSI (left) and MalL (right). Thin lines are individual runs, thick lines the average of ten runs. Residues for which the Cα RMSF difference between states are significant (p < 0.01 as determined by a two-sample t-test) are indicated by grey diamonds (full data in Supplementary Fig. 7). Middle: calculated partial \({\mathrm{\Delta }}C_P^\ddagger\) values for protein regions. Values including contribution from the ligand are indicated (*). Bottom: illustration of KSI (left) and MalL (right) coloured based on partial \({\mathrm{\Delta }}C_P^\ddagger\) regions from the middle pane. Standard deviations are indicated for the total \({\mathrm{\Delta }}C_P^\ddagger\) values; see Supplementary Table 5 for residue ranges and standard deviations for partial \({\mathrm{\Delta }}C_P^\ddagger\) values. Transition-state analogues are shown with space-filling spheres Full size image

KSI, as a dimer, offers the opportunity to assess the dynamical role of the monomer that is distal to substrate turnover. The distal monomer of KSI is the main contributor to reduce \({\mathrm{\Delta }}C_P^\ddagger\) at the TS (Fig. 3). Overall, the catalytic monomer contributes a positive \({\mathrm{\Delta }}C_P^\ddagger\); the N- and C-terminal regions forming the back of the active site and more remote regions contribute a positive \({\mathrm{\Delta }}C_P^\ddagger\), while helix 48–59 that closes over the active site opposite from the catalytic Asp38 rigidifies and contributes a negative \({\mathrm{\Delta }}C_P^\ddagger\). The finding that the non-catalytic chain is as a significant contributor to negative \({\mathrm{\Delta }}C_P^\ddagger\) points to an important role for the oligomer in the temperature dependence of the catalytic process. Enzyme oligomerization is common, indicating an evolutionary advantage25; however, the functional purpose of these quaternary interactions is not well understood. If interactions are optimised to allow global contributions from changes in the distribution of vibrational modes across the multimer26, oligomerization may provide a means to tune the temperature dependence of rates through global contributions to the overall C P change.

The active site of MalL (and TIM barrel enzymes in general27) sits displaced to one side above the TIM barrel core, interacting with a half of the barrel comprising β5–8 and α4–7 (Fig. 3). Analogous to the ligand bound chain of KSI, the catalytic half of the TIM barrel increases in C P at the TS, while more remote protein components, including the second TIM half barrel contribute to the overall negative \({\mathrm{\Delta }}C_P^\ddagger\). The lid domain, consisting of a helix–loop–helix extension above the barrel, contributes significantly to the overall reduction in C P at the TS, consistent with a role in shielding the active site from solvent at the catalytic step. The parallels between the KSI dimer and MalL barrel halves are especially noteworthy in that the TIM barrel is argued to have evolutionary origins as a dimer of (βα) 4 units13,27, the dynamical origins of which may still be discernible in the now fused structure.