Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

# Computer simulations explain the anomalous temperature optimum in a cold-adapted enzyme

## Abstract

Cold-adapted enzymes from psychrophilic species show the general characteristics of being more heat labile, and having a different balance between enthalpic and entropic contributions to free energy barrier of the catalyzed reaction compared to mesophilic orthologs. Among cold-adapted enzymes, there are also examples that show an enigmatic inactivation at higher temperatures before unfolding of the protein occurs. Here, we analyze these phenomena by extensive computer simulations of the catalytic reactions of psychrophilic and mesophilic α-amylases. The calculations yield temperature dependent reaction rates in good agreement with experiment, and also elicit the anomalous rate optimum for the cold-adapted enzyme, which occurs about 15 °C below the melting point. This result allows us to examine the structural basis of thermal inactivation, which turns out to be caused by breaking of a specific enzyme-substrate interaction. This type of behaviour is also likely to be relevant for other enzymes displaying such anomalous temperature optima.

## Introduction

Among the different types of extreme environments found on Earth, cold areas with permanent temperatures below 5 °C cover the major part of the planet, particularly due to the large oceans. Such a cold environment would be detrimental for biochemical processes in species that have the same internal temperature as their surroundings, had evolution not adapted them to this situation. A fundamental question here is how the enzymes of such psychrophilic organisms can overcome the exponential retardation of chemical reaction rates at low temperatures. That is, most mesophilic enzymes lose a large fraction of their maximal activity near the freezing point of water, while cold-adapted enzymes have managed to maintain reaction rates high enough to sustain life in the cold. There is now ample evidence from kinetic studies of numerous orthologous mesophilic–psychrophilic enzyme pairs that cold-adapted enzymes have shifted their thermodynamic activation parameters so that the activation enthalpy (ΔH) is reduced, which is partly counterbalanced by an increased activation entropy penalty (ΔS is more negative)1,2,3,4. This seemingly universal characteristic has the effect that the exponential dampening of rates at lower temperatures is less pronounced, since the reaction rate $$k_{{\mathrm{rxn}}} = C \cdot T{\mathrm{e}}^{ - \Delta G^\ddagger /RT} = C \cdot T{\mathrm{e}}^{ - \Delta H^\ddagger /RT}{\mathrm{e}}^{\Delta S^\ddagger /R}$$ according to standard transition-state theory (where $$\Delta G^\ddagger = \Delta H^\ddagger - T\Delta S^\ddagger$$ is the activation free energy). It is thus now generally accepted that the main adaptive feature of psychrophilic enzymes is the redistribution of the thermodynamic activation parameters compared with mesophilic and thermophilic orthologs, while the activation free energies are usually similar at room temperature1,2,3,4,5.

While the shifted activation enthalpy–entropy balance is what is directly responsible for the higher rates at low temperature, a lower melting temperature is another characteristic of cold-adapted enzymes. Hence, compared with mesophilic enzymes, one typically finds that Tm is shifted downward by about 5–20 °C2,3. It is usually the opposing effects of the exponential rate increase at higher temperatures, and the eventual melting of the enzyme that gives rise to a temperature optimum (Topt) of the catalytic rate. It should be noted here that there is another important difference between psychrophilic enzymes and mesophilic/thermophilic ones, in that the former generally work far away from Topt, while the latter operate closer to the optimum and the melting temperature. This means that the evolutionary pressure on protein stability at the physiological working temperature is bound to be considerably weaker for psychrophilic enzymes, which may be one reason for why their Tm has drifted toward lower values compared with mesophilic orthologs5,6.

The structural origins of both the change in thermodynamic activation parameters and the lower melting temperature have been linked to a higher protein flexibility of cold-adapted enzymes2,4, especially in surface loop regions4,5,7,8,9,10,11,12. In particular, computer simulations that directly evaluated free energy barriers for the catalyzed reactions9,10,11,12,13 have shown that the mobility of active-site residues is generally very similar in psychrophilic and mesophilic enzyme orthologs, while surface loop mobilities may differ considerably and strongly affect the balance between ΔH and ΔS. This finding is entirely in line with the fact that multiple sequence alignments between orthologous psychrophilic and mesophilic enzymes reveal that characteristic mutations are typically located in such loop regions9,10,11.

The most well-studied psychrophilic enzyme is probably the α-amylase from the Antarctic bacterium Pseudoalteromonas haloplanktis (AHA), largely due to the efforts of Gerday, Feller, and coworkers3,14,15,16,17. This multidomain enzyme (Fig. 1) has thus served as a model for cold adaptation, and it has been shown to be a much faster catalyst at low temperatures than its mesophilic and thermophilic orthologs from pig pancreas (PPA) and Bacillus amyloliquefaciens, respectively14. It is also found to be more heat labile with a melting temperature that is ~15 °C lower than the porcine pancreatic variant. Furthermore, it was suggested that the decreased stability of AHA originates from localized increases in flexibility (or local unfolding) of structural regions close to the active site, which may be connected to the higher Km values observed for similar substrates, compared with those in PPA17. It has also been shown that introduction of mutations into AHA, encompassing structural features typical for PPA, diminishes the psychrophilic properties of the former15. In addition, a few computational studies have addressed the catalytic mechanism of α-amylases with QM/MM simulations18,19, which corroborate the standard lysozyme-like mechanism20. Hence, the glycosylation step is predicted to be rate-limiting, where a carboxylate moiety of the enzyme (Asp174 in AHA) attacks the glycosidic bond to yield a covalent enzyme–substrate intermediate, concomitantly with leaving group protonation by a carboxylic acid moiety (Glu200 in AHA). Classical MD simulations comparing AHA and PPA and focusing on the dynamic properties of the apo-enzymes have also been reported21,22.

A particularly interesting feature of the cold-adapted α-amylase is that its rate optimum at 25–30 °C is well below the melting temperature of 44 °C, indicating that inactivation of AHA is due to something other than global unfolding. This is not the case for the mesophilic and thermophilic orthologs where the optimum essentially coincides with Tm14. While a coincidence of optimum rate and the onset of unfolding is also found for most psychrophilic enzymes23,24,25,26,27, there are other such examples where the optimum occurs significantly earlier than melting of the protein28,29. This has usually been interpreted in terms of local unfolding of the enzyme or an unstable active site14,28,29, but a more exotic hypothesis posits that there is a non-negligible heat capacity difference between the transition and ground state ($$\Delta C_P^\ddagger$$)30,31,32. This would give rise to a strong temperature dependence of ΔH and ΔS, leading to curved free energy and activity plots.

It is important here to emphasize that in order to computationally address the outstanding questions of how psychrophilic enzymes can alter the activation enthalpy–entropy balance, and how catalytic rate optima could arise far from the melting temperature, it is absolutely necessary to be able to directly calculate the temperature dependence of the free energy barriers and rates5. Thus, to establish the origins of these effects in AHA, we constructed an empirical valence bond (EVB) model33,34 for the rate-limiting glycosylation step in the α-amylases, based on density functional theory (DFT) quantum mechanical calculations. With this EVB model, we are able to carry out very extensive molecular dynamics (MD) free energy calculations over a wide temperature range, and uncover the origins of these interesting catalytic properties of the psychophilic α-amylase.

