# On the Physical and Evolutionary Limits to the Rates of Enzyme-Catalyzed Reactions

Authors: Carlo Canepa
Journal: Journal of Chemistry
 (2020)

---

## Abstract

An existing model for the rate coefficients of enzyme-catalyzed processes involves the regularized gamma function of Euler replacing the exponential dependence of the rate coefficient from the reaction barrier. The application of this model to experimental data, on one hand, validates the model by correctly describing the negative curvature of Eyring plots. On the other hand, this analysis evidences that enzymes never reach the maximum theoretical efficiency, a counterintuitive fact that requires an explanation. This work interprets this evolutionary limit in terms of the necessity of living systems to achieve and maintain homeostasis. Further validation of the expression for the rate coefficients comes from the analysis of the discrepancy between the theoretically predicted energy difference between reactants and products in a chemical equilibrium and the corresponding value obtained by regression to the classical expression for the equilibrium constant. The discrepancy is resolved by making use of the proposed model.


---

## Body

## 1. Introduction

In previous works, a functional form of the rate coefficient for enzyme-catalyzed reactions based on the coupling of the protein vibrational modes to the reaction coordinate has been proposed [1, 2]. The model was tested on a series of both catalyzed and uncatalyzed reactions, for which the temperature-dependent rate coefficients have been measured [3]. The outcome was that, in the majority of cases, the curvature of the Arrhenius plots was found to be negative, a feature that transition state theory (TST) in its various forms is not able to explain. On the other hand, the curvature computed according to the model is always negative, an indication that the new model is successful. We first give an outline of the model to make clear what the topic of this investigation is.The thermally averaged rate coefficient of a chemical reaction is expressed as the integral of the product of an energy-dependent rate coefficientkrε and the classical Boltzmann distribution of low-frequency vibrational modes:(1)krT=∫uTS∞krερεzRe−ε/kTdε.The lower limit of the integral is the threshold valueuTS of the potential energy that the substrate must exceed to be able to cross the reaction barrier. The functional form of the energy-dependent krε is expressed as the ratio of the number of states of the transition structure to the number of states of the reactant, multiplied by the intrinsic frequency for the energy transfer from active modes to the reaction coordinate (ν):(2)krε=νρTSεdερεdε.The total number of modess includes the subset of all the excited modes that can transfer energy to the reaction coordinate (indicated by a), implying ℏωi/kT≪1, and the classical expression for the vibrational density of states(3)ρ=εs−1s−1!∏i=1sℏωi,leading to the partition function(4)zv=∏i=1skTℏωi.The density of states of the transition structure at the total energyε (ρTS) is expressed as the convolution integral of the classical approximation for the vibrational density of states of the subsystem of a oscillators and the subsystem of the remaining s−a oscillators:(5)ρTS=∫0εε¯a−1a−1!∏i=1aℏωiε−ε¯s−a−1s−a−1!∏i=a+1sℏωidε¯=zTSkTxs−1Γs∫ζ/x1ua−1us−a−1du∫01ua−1us−a−1du=zTSkTxs−1ΓsBa,s−a,ζ/xBa,s−a,in the variables u=ε¯/ε and x=ε/kT, where the last equality defines the regularized beta function Ba,s−a,ζ/x/Ba,s−a. Taking the anharmonicity of the low-frequency vibrational modes into account, we would have an increase in the vibrational density of states of the transition structure with respect to equation (5) [4]. However, this effect is expected to be small, since the low-frequency modes would be far from the dissociation limit, having energy ≈kT. The thermally averaged rate coefficient in equation (1) thus becomes(6)krT=νzR∫uTS∞ρTSεe−ε/kTdε=νzTSzR1Γs∫ζ∞Ba,s−a,ζ/xBa,s−axs−1e−xdx,with ζ=uTS/kT. Making use of the theorem [5],(7)1Γs∫ζ∞Ba,s−a,ζ/xBa,s−axs−1e−xdx=1Γa∫ζ∞xa−1e−xdx,we obtain the expression, independent of s, for the rate coefficient in terms of the regularized gamma function Q:(8)kr=νzTSzR1Γa∫ζ∞xa−1e−xdx=νzTSzRΓa,ζΓa=νzTSzRQa,ζ.This expression is plotted in Figure1, and it exhibits an inflection point at ζ=a−1. The above expression for kr is able to account for the large average mass of catalytic proteins, pointed out by Pawlowski and Zielenkiewicz [6] in a seminal paper on the correlation between the molecular mass and the activation energy (according to the Eyring equation) of enzymes. Analyzing 6,915 cases, the most common value for the molecular mass of enzymes was found to be around 75kDa. In expression (8), the value of active mode a is a subset of the total number of modes s, and in order for a to reach the values necessary for enzyme catalysis, s needs to be large, hence the large mass.Figure 1
Plots ofQa,ζ versus a at various fixed values of ζ.For integer values ofa, we can integrate Q analytically, obtaining(9)kr=νzTSzRe−ζ∑i=0a−1ζii!.At the inflection point, the value of the functionQ is given by(10)Qa,a−1=e−a−1∑j=0a−1a−1jj!,and its values are plotted in Figure 2 as a function of a.Figure 2
Values ofQa,a−1 versus a.Figures3–5 report logarithmic plots of experimental data along with the best-fit curves for the rate coefficients of three reactions (one in solvent and two catalyzed by enzymes). For a detailed analysis of these and other cases [3], Figures 3–5 show that equation (8) gives a better fit to the experimental data with respect to the TST expression and, most importantly, reproduces the correct negative curvature of the plots. How anharmonic potentials affect the frequency ν is probably specific to each system, i.e., the details of the constraints on the system of coupled oscillators and its structure. An example of the behavior of an harmonic oscillator with structural constraints coupled to a second oscillator through an anharmonic potential is discussed in [7]. The analysis of experimental data for these and other cases of enzyme-catalyzed reactions indicates that a/ζ≃0.9, and the most common values of Qa,ζ are slightly less than 1/2, since a typical value of the reduced barrier ζ usually exceeds 40. This prompts the question of why the evolutionary pressure did not push the structure of enzymes to reach the maximum efficiency at Qa,ζ≃1.Figure 3
Eyring plots [3] according to TST (green) of rate coefficients for the neutral hydrolysis of methyl trifluoroacetate (MeTFA) in H2O/DMSO. The corresponding curve from equation (8), with a=35,νzTS/zR=0.50s−1,ζ=49.29,a/ζ=0.71, is shown in blue and the experimental data in red.
(a)(b)Figure 4
Eyring plots [3] according to TST (green) of rate coefficients for the hydrolysis of Suc-Ala-PhepNA catalyzed by α-chymotrypsin. The corresponding curve from equation (8), with a=133,νzTS/zR=1.619s−1,ζ=156.5,a/ζ=0.85, is shown in blue and the experimental data in red.
(a)(b)Figure 5
Eyring plots [3] according to TST (green) of rate coefficients for hydride transfer catalyzed by Escherichia coli dihydrofolate reductase in 33% glycerol. The corresponding curve from equation (8), with a=96,νzTS/zR=1061.7s−1,ζ=106.7,a/ζ=0.90, is shown in blue and the experimental data in red.
(a)(b)


