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

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. $e 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. $is 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. $e discrepancy is resolved by making use of the proposed model.


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]. e model was tested on a series of both catalyzed and uncatalyzed reactions, for which the temperature-dependent rate coefficients have been measured [3]. e 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. e thermally averaged rate coefficient of a chemical reaction is expressed as the integral of the product of an energy-dependent rate coefficient k r (ε) and the classical Boltzmann distribution of low-frequency vibrational modes: e lower limit of the integral is the threshold value u TS of the potential energy that the substrate must exceed to be able to cross the reaction barrier. e functional form of the energy-dependent k r (ε) 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 (]): e total number of modes s includes the subset of all the excited modes that can transfer energy to the reaction coordinate (indicated by a), implying (Zω i /kT) ≪ 1, and the classical expression for the vibrational density of states leading to the partition function e 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: in the variables u � (ε/ε) and x � (ε/kT), where the last equality defines the regularized beta function (B(a, s − a, (ζ/x))/B(a, s − a)). Taking the anharmonicity of the lowfrequency 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. e thermally averaged rate coefficient in equation (1) thus becomes with ζ � (u TS /kT). Making use of the theorem [5], we obtain the expression, independent of s, for the rate coefficient in terms of the regularized gamma function Q: is expression is plotted in Figure 1, and it exhibits an inflection point at ζ � a − 1. e above expression for k r 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 75 kDa. 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.
For integer values of a, we can integrate Q analytically, obtaining At the inflection point, the value of the function Q is given by and its values are plotted in Figure 2 as a function of a. Figures 3-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],  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]. e analysis of experimental data for these and other cases of enzyme-catalyzed reactions indicates that (a/ζ) ≃ 0.9, and the most common values of Q(a, ζ) are slightly less than 1/2, since a typical value of the reduced barrier ζ usually exceeds 40. is prompts the question of why the evolutionary pressure did not push the structure of enzymes to reach the maximum efficiency at Q(a, ζ) ≃ 1.

e 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 with with v s being the speed of sound in a crystal of volume V with N atoms. For a protein with an average molecular weight of 125 kDa [6], we obtain 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. is can be exemplified by a simple system of two linear coupled oscillators with displacements x and y. An initial energy E 0 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.
We observe [8] that the energy transfer between oscillators takes place at ω 2 t ≈ 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. us, 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. e derivatives of Q with respect to a, ζ, and T are, respectively,   Journal of Chemistry Q(a, ζ).

(14)
Plots of the above derivatives are shown in Figure 8. e maximum of the temperature dependence of Q is obtained by setting the second derivative: equal to zero. e digamma function d a ln Γ(a) can be approximated as and the solution of z 2 aT Q(a, ζ) � 0 gives the optimal value of a at and the ratio (a/ζ) ≃ 1 + (1/2ζ) to be compared to the observed (a/ζ) ≃ 0.9. We present an interpretation of the experimental observation (a/ζ) ≃ 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.

Homeostasis.
Enzymes operating in an autotrophic organism with a ≫ ζ would synthesize a sugar with concentration x with the rate Φ〈sin 2 χ〉, where the parameter Φ expresses both the solar energy flux and the efficiency of the photosynthetic process. e average of the sun elevation 〈sin 2 χ〉 is time-dependent and averaged over one-day intervals. e same compound would be processed at a rate k r x, resulting in the overall evolution of x given by the equation e sun elevation is given by and sin χ � 0 otherwise. e angles H and δ in equations (19) and (20) depend on time and represent the hour angle of the sun and its declination, respectively. e latitude φ is a constant. e solution of equation (18) is 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 which has the solution also shown in Figure 9 (blue curve). us, 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.
is temperature dependence, along with homeostasis, would be lost at high values (a ≫ ζ) of the number of coupled vibrational mode a.

e Equilibrium Constant.
We can find further support to the validity of the model expressed by equation (8)  analyzing the discrepancy between the computed energy difference between reactants and products in the equilibrium in water and the corresponding quantity obtained by regression to the thermodynamic expression Neglecting the weak temperature dependence of the ratio of the partition functions, the temperature dependence of K eq would be us, the plot of K eq 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.79 kcal mol − 1 (green line (z P /z R ) � 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− 4 C 1 of α − D − glucopyranose and the β conformation with the lowest energy being the β − gg− 4 C 1 with Δu � 0.855 kcal mol − 1 [9]. ere is thus a discrepancy between this computational result and the negative Δu obtained by the regression of equation (25) to experimental data.
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 a 1 active oscillators and reduced barrier ζ 1 ) to the corresponding reverse reaction (with a 2 active oscillators and reduced barrier ζ 2 ), its form would be with the temperature dependence 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 a 1 � 5 and a 2 � 1, affording ] � 0.153 and Δu � 0.989 kcal mol − 1 , in good agreement with the result given by theoretical calculations.

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. 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.

Data Availability
e data used to support this study are included within the article.