## Results

### EVB and DFT modeling of the reaction of α-amylases

The predominant enzymatic activity of α-amylases is the hydrolysis of α-1,4 bonds within starch molecules. The catalytic reaction first involves a glycosylation step resulting in covalent enzyme–substrate intermediate, which is subsequently hydrolyzed by a water molecule in a deglycosylation step, thereby restoring the enzyme-active site35,36. The glycosylation step has been shown to be rate-limiting for α-amylases37, and was described here by the EVB method33,34, using a five-residue glucose oligomer as the substrate in complex with both the psychrophilic P. haloplanktis α-amylase (AHA)36 and the mesophilic porcine pancreatic α-amylase (PPA)38. To calibrate the EVB model, we employed DFT calculations with a continuum solvent model for a small reference system consisting of two carboxylic sidechains, representing the nucleophile and general acid of the reaction (Asp174/197 and Glu200/233 in AHA/PPA, respectively), set to cleave the linkage between two glucose residues. Here, the geometry of reacting moieties represents that found in typical α-amylase active sites (see “Methods”).

The reaction energetics in water for this reference system was obtained by M06-2X/6-311 + G(2d,2p) DFT calculations with the SMD solvent model39,40. Geometry optimization and transition-state (TS) search yielded a concerted oxocarbenium-like transition state with a sole imaginary frequency, where the general acid is essentially deprotonated, and the anomeric carbon distances to the nucleophile and leaving oxygen are 2.31 and 2.33 Å, respectively (Fig. 1c). This type of TS is comparable with that found in several studies reported earlier18,20. The activation and reaction free energy at 25 °C for the glycosylation reaction in water were thus calculated as $$\Delta G_{{\mathrm{ref}}}^\ddagger = 23.9$$ and $$\Delta G_{{\mathrm{ref}}}^0 = - 1.1$$ kcal mol−1, respectively (Supplementary Data 1), notably with a free energy barrier that is about 10 kcal mol−1 higher than that observed for the corresponding enzyme reaction15,16. These results were used to parameterize an (uncatalyzed) EVB reference free energy surface in water, which can then be used to analyze the catalytic effect of the rest of the enzyme, apart from the groups directly involved in the chemistry. Hence, the resulting EVB model exactly reproduces the energetics of the solution reaction, as obtained by the DFT calculations (Fig. 2). The effect of the surrounding enzyme is then obtained by carrying out EVB/MD free energy simulations of the same reaction in the enzyme-active site, where the entire protein is solvated in water. As noted above, it is the efficiency of the EVB representation of the reaction surface that allows for extensive sampling by MD and evaluation of thermodynamic activation parameters. It should also be noted here that the general acid corresponding to Glu200/233 must be treated as protonated in the reactant state since experimental pH-rate profiles unambiguously show that this is the case for the enzyme–substrate complex41,42.

### Origin of the catalytic effect in the α-amylases

The results of the EVB/MD simulations in the psychrophilic (AHA) and mesophilic (PPA) α-amylases are summarized in Fig. 2 in terms of average reaction free energy profiles, each from ~300 independent simulations at 25 °C. It can immediately be seen that the two enzymes lower the activation barrier by ~10 kcal mol−1 compared with the reference reaction in water. Moreover, the psychrophilic enzyme is predicted to be faster than the mesophilic at room temperature, and the calculated free energy barriers are 13.3 ± 0.08 and 15.1 ± 0.14 kcal mol−1 for AHA and PPA, respectively (error bars ± 1 s.e.m. from ~300 replicate simulations). These results agree very well with the available experimental data, which imply free energy barriers of ~14 kcal mol−1 for typical short substrates and an approximately threefold higher kcat value for AHA compared with PPA14,16. The first question that arises from our calculations is where the 10 kcal mol−1 catalytic effect of the α-amylases, compared with the solution reaction, comes from. Actually, a survey of the literature regarding computational modeling of glucosidases indicates that the question of what causes their high catalytic rates basically remains unanswered. That is, computational studies in the field have mainly been focused on details of the catalytic mechanisms (inverting, retaining, concerted, stepwise etc.) and on the resulting free energy barriers, without comparison with any uncatalyzed reference reaction18,19,20,43. A notable exception is the pioneering work of Warshel on lysozyme33,44, who ascribed the catalytic effect mainly to electrostatic stabilization of a proposed ionic oxocarbenium intermediate by the negatively charged nucleophilic carboxylate group. However, such a mechanism can now be ruled out as the covalent enzyme-bound intermediate has been observed experimentally45,46, and reliable QM/MM calculations have provided convincing results in support of the covalent reaction pathway20.

As virtually all computational studies of glucosidase reactions have employed QM/MM methods, where reference reactions in solution are generally avoided, there have been no uncatalyzed reactions to compare with. However, in our case, the results clearly show that just putting the nucleophilic group, general acid, and substrate in the same geometry as in the α-amylases, but surrounded by water, does not at all suffice for achieving the low activation barriers found in the enzymes. It is also clear that while part of the 10 kcal mol−1 barrier reduction originates from stabilization of the covalent intermediate, this does not explain the entire effect (Fig. 2a; Supplementary Fig. 1). Here, analysis of the calculated energetics shows that it is primarily the reorganization free energy33,34,47 that has been significantly reduced in the enzymes compared with the solution reaction (Fig. 2b). This type of transition-state stabilization is one of the major features of enzyme-catalyzed reactions and reflects the preorganization of the active site48, which diminishes the energetic cost of reorienting polar groups surrounding the substrate as its charge distribution changes along the reaction path47,48,49. Such (solvent) reorganization is otherwise a major contributor to the free energy barrier of uncatalyzed reactions in water. The magnitude of the reorganization energy is also directly reflected by the energy gap between the reactant and product states at the reactant minimum. As can be seen from Fig. 2, the reactant minimum in the enzymes is shifted toward smaller absolute values of this energy gap than in the reference reaction, demonstrating that the required reorganization between reactant and products (covalent intermediate) is significantly smaller in the enzymes.

As the negative charge on Asp174/197 in the reactant state ends up on Glu200/233 in the intermediate (Fig. 2a), this migration (delocalization) of charge requires a substantial reorientation of water molecules in the solution reaction, with a large associated free energy cost. In the enzymes, on the other hand, the two carboxylates are largely shielded from solvent, and their negative charge is instead stabilized by interactions with polar groups of the protein and substrate. In particular, the sidechain of Arg172/195 is located between Asp174/197 and Glu200/233, and can thus form ionic interactions with the charged carboxylate, both in the reactant and intermediate states, with only a small reorganizational cost (Fig. 1b).

### Temperature dependence of the AHA and PPA reactions