## 2. Results and Discussion

### 2.1. The Physical Limit to Enzyme Catalysis

We begin by clarifying the meaning of the frequency factorν by approximating the density of vibrational states of the catalytic protein with the Debye model that gives(11)ρωdω=9Nω2ωD3dω,with(12)ωD=6π2NV1/3vs=1.21×1013Hz,with vs being the speed of sound in a crystal of volume V with N atoms. For a protein with an average molecular weight of 125kDa [6], we obtain(13)∫0ω¯ρωdω=3Nω¯ωD3=8.17×10−21≪1,an indication that ν is not the frequency of any vibrational mode, but it is a reflection of the timescale ν−1 of the exchange of energy between active modes and the reaction coordinate. This can be exemplified by a simple system of two linear coupled oscillators with displacements x and y. An initial energy E0 is placed in the first oscillator, and the system evolves with the displacements shown in Figure 6, while plots of the potential energy of the two bonds are reported in Figure 7.Figure 6
Plots ofxt and yt versus time.Figure 7
Plots ofu1/E0 and u2/E0 versus time.We observe [8] that the energy transfer between oscillators takes place at ω2t≈25, clearly indicating the physical significance of the frequency ν as the inverse of the timescale for the energy transfer from the reaction medium (catalytic protein in the case of enzyme catalysis) to the substrate. Thus, the time dependence of the energy transfer among the oscillators is a complex function of the parameters of the system, i.e., the coupling constants, the masses, and the structure of the system with its geometrical constraints. All the physical parameters contribute to the frequency ν, a fundamental quantity for each enzyme that limits its catalytic efficiency.The derivatives ofQ with respect to a, ζ, and T are, respectively,(14)∂aQa,ζ=1Γa∫ζ∞xa−1lnxe−xdx−QΓa∫0∞xa−1lnxe−xdx,∂ζQa,ζ=−ζa−1Γae−ζ,∂TQa,ζ=−ζT∂ζQa,ζ.Plots of the above derivatives are shown in Figure8. The maximum of the temperature dependence of Q is obtained by setting the second derivative:(15)∂aT2Qa,ζ=∂TQa,ζlnζ−1Γa∫0∞xa−1lnxe−xdx=∂TQa,ζlnζ−dalnΓa,equal to zero. The digamma function dalnΓa can be approximated as(16)dalnΓa≃lna−12,and the solution of ∂aT2Qa,ζ=0 gives the optimal value of a at(17)a≃ζ+12,and the ratio a/ζ≃1+1/2ζ to be compared to the observed a/ζ≃0.9.Figure 8
Plots of∂aQa,ζ, ∂TQa,ζ, and ∂aT2Qa,ζ versus a−1/ζ.We present an interpretation of the experimental observationa/ζ≃0.9 in terms of the inability of a living system to maintain homeostasis for higher values of a, followed by a substantiation of equation (8) based on its consequence on the functional form of the chemical equilibrium constant.