We calculated the temperature dependence of the glycosylation step in both AHA and PPA by performing ~300 independent EVB/MD calculations of free energy profiles for each enzyme at eight different temperatures between 5 and 40 °C. From these simulations, Arrhenius plots of ΔG/T vs. 1/T were constructed in order to obtain the thermodynamic components of the free energy barrier (Fig. 3). It can immediately be seen from the Arrhenius plots that while PPA essentially shows a linear relationship over the entire temperature range (R2 = 0.94), there is evidently a break in the plot for AHA. For PPA, the linear relation gives values of ΔH = 10.8 and TΔS = −4.3 kcal mol−1 at 15 °C. In the case of AHA, the corresponding values are ΔH = 5.2 and TΔS = −7.8 kcal mol−1, if the entire lower temperature range 5–25 °C is considered, and ΔH = 6.5 and TΔS = −6.6 kcal mol−1 if one restricts the range to 10–25 °C (as done experimentally for AHA16). The predicted values above agree reasonably well with those reported from experiments, where ΔH = 11.1 and TΔS = −2.9 kcal mol−1 for PPA and ΔH = 8.3 and TΔS = −5.1 kcal mol−1 for AHA at 15 °C16. Moreover, it is clear that the computer simulations capture the shift of the activation enthalpy–entropy balance that is observed experimentally both for this and other orthologous mesophilic–psychrophilic enzyme pairs. That is, AHA shows a smaller ΔH and a more negative TΔS than PPA. For comparison, we also calculated the corresponding Arrhenius plot for the uncatalyzed reference reaction, which yields ΔH = 14.8 and TΔS = −8.9 kcal mol−1 at 15 °C (Fig. 3b). This shows that the two enzymes are able to reduce both ΔH and −TΔS, but that the faster psychrophilic enzyme has a larger effect on the enthalpy term. Hence, the transition-state stabilization evidently has both enthalpic and entropic contributions.

It is noteworthy that the break in the Arrhenius plot for AHA occurs precisely where the elusive temperature optimum is observed for the psychrophilic α-amylase (~25 °C) which, as noted above, does not reflect thermal unfolding of the enzyme14. Hence, if we convert our free energy barriers to rate constants (kcat), it can be seen that the calculations predict a rate optimum that coincides with that observed in kinetic experiments (Fig. 3c). This is the first time that such a nontrivial reaction rate behavior has been seen from computer simulations, and this unexpected finding now allows us to explore the cause of the anomalous thermal inactivation in a well-defined enzyme case.

### Origin of the thermal inactivation of psychrophilic α-amylase

It is immediately clear from the present EVB/MD simulations that the cause of the rate decline for AHA above 25 °C does not reflect unfolding of the protein. This is as expected since the timescale of the simulations (sub-μs) would not be able to capture any global melting transition. In order to monitor possible temperature-dependent structural changes near the active site, we also calculated plain MD trajectories in the reactant state at 10, 25, and 40 °C (300 ns at each temperature) for the two enzymes. While the backbone mobility of AHA and PPA around the catalytic residues Asp174/197 and Glu200/233 is found to be very similar (Fig. 4), it turns out that the L7 surface loop following the short helical motif carrying the active-site residues His263/299 and Asp264/300 has significantly higher flexibility in the psychrophilic enzyme, and is, in fact, the most flexible part of AHA (Supplementary Fig. 2). This effect was also seen in earlier MD simulations of the two enzymes without any bound substrate21, and the flexibility of this loop was indeed deemed important in earlier QM/MM simulations19 (where it was denoted L2).

The L7 loop (269–277/305–314) has an alanine insertion in PPA (Ala307), and its conformation becomes partly helical that apparently stabilizes its structure compared with AHA. The strictly conserved Asp264/300 upstream of the L7 loop is a critical residue that anchors the −1 position of the substrate via H bonds to the sugar 2- and 3-OH groups, and it has also been proposed to stabilize the protonated form of the general acid/base Glu200/23337,50. In addition, His263/299 also makes a H bond to the 3-OH group of the substrate, and is likewise strictly conserved. The effect of the increased backbone mobility of the L7 loop propagates upstream to Asp264/300 and His263/299, thereby significantly affecting also their sidechain mobilities, which are found to be higher in AHA than in PPA at all temperatures (Fig. 4). Since the global dynamics of the two enzymes is very similar, apart from the mobility of the L7 loop (Supplementary Fig. 2), it appears that the sequence change in the middle of this loop (AHA: HG−GAGNVIT vs. PPA: HGAGGSSILT) may be the main determinant of its different behavior in the two cases.

Further analysis of the strong Asp264/300 interaction with the substrate reveals that the reactant-state probability distribution of the Asp−COO···O2/O3 distances is considerably wider in AHA and also strongly temperature dependent (Fig. 5a, b). In particular, it can be seen that somewhere above room temperature, there is a notable shift in the probability distribution for AHA toward larger Asp interaction distances, which is not observed in the mesophilic enzyme. This clearly indicates that above 25 °C, where both the calculated and observed rate optima are found, the psychrophilic enzyme starts to populate a less favorable reactant state, leading to a higher free energy barrier for the reaction. In order to rigorously test this hypothesis, we repeated the EVB/MD simulations for AHA at 10, 25, 30, 35, and 40 °C with distance restraints applied to the Asp−COO··· substrate interaction. As intended, these restraints change the probability distributions toward those observed in PPA (Fig. 5b), and are found to abolish the rate decay above 25 °C in the EVB/MD calculations (Fig. 5c). Hence, this computational experiment clearly suggests that it is the different temperature dependence of the substrate interaction with Asp264/300 that causes the rate optimum observed for AHA. Moreover, the Arrhenius plot for AHA with the above restraints now yields a nice linear fit (R2 = 0.98) over the entire temperature range 10–40 °C with computed values of ΔH = 6.8 and TΔS = −6.0 kcal mol−1 at 15 °C (Fig. 5d). These values are thus very similar to those obtained for the unrestrained psychrophilic enzyme below its temperature optimum.

### Kinetic model for the catalytic reaction of AHA

The simplest kinetic model to describe the above situation would be that with a dead-end inhibitory state in equilibrium with the Michaelis complex leading to products

$$\begin{array}{ll}{\mathrm{E}} + {\mathrm{S}} \rightleftharpoons & {\mathrm{ES}}\mathop{\longrightarrow}\limits^{{k_3}}{\mathrm{E}} + P\\ & \upharpoonleft \downharpoonright {\mathrm{K}}_{{\mathrm{eq}}} \\ & {\mathrm{ES}}^{\prime} \\ \end{array}$$
(1)

This two-state model yields $$k_{\mathrm{cat}}=k_3/(1+K_{\mathrm{eq}})$$ and its temperature dependence is thus determined by the thermodynamic parameters $$\Delta H_3^\ddagger ,\Delta S_3^\ddagger ,\Delta H_{{\mathrm{eq}}}$$, and ΔSeq. Fitting the rate curve (Fig. 3c) resulting from our EVB/MD calculations to this model yields $$\Delta H_3^\ddagger = 10.2$$ kcal mol−1, $$\Delta S_3^\ddagger = - 0.00925$$ kcal mol−1 K−1, ΔHeq = 32.0 kcal mol−1, and ΔSeq = +0.10746 kcal mol−1 K−1, and the corresponding curve is shown in Fig. 6a. It is clear that with only eight temperature points, the enthalpy and entropy values resulting from the fit are not very accurate but, nevertheless, give a very informative qualitative picture. That is, the kinetic model requires a strongly temperature-dependent equilibrium with the inactive state in order to capture the rate optimum. Hence, the large enthalpy cost of moving to ES’ is balanced with the entropy gain at precisely around 25 °C, where −TΔSeq = −32 kcal mol−1. At lower temperatures, the system will primarily reside in ES and is higher in ES’. The magnitudes of ΔHeq and ΔSeq appear very reasonable since they primarily would reflect the breaking of ionic H bonds between the substrate and Asp264, which are bound to be strong in terms of enthalpy. There is also always a second solution to Eq. (1) obtained by changing the sign of ΔHeq and ΔSeq, but this can be discarded in the present case as it has a large negative activation enthalpy for the chemical step and the wrong sign of ΔHeq and ΔSeq.

Due to the increasing population of ES’ at higher temperatures, the apparent activation enthalpy $$\Delta H_{{\mathrm{app}}}^\ddagger$$ will become smaller and eventually negative as it becomes dominated by −ΔHeq, the enthalpy associated with going back to ES. At that point, however, the TΔSeq penalty is larger and will dominate the apparent free energy barrier, which is what causes the rate decay in the high-temperature regime. The apparent activation enthalpy and entropy are approximately given by

$$\Delta H_{{\mathrm{app}}}^\ddagger \approx \Delta H_3^\ddagger - \frac{{K_{{\mathrm{eq}}}}}{{1 + K_{{\mathrm{eq}}}}}\Delta H_{{\mathrm{eq}}}$$
(2)

and

$$\Delta S_{{\mathrm{app}}}^\ddagger \approx \Delta S_3^\ddagger - \frac{{K_{{\mathrm{eq}}}}}{{1 + K_{{\mathrm{eq}}}}}\Delta S_{{\mathrm{eq}}}$$
(3)

where the second term in the two equations corresponds to the fractional population of ES’, for which −ΔHeq and −ΔSeq will add to the activation parameters of the chemical step. These expressions give the temperature-dependent apparent activation parameters shown in Fig. 6c. In the low-temperature regime around 15 °C, the model predicts a value of $$\Delta H_{{\mathrm{app}}}^\ddagger$$ around 6 kcal mol−1, in reasonable agreement with our original Arrhenius plot (Fig. 3a). It is important to note, however, that the above two-state kinetic model predicts that the apparent parameters have limiting values of $$\Delta H_{{\mathrm{app}}}^\ddagger = \Delta H_3^\ddagger$$ and $$\Delta S_{{\mathrm{app}}}^\ddagger = \Delta S_3^\ddagger$$ at low temperatures, and $$\Delta H_{{\mathrm{app}}}^\ddagger = \Delta H_3^\ddagger - \Delta H_{{\mathrm{eq}}}$$ and $$\Delta S_{{\mathrm{app}}}^\ddagger = \Delta S_3^\ddagger - \Delta S_{{\mathrm{eq}}}$$ at high temperatures.

An alternative model that has been invoked to explain the occurrence of enzyme temperature optima (that do not reflect unfolding) assumes that there is a constant difference in heat capacity between the reactant and transition states ($$\Delta C_p^\ddagger \ne 0$$)30,31,32. This idea is apparently inspired from studies of protein folding, where it is well-known that there is a considerable change in heat capacity between the folded and unfolded states51. Hence, it is hypothesized that there may also be a non-negligible heat capacity difference between the rate-limiting TS and the Michaelis complex (ES). If so, neither the activation enthalpy nor entropy will be constant, but depend on temperature according to

$$\Delta H^\ddagger \left( T \right) = \Delta H_0^\ddagger + \Delta C_p^\ddagger (T - T_0)$$
(4)
$$\Delta S^\ddagger \left( T \right) = \Delta S_0^\ddagger + \Delta C_p^\ddagger \ln \left( {\frac{T}{{T_0}}} \right)$$
(5)

where the subscript 0 denotes ΔH and ΔS values at a reference temperature T0. Due to the temperature dependence of ΔH and ΔS, this model will yield curved Arrhenius plots for kcat even for a simple one-state model, $${\mathrm{E}} + {\mathrm{S}} \rightleftharpoons {\mathrm{ES}}\mathop { \to }\limits^{k_{{\mathrm{cat}}}} {\mathrm{E}} + {\mathrm{P}}$$, if $$\Delta C_p^\ddagger \ne 0$$. For temperature optima not related to unfolding, a negative value of $$\Delta C_p^\ddagger$$ would be required to fit the rate curve, and this has been interpreted such that the TS is stiffer (lower heat capacity) than the ground state30,31. The only experimental evidence for this type of model, however, comes precisely from curve fitting31,32 since the heat capacity change along an enzymatic reaction coordinate would be very difficult to measure directly.

It is again relatively straightforward to fit our calculated rate vs. temperature curve to the heat capacity model, although the fitted values of $$\Delta H_0^\ddagger$$, $$\Delta S_0^\ddagger$$, and $$\Delta C_p^\ddagger$$ are not very accurate since we only have eight temperature points, as noted above. However, it is an interesting qualitative exercise and the resulting fit is shown in Fig. 6b. This yields $${\Delta} H_{0}^{\ddagger}=-6.2$$ kcal mol−1, $$\Delta S_0^\ddagger$$=−0.06576 kcal mol−1 K−1, and $$\Delta C_p^\ddagger = - 1.13$$ kcal mol−1 K−1, with T0 taken as 298 K. Comparing the above two-state dead-end model with the one-state heat capacity model (Fig. 6a, b), we see that the two theoretical curves are essentially indistinguishable. This demonstrates, beyond any doubt, that curve fitting to these kinetic models cannot be used to determine the actual mechanism behind the curved Arrhenius plot and temperature optimum. However, in our case, we now know from the analysis of conformational populations and our restrained EVB/MD experiments that it is the two-state model that is the correct one, and we have also been able to identify the nature of the ES and ES’ states. It is further interesting to note that the one-state heat capacity model predicts a very strange behavior of ΔH(T) and ΔS(T) as shown in Fig. 6d. The activation enthalpy is thus predicted to grow linearly at low temperatures toward very high values and, conversely, become large and negative at high temperatures, while the −TΔS term behaves in an opposite manner. This is quite different from the prediction of the two-state model where $$\Delta H_{{\mathrm{app}}}^\ddagger$$ and $$\Delta S_{{\mathrm{app}}}^\ddagger$$ reach limiting values at low and high values of T. Hence, the physical basis for why the activation enthalpy and entropy of an enzyme reaction would behave as predicted by the heat capacity model (Fig. 6d) is rather obscure.