### 2.2. Homeostasis

Enzymes operating in an autotrophic organism witha≫ζ would synthesize a sugar with concentration x with the rate Φsin2χ, where the parameter Φ expresses both the solar energy flux and the efficiency of the photosynthetic process. The average of the sun elevation sin2χ is time-dependent and averaged over one-day intervals. The same compound would be processed at a rate krx, resulting in the overall evolution of x given by the equation(18)x˙=Φsin2χ−krx.The sun elevation is given by(19)sinχ=cosHcosφcosδ+sinφsinδ,if(20)cosH+tanφtanδ>0,and sinχ=0 otherwise. The angles H and δ in equations (19) and (20) depend on time and represent the hour angle of the sun and its declination, respectively. The latitude φ is a constant. The solution of equation (18) is(21)x=x0e−krt+Φe−krt∫0tsin2χekrudu,and its plot is shown in Figure 9 (green curve). In the presence of homeostasis, the metabolism rate would also depend on the solar irradiation through temperature, and we would have(22)x˙h=Φsin2χ−krsin2χxh,which has the solution(23)xh=Φkr+x0−Φkrexp−kr∫0tsin2χdu,also shown in Figure 9 (blue curve).Figure 9
Plots ofxt (green) and xht (blue) versus time.Thus, a temperature-dependent rate coefficient for metabolism, that in turn is a consequence of the temperature dependence of the rate coefficient for enzyme catalysis, entails a more stable concentration of the nutrient that would not attain the low values endangering the survival of the organism. This temperature dependence, along with homeostasis, would be lost at high values (a≫ζ) of the number of coupled vibrational mode a.


### 2.3. The Equilibrium Constant

We can find further support to the validity of the model expressed by equation (8) by analyzing the discrepancy between the computed energy difference between reactants and products in the equilibrium(24)α−D−glucopyranose⇌k2k1β−D−glucopyranose,in water and the corresponding quantity obtained by regression to the thermodynamic expression(25)Keq=zPzRe−Δu/kT.Neglecting the weak temperature dependence of the ratio of the partition functions, the temperature dependence ofKeq would be(26)∂TlnKeq≃ΔukT2.Thus, the plot ofKeq versus temperature would exhibit a slope with the same sign as Δu. Hence, the negative slope of the plot in Figure 10 for the equilibrium constant of the reaction α−D−glucopyranose⇌β−D−glucopyranose in water implies the negative Δu=−0.79kcal mol−1 (green line zP/zR=0.40). On the other hand, the potential energies of all the conformations of β−D−glucopyranose calculated at the B3LYP/6-311 + + G∗∗ level of theory exhibit positive differences with respect to the ground state conformation α−gt−4C1 of α−D−glucopyranose and the β conformation with the lowest energy being the β−gg−4C1 with Δu=0.855kcal mol−1 [9]. There is thus a discrepancy between this computational result and the negative Δu obtained by the regression of equation (25) to experimental data.Figure 10
Plots ofKeq versus temperature according to equation (25) (green line) and (27) (blue curve).If we consider, instead, the expression for the equilibrium constant according to equation (8), i.e., the ratio of the rate coefficient for the forward reaction (with a1 active oscillators and reduced barrier ζ1) to the corresponding reverse reaction (with a2 active oscillators and reduced barrier ζ2), its form would be(27)Keq=ν1ν2zPzRQa1,ζ1Qa2,ζ2=ν^Q1Q2,with the temperature dependence(28)∂TlnKeq≃∂TlnQ1−∂TlnQ2,that can be negative also if Δu>0. In fact, the blue curve in Figure 10 is the best fit of equation (27) to the experimental data with a1=5 and a2=1, affording ν¯=0.153 and Δu=0.989kcal mol−1, in good agreement with the result given by theoretical calculations.