## Discussion

The increasing interest in enzymes from extremophilic species and, in particular those from psychrophiles, has raised a number of fundamental questions regarding how evolution has shaped the temperature dependence of their catalytic activity and thermal stability. A major discovery in the field was the observation by Somero and coworkers1 that the thermodynamic activation parameters of cold-adapted enzymes have shifted toward a lower value of ΔH and a more negative ΔS, compared with mesophilic orthologs. This finding has since then been confirmed for many different classes of enzymes2,5, and appears to be a universal feature of cold adaptation. The other universal feature is the reduced thermal stability that psychrophilic enzymes display2. The α-amylase from P. haloplanktis (AHA) obeys both these rules of psychrophilicity, and also shows the enigmatic behavior of becoming inactivated before protein melting sets in, as the temperature is increased14.

The present computer simulations reveal several hitherto hidden properties of AHA, which explain its catalytic behavior. First, it is evident from our comparison to an uncatalyzed reference reaction in solution, containing exactly the same catalytic groups as the enzyme, that a major contribution to the low activation barrier in the α-amylases originates from reduction of the reorganization free energy of the chemical reaction. This effect is often simply referred to as active-site preorganization, but can, however, be quantified, and has been seen for many other types of enzymes as well33,47,49. In fact, it may well be a characteristic of the glucosidases utilizing two enzyme carboxylic sidechains as their nucleophile and general acid/base for cleaving glycosidic bonds. In these reactions, there is a migration of a net negative charge over a distance of about 5 Å, which is associated with a large reorganization free energy penalty in solution. Hence, if the enzyme can provide stabilizing interaction for the negative charge at both of its carboxylate locations during the reaction, the reorganization energy can be reduced significantly. It appears that the reason why this major catalytic effect has not been revealed earlier for glucosidases, is that no comparisons to uncatalyzed water reactions have been reported in previous computational studies.

The second major result from our EVB/MD simulations at different temperatures is that, not only the free energy barriers, but also the experimental values16 of ΔH and TΔS reported at 15 °C for AHA and PPA are reproduced remarkably well, with shifts of several kcal mol−1 in the expected direction between the two enzymes. While the 3D structures of AHA and PPA are very similar, with the same overall fold (Fig. 1a), their sequence identity is, however, only 47% with several insertions/deletions. The fact that these enzymes are also rather large, with almost 500 amino acids, thus makes it difficult to identify specific protein regions responsible for the cold adaptation of AHA. It should also be noted that the shifts in ΔH and TΔS are smaller for the α-amylases than, e.g., endonucleases27, xylanases29, and trypsins52, where they are on the order of 10 kcal mol−1. Experimental work on AHA mutants has, however, shown that a combination of five or six AHA → PPA substitutions are sufficient for rendering mesophilic characteristics in terms of kcat, Km, ΔH, and TΔS to the psychrophilic enzyme15. Neither of the corresponding single mutations, on the other hand, were able to cause these changes, demonstrating that multiple contributions are at play. Interestingly, all but one (K300R) of the selected mutations were located on the protein surface, 20–25 Å away from the reaction center, which reinforces our earlier findings regarding the importance of the protein surface in cold adaptation9,10,11,12. However, since it is probably easier to destroy psychrophilic properties than to create them, a more critical test would be to instead identify minimal sets of mutations that confer cold adaptation to the mesophilic enzyme.

The most noteworthy result from the present computer simulations is the break in the Arrhenius plot for AHA and the corresponding calculated rate maximum at around 25 °C. This type of optimum has never been predicted by computational enzyme studies before. The optimum also occurs precisely where it is found experimentally, which is about 20 °C below the melting temperature of the enzyme14. We were thus able here to analyze the origin of such anomalous temperature optima in the case of AHA. Our results unambiguously show that it is primarily the Asp264 interaction with the substrate that breaks above room temperature and causes the rate decline. The evidence comes from analysis of backbone and sidechain atomic mobilities, conformational populations, and computational experiments where the interaction is maintained at higher temperatures through restraints. Hence, the kinetics of the psychrophilic enzyme is best described by a two-state model in which ES is in equilibrium with a dead-end state ES’. Such a macroscopic model yields very reasonable values for the apparent activation enthalpy and entropy as a function of temperature, and predicts that the ES  ES′ equilibrium is characterized by a large enthalpy (~30 kcal mol−1) that is counterbalanced by an equally large entropy term at room temperature, indicative of breaking/forming a strong interaction. The two-state dead-end model is thus similar to that proposed by Danson for explaining unusual curved Arrhenius plots in enzyme kinetics53. That model also adds a slow irreversible inactivation of the enzyme from the ES’ state, and describes how the rate vs. temperature curve changes with time.

We further showed here that the heat capacity model30,31,32 can fit our calculated rate curve equally well as the two-state model, although it is clearly incorrect in this case. This strongly suggests that curve fitting alone cannot be used to prove the somewhat controversial idea that $$\Delta C_p^\ddagger$$ could be nonzero for an elementary enzyme reaction step, and that this would be the underlying reason for curved Arrhenius plots and unusual temperature optima of the type discussed here. Moreover, we show that the two alternative models make very different predictions regarding the temperature dependence of the apparent activation enthalpies and entropies. While $$\Delta H_{{\mathrm{app}}}^\ddagger$$ and $$\Delta S_{{\mathrm{app}}}^\ddagger$$ in the two-state model reach limiting values at low and high temperatures, determined by the chemical step and inactivation equilibrium, the heat capacity model predicts a seemingly unphysical behavior of these quantities. Hence, $$\Delta H_{{\mathrm{app}}}^\ddagger$$ is predicted to be inversely proportional to T over the entire temperature range, and reaches a value of about 330 kcal mol−1 at 0 K for the AHA reaction. What the microscopic origin of such an enormous energy barrier would be is difficult to understand. Likewise, a large negative activation enthalpy is predicted at high temperatures for an elementary chemical reaction step, for which the physical explanation is also unclear. It should, however, be noted that the equilibrium model of Eq. (1) predicts an increased apparent heat capacity of the combined reactant state when both ES and ES’ are significantly populated, since the temperature derivative of the second term of Eq. (2) is not zero. This causes a negative dip of the apparent $$\Delta C_p^\ddagger$$ in the temperature region where the population transition occurs, while it goes to zero at low and high temperatures (Supplementary Fig. 3). However, in this case, the negative dip in $$\Delta C_p^\ddagger$$ is a consequence of the ES  ES′ equilibrium, and $$\Delta C_p^\ddagger$$ is not a temperature-independent constant.

Finally, it is interesting to note that, in contrast to earlier computational studies of psychrophilic enzymes, the α-amylases show distinct differences in the mobilities of the conserved active-site residues Asp264/300 and His263/299. That is, earlier work on a number of different enzymes has revealed that mutations and flexibility differences are generally located at surface loops of the protein, while the conserved active-site residues showed the same low mobility in both psychrophilic and mesophilic orthologs. Moreover, increased surface loop mobility in psychrophilic enzymes could be directly related to the changes in activation enthalpies and entropies9,10,11,12. The influence of protein surface residues far away from the active site on thermodynamic activation parameters is further illustrated by the multiple AHA mutations analyzed by Feller and coworkers15. In AHA and PPA, the catalytic groups and core of the active site do indeed show very similar and low mobility (Fig. 4), but Asp264/300 and His263/299 are clearly more flexible in the psychrophilic enzyme although the residues are totally conserved among the α-amylases54. The reason for this is to be found in the L7 loop downstream of the two residues, which is the most mobile part of AHA and much more flexible than in PPA. Hence, it again appears that it is mutation in a surface loop that is ultimately responsible for the different behavior of the two enzymes with regard to their temperature optima.

## Methods

### MD simulations

Models for the psychrophilic P. haloplanktis (AHA) and mesophilic porcine (Sus scrofa) pancreatic α-amylase (PPA) were based on structural information accessible in the Protein Data Bank (PDB). The highest X-ray resolution structure available for each enzyme complexed with an acarbose inhibitor was chosen, 1G9436 (1.74-Å resolution) for AHA and 1HX038 (1.38 Å) for PPA. Acarbose was modified to represent a five-residue glucose oligomer, linked by α1–4 glycosidic bonds. Crystallographic waters were retained, except for those clashing with atoms of the modified substrate. The AHA and PPA systems were prepared for EVB/MD simulations33,34 by creating the enzyme–substrate complexes based on the crystallographic structures and solvating them with a 45-Å-radius spherical droplet of TIP3P water55 with origin at the center of mass of each complex, thereby covering the entire protein. Water molecules close to the spherical boundary were subjected to radial and polarization restraints according to the SCAAS model56,57. To construct a reference system for calibrating the EVB potential surface, a simplified model was used based on the conserved α-amylase active-site structure. The model comprised the sidechains of the nucleophilic group (Asp174/197) and the protonated general acid (Glu200/233), both truncated at the Cβ carbon, together with a maltose disaccharide molecule representing the substrate at the −1 and +1 positions. The reference system was solvated in a 18-Å-radius sphere of TIP3P water. The OPLS-AA/M force field58 was used in all calculations. Nonbonded interactions beyond a 10-Å cutoff were treated by local reaction field multipole expansion method59, except for the reacting groups for which all interactions were explicitly calculated. All EVB/MD calculations were carried out with the Q program57, using a 1-fs time step as described earlier9.

To calculate average backbone and sidechain-positional atomic fluctuations (RMSF), the enzyme–substrate complexes were solvated by a cubic TIP3P water box with side length 100 Å using periodic boundary conditions (PBC) in GROMACS60. Here, we carried out longer MD simulations of the reactant state in AHA and PPA, again using the OPLS-AA/M force field. The standard GROMACS PBC simulations involved 6 replicas of 50 ns each for both AHA and PPA at 283, 298, and 313 K and used a 2-fs time step.

### DFT cluster model of reference reaction

The DFT cluster model was used to determine the structure and energetics of the rate-limiting glycosylation transition state in a continuum solvent model. The model consisted of the same chemical groups as those used in the EVB reference system, except that the amino acid sidechains were truncated at the Cα carbon. The glucose ring at the −1 position was proposed to be in half-chair conformation in the initial TS structure guess (see Fig. 2c for the resulting optimized geometry). The energetics of this cluster model was evaluated with the M06-2X functional39 using the Gaussian09 program61. The transition state was first optimized, and then geometry optimization was carried out toward the reactant and intermediate states, using the 6–31 G(d,p) basis set. The α-carbons of the two amino acid sidechains were kept fixed during the optimizations. Electronic energies were then recalculated for the optimized geometries with the larger 6–311+G(2d,2p) basis set, with solvent effects from surrounding water estimated using SMD model40. Zero-point energies were obtained from frequencies calculated at the same level of theory as geometry optimization. Free energy estimates were thus obtained by adding the zero-point energy and solvent contributions to the electronic energies obtained by the large basis set.

### EVB/MD calculations of free energy profiles

The systems were first equilibrated for 246 ps at increasing temperatures, starting from 1 K, and heated in a stepwise manner toward the final desired temperature for EVB/MD production runs. EVB/MD free energy perturbation simulations11,34 were carried out with 51 discrete λ windows, starting at λ1=λ2 = 0.5, near the transition state, and propagated toward reactants (λ1 = 1) and intermediate (λ2 = 1). To calibrate the EVB reaction potential energy surface, the average free energy profile from 250 independent simulations for the reference reaction system in water was fitted to the corresponding free energy profile obtained by the DFT cluster calculations at 298 K. The temperature dependence of the reference reaction was analyzed, after calibration, by repeating the free energy profile calculations also at 283, 290.5, 305.5, and 313 K, again using 250 replicas at each temperature. Each individual free energy profile encompassed 510 ps, and the total simulation time for obtaining the reference system data was thus 0.64 μs.