## 2.1. The Physical Limit to Enzyme Catalysis

We begin by clarifying the meaning of the frequency factorν by approximating the density of vibrational states of the catalytic protein with the Debye model that gives(11)ρωdω=9Nω2ωD3dω,with(12)ωD=6π2NV1/3vs=1.21×1013Hz,with vs being the speed of sound in a crystal of volume V with N atoms. For a protein with an average molecular weight of 125kDa [6], we obtain(13)∫0ω¯ρωdω=3Nω¯ωD3=8.17×10−21≪1,an indication that ν is not the frequency of any vibrational mode, but it is a reflection of the timescale ν−1 of the exchange of energy between active modes and the reaction coordinate. This can be exemplified by a simple system of two linear coupled oscillators with displacements x and y. An initial energy E0 is placed in the first oscillator, and the system evolves with the displacements shown in Figure 6, while plots of the potential energy of the two bonds are reported in Figure 7.Figure 6
Plots ofxt and yt versus time.Figure 7
Plots ofu1/E0 and u2/E0 versus time.We observe [8] that the energy transfer between oscillators takes place at ω2t≈25, clearly indicating the physical significance of the frequency ν as the inverse of the timescale for the energy transfer from the reaction medium (catalytic protein in the case of enzyme catalysis) to the substrate. Thus, the time dependence of the energy transfer among the oscillators is a complex function of the parameters of the system, i.e., the coupling constants, the masses, and the structure of the system with its geometrical constraints. All the physical parameters contribute to the frequency ν, a fundamental quantity for each enzyme that limits its catalytic efficiency.The derivatives ofQ with respect to a, ζ, and T are, respectively,(14)∂aQa,ζ=1Γa∫ζ∞xa−1lnxe−xdx−QΓa∫0∞xa−1lnxe−xdx,∂ζQa,ζ=−ζa−1Γae−ζ,∂TQa,ζ=−ζT∂ζQa,ζ.Plots of the above derivatives are shown in Figure8. The maximum of the temperature dependence of Q is obtained by setting the second derivative:(15)∂aT2Qa,ζ=∂TQa,ζlnζ−1Γa∫0∞xa−1lnxe−xdx=∂TQa,ζlnζ−dalnΓa,equal to zero. The digamma function dalnΓa can be approximated as(16)dalnΓa≃lna−12,and the solution of ∂aT2Qa,ζ=0 gives the optimal value of a at(17)a≃ζ+12,and the ratio a/ζ≃1+1/2ζ to be compared to the observed a/ζ≃0.9.Figure 8
Plots of∂aQa,ζ, ∂TQa,ζ, and ∂aT2Qa,ζ versus a−1/ζ.We present an interpretation of the experimental observationa/ζ≃0.9 in terms of the inability of a living system to maintain homeostasis for higher values of a, followed by a substantiation of equation (8) based on its consequence on the functional form of the chemical equilibrium constant.


## 2.2. Homeostasis

Enzymes operating in an autotrophic organism witha≫ζ would synthesize a sugar with concentration x with the rate Φsin2χ, where the parameter Φ expresses both the solar energy flux and the efficiency of the photosynthetic process. The average of the sun elevation sin2χ is time-dependent and averaged over one-day intervals. The same compound would be processed at a rate krx, resulting in the overall evolution of x given by the equation(18)x˙=Φsin2χ−krx.The sun elevation is given by(19)sinχ=cosHcosφcosδ+sinφsinδ,if(20)cosH+tanφtanδ>0,and sinχ=0 otherwise. The angles H and δ in equations (19) and (20) depend on time and represent the hour angle of the sun and its declination, respectively. The latitude φ is a constant. The solution of equation (18) is(21)x=x0e−krt+Φe−krt∫0tsin2χekrudu,and its plot is shown in Figure 9 (green curve). In the presence of homeostasis, the metabolism rate would also depend on the solar irradiation through temperature, and we would have(22)x˙h=Φsin2χ−krsin2χxh,which has the solution(23)xh=Φkr+x0−Φkrexp−kr∫0tsin2χdu,also shown in Figure 9 (blue curve).Figure 9
Plots ofxt (green) and xht (blue) versus time.Thus, a temperature-dependent rate coefficient for metabolism, that in turn is a consequence of the temperature dependence of the rate coefficient for enzyme catalysis, entails a more stable concentration of the nutrient that would not attain the low values endangering the survival of the organism. This temperature dependence, along with homeostasis, would be lost at high values (a≫ζ) of the number of coupled vibrational mode a.


## 2.3. The Equilibrium Constant

We can find further support to the validity of the model expressed by equation (8) by analyzing the discrepancy between the computed energy difference between reactants and products in the equilibrium(24)α−D−glucopyranose⇌k2k1β−D−glucopyranose,in water and the corresponding quantity obtained by regression to the thermodynamic expression(25)Keq=zPzRe−Δu/kT.Neglecting the weak temperature dependence of the ratio of the partition functions, the temperature dependence ofKeq would be(26)∂TlnKeq≃ΔukT2.Thus, the plot ofKeq versus temperature would exhibit a slope with the same sign as Δu. Hence, the negative slope of the plot in Figure 10 for the equilibrium constant of the reaction α−D−glucopyranose⇌β−D−glucopyranose in water implies the negative Δu=−0.79kcal mol−1 (green line zP/zR=0.40). On the other hand, the potential energies of all the conformations of β−D−glucopyranose calculated at the B3LYP/6-311 + + G∗∗ level of theory exhibit positive differences with respect to the ground state conformation α−gt−4C1 of α−D−glucopyranose and the β conformation with the lowest energy being the β−gg−4C1 with Δu=0.855kcal mol−1 [9]. There is thus a discrepancy between this computational result and the negative Δu obtained by the regression of equation (25) to experimental data.Figure 10
Plots ofKeq versus temperature according to equation (25) (green line) and (27) (blue curve).If we consider, instead, the expression for the equilibrium constant according to equation (8), i.e., the ratio of the rate coefficient for the forward reaction (with a1 active oscillators and reduced barrier ζ1) to the corresponding reverse reaction (with a2 active oscillators and reduced barrier ζ2), its form would be(27)Keq=ν1ν2zPzRQa1,ζ1Qa2,ζ2=ν^Q1Q2,with the temperature dependence(28)∂TlnKeq≃∂TlnQ1−∂TlnQ2,that can be negative also if Δu>0. In fact, the blue curve in Figure 10 is the best fit of equation (27) to the experimental data with a1=5 and a2=1, affording ν¯=0.153 and Δu=0.989kcal mol−1, in good agreement with the result given by theoretical calculations.


## 3. Conclusions

(1)
A proposed model for the expression of the rate coefficient of enzyme-catalyzed reactions is in agreement with the observed negative curvature of the Eyring plots for these reactions.(2)
An explanation for the observed relation between the number of active modesa and the reduced barrier ζ (a/ζ≃0.9) is given in terms of the necessity of maintaining homeostasis in the early autotrophic organisms.(3)
An analysis of the experimental measurements of the equilibrium constant at various temperatures for the reversible reactionα−D−glucopyranose⇌β−D−glucopyranose in water affords a negative slope. The interpretation of this negative slope with the conventional expression for the equilibrium constant affords a negative potential energy difference between reactants and products, at odds with the result of theoretical calculations. On the other hand, consistently with the prediction of theoretical calculations, the expression for the equilibrium constant derived from the proposed model gives a positive potential energy difference between reactants and products even in the presence of a negative slope of the equilibrium constant at various temperatures.



---

*Source: 2524862-2020-06-27.xml*