The enzyme EVB/MD simulations for AHA and PPA were run at eight different temperatures from 278 to 313 K, with 5-K intervals. At each temperature, about 300 randomized individual replicate simulations were carried out to obtain the free energy profiles, using the same free energy perturbation scheme11 as above. In addition, a system with added distance restraints was analyzed for AHA. Weak harmonic restraints between the His263 and Asp264 sidechains of the enzyme and the O2 and O3 hydroxyl groups of the substrate glucose in the –1 position were then imposed during the EVB/MD calculations. This restrained system was simulated at 283, 298, 303, 308, and 313 K, again with 300 randomized replicas per temperature. The total simulation time of all enzyme EVB/MD calculations amounted to 3.2 μs. Fitting of enzyme activity curves to kinetic models was done with Gnuplot (http://www.gnuplot.info).

### Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

## Data availability

The data that support this study are available from the authors upon reasonable request.

## Code availability

Source code for the Q program57 is available at https://github.com/qusers/Q6, or upon request from J.Å.

## References

1. 1.

Low, P. S., Bada, J. L. & Somero, G. N. Temperature adaptation of enzymes: roles of the free energy, the enthalpy, and the entropy of activation. Proc. Natl Acad. Sci. USA 70, 430–432 (1973).

2. 2.

Siddiqui, K. S. & Cavicchioli, R. Cold-adapted enzymes. Annu. Rev. Biochem. 75, 403–433 (2006).

3. 3.

Feller, G. & Gerday, C. Psychrophilic enzymes: hot topics in cold adaptation. Nat. Rev. Microbiol. 1, 200–208 (2003).

4. 4.

Fields, P. A. & Somero, G. N. Hot spots in cold adaptation: localized increases in conformational flexibility in lactate dehydrogenase A4 orthologs of Antarctic notothenioid fishes. Physiology 95, 11476–11481 (1998).

5. 5.

Åqvist, J., Isaksen, G. V. & Brandsdal, B. O. Computation of enzyme cold adaptation. Nat. Rev. Chem. 1, 0051 (2017).

6. 6.

Harms, M. J. & Thornton, J. W. Evolutionary biochemistry: revealing the historical and physical causes of protein properties. Nat. Rev. Genet. 14, 559–571 (2013).

7. 7.

Wallon, G. et al. Sequence and homology model of 3-isopropylmalate dehydrogenase from the psychrotrophic bacterium Vibrio sp. I5 suggest reasons for thermal instability. Protein Eng. Des. Sel. 10, 665–672 (1997).

8. 8.

Alvarez, M. et al. Triose-phosphate isomerase (TIM) of the psychrophilic bacterium Vibrio marinus: kinetic and structural properties. J. Biol. Chem. 273, 2199–2206 (1998).

9. 9.

Isaksen, G. V., Åqvist, J. & Brandsdal, B. O. Protein surface softness is the origin of enzyme cold-adaptation of trypsin. PLoS Comput. Biol. 10, e1003813 (2014).

10. 10.

Åqvist, J. Cold adaptation of triosephosphate isomerase. Biochemistry 56, 4169–4176 (2017).

11. 11.

Sočan, J., Kazemi, M., Isaksen, G. V., Brandsdal, B. O. & Åqvist, J. Catalytic adaptation of psychrophilic elastase. Biochemistry 57, 2984–2993 (2018).

12. 12.

Sočan, J., Isaksen, G. V., Brandsdal, B. O. & Åqvist, J. Towards rational computational engineering of psychrophilic enzymes. Sci. Rep. 9, 19147 (2019).

13. 13.

Bjelić, S., Brandsdal, B. O. & Åqvist, J. Cold adaptation of enzyme reaction rates. Biochemistry 47, 10049–10057 (2008).

14. 14.

D’Amico, S., Marx, J.-C., Gerday, C. & Feller, G. Activity-stability relationships in extremophilic enzymes. J. Biol. Chem. 278, 7891–7896 (2003).

15. 15.

D’Amico, S., Gerday, C. & Feller, G. Temperature adaptation of proteins: engineering mesophilic-like activity and stability in a cold-adapted α-amylase. J. Mol. Biol. 332, 981–988 (2003).

16. 16.

D’Amico, S., Gerday, C. & Feller, G. Dual effects of an extra disulfide bond on the activity and stability of a cold-adapted α-amylase. J. Biol. Chem. 277, 46110–46115 (2002).

17. 17.

D’Amico, S., Sohier, J. S. & Feller, G. Kinetics and energetics of ligand binding determined by microcalorimetry: insights into active site mobility in a psychrophilic α-amylase. J. Mol. Biol. 358, 1296–1304 (2006).

18. 18.

Pinto, G. P. et al. establishing the catalytic mechanism of human pancreatic α-amylase with QM/MM methods. J. Chem. Theory Comput. 11, 2508–2516 (2015).

19. 19.

Kosugi, T. & Hayashi, S. Crucial role of protein flexibility in formation of a stable reaction transition state in an α-amylase catalysis. J. Am. Chem. Soc. 134, 7045–7055 (2012).

20. 20.

Bowman, A. L., Grant, I. M. & Mulholland, A. J. QM/MM simulations predict a covalent intermediate in the hen egg white lysozyme reaction with its natural substrate. Chem. Commun. 37, 4425–4427 (2008).

21. 21.

Pasi, M., Riccardi, L., Fantucci, P., De Gioia, L. & Papaleo, E. Dynamic properties of a psychrophilic α-amylase in comparison with a mesophilic homologue. J. Phys. Chem. B 113, 13585–13595 (2009).

22. 22.

Papaleo, E., Pasi, M., Tiberti, M. & De Gioia, L. Molecular dynamics of mesophilic-like mutants of a cold-adapted enzyme: insights into distal effects induced by the mutations. PLoS ONE 6, e24214 (2011).

23. 23.

Lanes, O., Leiros, I., Smalås, A. & Willassen, N. Identification, cloning, and expression of uracil-DNA glycosylase from Atlantic cod (Gadus morhua): characterization and homology modeling of the cold-active catalytic domain. Extremophiles 6, 73–86 (2002).

24. 24.

Assefa, N. G., Niiranen, L., Willassen, N. P., Smalås, A. & Moe, E. Thermal unfolding studies of cold adapted uracil-DNA N-glycosylase (UNG) from Atlantic cod (Gadus morhua). A comparative study with human UNG. Comp. Biochem. Physiol. B Biochem. Mol. Biol. 161, 60–68 (2012).

25. 25.

Xu, Y., Feller, G., Gerday, C. & Glansdorff, N. Moritella cold-active dihydrofolate reductase: are there natural limits to optimization of catalytic efficiency at low temperature? J. Bacteriol. 185, 5519–5526 (2003).

26. 26.

Coquelle, N., Fioravanti, E., Weik, M., Vellieux, F. & Madern, D. Activity, stability and structural studies of lactate dehydrogenases adapted to extreme thermal environments. J. Mol. Biol. 374, 547–562 (2007).

27. 27.

Altermark, B., Niiranen, L., Willassen, N. P., Smalås, A. O. & Moe, E. Comparative studies of endonuclease I from cold-adapted Vibrio salmonicida and mesophilic Vibrio cholerae. FEBS J. 274, 252–63 (2007).

28. 28.

Georlette, D. et al. Structural and functional adaptations to extreme temperatures in psychrophilic, mesophilic, and thermophilic DNA ligases. J. Biol. Chem. 278, 37015–37023 (2003).

29. 29.

Collins, T., Meuwis, M.-A., Gerday, C. & Feller, G. Activity, stability and flexibility in glycosidases adapted to extreme thermal environments. J. Mol. Biol. 328, 419–428 (2003).

30. 30.

Hobbs, J. K. et al. Change in heat capacity for enzyme catalysis determines temperature dependence of enzyme catalyzed rates. ACS Chem. Biol. 8, 2388–2393 (2013).

31. 31.

van der Kamp, M. W. et al. dynamical origins of heat capacity changes in enzyme-catalysed reactions. Nat. Commun. 9, 1177 (2018).

32. 32.

Nguyen, V. et al. Evolutionary drivers of thermoadaptation in enzyme catalysis. Science 355, 289–294 (2017).

33. 33.

Warshel, A. Computer Modeling of Chemical Reactions in Enzymes and Solutions (John Wiley and Sons, New York, 1991).

34. 34.

Åqvist, J. & Warshel, A. Simulation of enzyme reactions using valence bond force fields and other hybrid quantum/classical approaches. Chem. Rev. 93, 2523–2544 (1993).

35. 35.

Rydberg, E. H. et al. Mechanistic analyses of catalysis in human pancreatic α-amylase: detailed kinetic and structural studies of mutants of three conserved carboxylic acids. Biochemistry 41, 4492–4502 (2002).

36. 36.

Aghajari, N., Roth, M. & Haser, R. Crystallographic evidence of a transglycosylation reaction: ternary complexes of a psychrophilic α-amylase. Biochemistry 41, 4273–4280 (2002).

37. 37.

Brayer, G. D. et al. Subsite mapping of the human pancreatic α-amylase active site through structural, kinetic, and mutagenesis techniques. Biochemistry 39, 4778–4791 (2000).

38. 38.

Qian, M. et al. Enzyme-catalyzed condensation reaction in a mammalian α-amylase. high-resolution structural analysis of an enzyme−inhibitor complex. Biochemistry 40, 7700–7709 (2001).

39. 39.

Zhao, Y. & Truhlar, D. G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four m06-class functionals and 12 other functionals. Theor. Chem. Acc. 120, 215–241 (2008).

40. 40.

Marenich, A. V., Cramer, C. J. & Truhlar, D. G. Universal solvation model based on solute electron density and on a continuum model of the solvent defined by the bulk dielectric constant and atomic surface tensions. J. Phys. Chem. B 113, 6378–6396 (2009).

41. 41.

Ishikawa, K., Matsui, I., Honda, K. & Nakatani, H. Substrate-dependent shift of optimum pH in porcine pancreatic α-amylase-catalyzed reactions. Biochemistry 29, 7119–7123 (1990).

42. 42.

Feller, G., le Bussy, O., Houssier, C. & Gerday, C. Structural and functional aspects of chloride binding to Alteromonas haloplanctis α-amylase. J. Biol. Chem. 271, 23836–23841 (1996).

43. 43.

Liu, J., Wang, X. & Xu, D. QM/MM study on the catalytic mechanism of cellulose hydrolysis catalyzed by cellulase Cel5A from Acidothermus cellolyticus. J. Phys. Chem. B 114, 1462–1470 (2010).

44. 44.

Warshel, A. Calculations of enzymatic reactions: calculations of pKa, proton transfer reactions, and general acid catalysis reactions in enzymes. Biochemistry 20, 3167–3177 (1981).

45. 45.

Vocadlo, D. J., Davies, G. J., Laine, R. & Withers, S. G. Catalysis by hen egg-white lysozyme proceeds via a covalent intermediate. Nature 412, 835–838 (2001).

46. 46.

Zhang, R. et al. Directed “in situ” inhibitor elongation as a strategy to structurally characterize the covalent glycosyl-enzyme intermediate of human pancreatic α-amylase. Biochemistry 48, 10752–10764 (2009).

47. 47.

Warshel, A., Hwang, J. K. & Åqvist, J. Computer simulations of enzymatic reactions: examination of linear free-energy relationships and quantum-mechanical corrections in the initial proton-transfer step of carbonic anhydrase. Faraday Discuss. 93, 225–238 (1992).

48. 48.

Warshel, A. Energetics of enzyme catalysis. Proc. Natl Acad. Sci. USA 75, 5250–5254 (1978).

49. 49.

Feierberg, I. & Åqvist, J. Computational modeling of enzymatic keto-enol isomerization reactions. Theor. Chem. Acc. 108, 71–84 (2002).

50. 50.

Aghajari, N., Feller, G., Gerday, C. & Haser, R. structural basis of α-amylase activation by chloride. Protein Sci. 11, 1435–1441 (2002).

51. 51.

Prabhu, N. V. & Sharp, K. A. Heat capacity in proteins. Annu. Rev. Phys. Chem. 56, 521–548 (2005).

52. 52.

Lonhienne, T., Gerday, C. & Feller, G. Psychrophilic enzymes: revisiting the thermodynamic parameters of activation may explain local flexibility. Biochim. Biophys. Acta 1543, 1–10 (2000).

53. 53.

Daniel, R. M. & Danson, M. J. A new understanding of how temperature affects the catalytic activity of enzymes. Trends Biochem. Sci. 35, 584–591 (2010).

54. 54.

D’Amico, S., Gerday, C. & Feller, G. Structural similarities and evolutionary relationships in chloride-dependent α-amylases. Gene 253, 95–105 (2000).

55. 55.

Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W. & Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79, 926–935 (1983).

56. 56.

King, G. & Warshel, A. A surface constrained all‐atom solvent model for effective simulations of polar solutions. J. Chem. Phys. 91, 3647–3661 (1989).

57. 57.

Marelius, J., Kolmodin, K., Feierberg, I. & Åqvist, J. Q: a molecular dynamics program for free energy calculations and empirical valence bond simulations in biomolecular systems. J. Mol. Graph. Model. 16, 213–225 (1998).

58. 58.

Robertson, M. J., Tirado-Rives, J. & Jorgensen, W. L. Improved peptide and protein torsional energetics with the OPLS-AA force field. J. Chem. Theory Comput. 11, 3499–3509 (2015).

59. 59.

Lee, F. S. & Warshel, A. A local reaction field method for fast evaluation of long‐range electrostatic interactions in molecular simulations. J. Chem. Phys. 97, 3100–3107 (1992).

60. 60.

Abraham, M. J. et al. GROMACS: high performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 1–2, 19–25 (2015).

61. 61.

Frisch, M. J. et al. Gaussian 09, Revision D.01. (Gaussian Inc., Wallingford, CT, 2009).

## Acknowledgements

Support from the Swedish Research Council (VR), the Knut and Alice Wallenberg Foundation, and the Swedish National Infrastructure for Computing (SNIC) is gratefully acknowledged. Open access funding provided by Uppsala University.

## Author information

Authors

### Contributions

J.S. performed experiments. J.Å. designed the study. J.S., M.P., and J.Å. analyzed the data and wrote the paper.

### Corresponding author

Correspondence to Johan Åqvist.

## Ethics declarations

### Competing interests

The authors declare no competing interests.

Peer review information Nature Communications thanks Guillaume Lamoureux, and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Rights and permissions

Reprints and Permissions

Sočan, J., Purg, M. & Åqvist, J. Computer simulations explain the anomalous temperature optimum in a cold-adapted enzyme. Nat Commun 11, 2644 (2020). https://doi.org/10.1038/s41467-020-16341-2

• Accepted:

• Published:

• ### Hexamerization and thermostability emerged very early during geranylgeranylglyceryl phosphate synthase evolution

• Cosimo Kropp
• , Kristina Straub
• , Mona Linde
•  & Patrick Babinger

Protein Science (2021)

• ### Cold-Active β-Galactosidases: Insight into Cold Adaptation Mechanisms and Biotechnological Exploitation

• Marco Mangiagalli
•  & Marina Lotti

Marine Drugs (2021)

• ### Biomolecular QM/MM Simulations: What Are Some of the “Burning Issues”?

• Qiang Cui
• , Tanmoy Pal
•  & Luke Xie

The Journal of Physical Chemistry B (2021)

• ### A meta-analysis of the activity, stability, and mutational characteristics of temperature-adapted enzymes

• Stewart Gault
• , Peter M. Higgins
• , Charles S. Cockell
•  & Kaitlyn Gillies

Bioscience Reports (2021)

• ### Hidden Conformational States and Strange Temperature Optima in Enzyme Catalysis

• Johan Åqvist
• , Jaka Sočan
•  & Miha Purg

Biochemistry (2020)