Triosephosphate Isomerase (TPI): The Favorite Enzyme for Computational Optimization of Michaelis-Menten Type Kinetics
Triosephosphate isomerase (TPI, EC 5.3.1.1) is an essential enzyme in glycolysis [
51,
52]. Its central housekeeping role is very fast catalytic interconversion of dihydroxyacetone phosphate (DHAP) and glyceraldehyde-3-phosphate (GAP). There would be no net yield of ATP from anaerobic glucose metabolism without the TPI forward activity (DHAP → GAP). Of all enzyme-catalyzed reactions, the free-energy profile was first determined for TPI [
53]. The seminal works of Jeremy Knowles [
32], John Albery [
54], and other authors described the TPI as a perfect enzyme in the sense that it is the perfectly evolved enzyme with catalytic efficiency close to the diffusion limit. In 1984, John Richard [
55] estimated that k
cat/K
M for TPI increased 3∙10
10 times compared to the inorganic DHAP to GAP conversion. Enzyme efficiency inside the diffusion limit was confirmed for the wild-type TPI enzymes isolated from many species [
56].
As a reversible enzyme working close to thermodynamic equilibrium, the TPI can be easily induced to work in the backward direction (GAP → DHAP). Its central physiological role is maintaining the delicate balance between glycolysis and gluconeogenesis. However, since TPI belongs to the most ancient enzymes [
57], the biological evolution involved it in the pentose phosphate pathway, triacylglyceride accumulation, and many other moonlighting functions [
58,
59]. With such a broad spectrum of activities and functions, it is not surprising that the TPI enzyme has attracted the medical community's interest. TPI inhibitors are promising as antiprotozoal drugs for the treatment of diseases caused by
Trypanosoma cruzi, Trypanosoma brucei, Plasmodium falciparum, Giardia lamblia, Leishmania mexicana, Trichomonas vaginalis, and
Entamoeba histolytica [
60]. The upregulation of the TPI gene is common in many cancers [
61]. At the same time, TPI deficiency or reduced activity causes the accumulation of DHAP connected to severe diseases, such as hemolytic anemia, recurrent infections, cardiomyopathy, and fatal neuromuscular dysfunction [
62].
Stephen Blacklow asserted [
63] that TPI enzyme "can improve no further as a catalyst, " assuming constraints of free diffusion and in vivo levels of its substrates. In the meantime, researchers proposed electrostatic screening [
64,
65], TPI oligomerization [
66], elevated temperature for the TPI from thermophilic cells [
67], and other mechanisms [
68] how TPI catalytic efficiency can be increased above observed values. Ideally, the mutations or modifications making TPI more resistant to oxidative damage and a more efficient catalyst can help prevent and treat Alzheimer's disease [
68,
69].
We stressed in our previous contributions [
70,
71] that increasing the TPI catalytic turnover and efficiency above observed "perfect "values is theoretically possible when enzyme kinetics is connected to the maximal partial entropy production principle from irreversible thermodynamics [
2]. In this chapter, we shall attempt to answer the following questions: a) Does TPI performance change after noise is taken into account? b) If it does change, is it possible to find the combination of microscopic rate constants resulting in at least ten-fold increased performance regarding the k
cat/K
M value calculated from the experimental data? c) How is the entropy production by TPI related to corresponding enzyme efficiency values? d) Are any published optimization methods better at finding high forward k
cat/K
M values than different means of noise introduction?
Let us first present observed values for TPI kinetic parameters [
32,
70] to easily compare all our simulations with the experimental values (
Table 1). Triosephosphate isomerase can be found in four functional states [
54]. According to
Figure 1c) 1 is the free enzyme (E), 2 is the enzyme-substrate bound complex (ES), 3 is a transition state intermediate (EZ), and 4 is the enzyme-product bound complex (EP). The reference steady state [
54] is such that the concentration of substrate is [S] = 40 μM and the concentration of product is [P] = 0.064 μM. The values of the kinetic constants k
1 and k
8 in
Table 1 are obtained respectively from expressions k
1= k
1*∙[S] and k
8 = k
8*∙[P], where second-order rate constants k
1* and k
8* are measured in (Ms)
-1.
The initial TPI concentration in our simulations ranged from 10 to 50 nM. Mass conservation for all enzyme conformations is always taken into account in all simulations. All NetLogo programs also required the mass conservation of ligands (substrates, products, and their intermediate TPI-bound forms). That requirement entered the FORTRAN programs as the [S]+[P] = constant condition when we allowed for changes in the concentrations of ligands. The concentration of bound ligands [ES]+[EZ]+[EP] is always much smaller than [S]initial + [P]initial concentration because bound ligands concentration cannot exceed the initial low concentration of free TPI enzymes. Thus, the mass conservation of ligands is considered a good approximation in those FORTRAN programs that examined how different parameters change after changes in the substrate and product concentrations.
Stepwise increases of rate constants from the product-release transition
Let us first consider how catalytic efficiency depends on overall entropy production in a deterministic manner when the implicit assumption is that noise does not exist. FORTRAN program is convenient to use for such a study. For constant temperature, the dissipation function φ and total entropy production
P have the absolute temperature T as the proportionality factor: φ = T∙
P. Assuming that
P is not the conseqeunce but the cause for the catalytic efficiency (see Introduction), we use either entropy production or dissipation term to plot the functional relationship between k
cat/K
M values at the y-axis and the dissipation/RT values at the x-axis. The first such plot (
Figure 2) illustrates how TPI efficiency changes after the stepwise increase in the microscopic rate constant k
7. All other rate constants and equilibrium constants K
1, K
2, and K
3 are kept at their observed values (see the calculated values of rate constants and the values for the initial concentrations of substrates and products from
Table 1. Since the equilibrium constant K
4 = k
7/k
8 also goes through the stepwise increase, the expected outcome of the first simulation scenario is a regular increase in the chemical affinity or force (expressed as X
tot/RT values) from negative to positive values.
Negative force values correspond to negative backward flux (GAP → DHAP) and positive dissipation, while positive force values correspond to positive forward flux (DHAP → GAP) and positive dissipation. Both limits in the force range, negative and positive, are associated with the high dissipation. Still, only the positive limit corresponds to the maximal enzyme efficiency value of 1.25∙10
6 M
-1s
-1 (
Figure 2). That result is an encouraging 1.59-fold increase over the observed value of 7.9∙10
5 M
-1s
-1 (corresponding to the X
tot/RT = 0.685), but not the significant improvement over the 1.13∙10
6 M
-1s
-1 value we obtained in the 2017 [
70].
From the output of the same FORTRAN program, it is easy to select only the positive X
tot/RT values. The resulting efficiency dependence on dissipation is then well correlated (R
2 = 0.944) with the straight-line proportionality (
Figure 3) Thus, from zero forward catalytic efficiency and vanishing entropy production in the thermodynamic equilibrium, there must be an obligatory increase in the dissipation, which is tightly coupled to the catalytic efficiency increase. We did not ask how one can achieve the increase in only the chosen kinetic constant k
7 in practice without any other change. It is unlikely that random or intentional mutations can ever do it. However, fine-tuning microwave irradiation may potentially produce the non-thermal effect of significantly accelerating the product release catalytic step. It is easier to answer why simulations presented in
Figure 2 and
Figure 3 dealt with the k
7 stepwise increase. We assumed, as in [
70], that the product release rate limits the TPI catalytic power. In our notation for rate constants (see
Figure 1c), the k
7 is the first-order rate constant, determining the product release rate.
We can also explore the stepwise increase of k
7 and k
8 when all equilibrium constants and all other rate constants maintain their observed values (see
Table 1). Almost perfect proportionality is obtained between k
cat/K
M and corresponding entropy production values (
Figure 4).
The next task is to answer how thermodynamic and kinetic parameters change for the TPI catalytic cycle when we limit deterministic changes to decreasing substrate and increasing product concentrations. The answer is provided in
Figure 5, which illustrates how net flux and overall dissipation vary with the force changes.
The catalytic activity optimizations in the forward direction, when the substrate is converted into the product, are better connected with the physiological role of the TPI enzyme in glycolysis. We published one example of such optimization in 2017 [
70]. It was for the fixed positive force (chemical affinity) corresponding to X
tot/RT = 0.685 (the vertical line in the insert of
Figure 5). The optimization example for the reverse process (product-to-substrate conversion) leads to decreased catalytic efficiency for the forward process. The dissipation and net flux for the reverse process increase by several orders of magnitude when the applied force has a high negative value (the vertical line in the main figure). It is a pathological situation with no connection to TPI's role in cellular metabolism.
Šterk et al. result [
72], after maximizing the total entropy production density, was J = - 1272 s
-1 (blue point at the vertical line in the main figure). It is about 100 times higher net reaction flow in the reverse GAP→DHAP direction when compared to experimentally observed reaction rate J = 14 s
-1 facilitating glycolysis [
32,
70]. The maximal entropy production of Šterk et al. [
72] is σ/nR = 5685 s
-1 (the red point at the vertical line in the main figure). Klipp and Heinrich [
73] obtained even higher net reaction flow in the backward GAP→DHAP direction ranging from J = -1620 (for experimental rate constants values when [DHAP] = [GAP] = 40 μM) to J = - 4010 s
-1 (for the separate limit optimization model), the result that was verified and commented by Bish and Mavrovouniotis [
74]. These optimizations for highly negative flux and negative force can only ensure the non-physiological operation of the TPI enzyme and the loss of its primary function of balancing glycolysis and gluconeogenesis. For instance, optimized k
cat in the forward direction of Šterk et al. [
72] is k
cat = 222 s
-1, which is worse then k
cat = 432 s
-1 (experimental data [
32]), while our optimized k
cat = 686 s
-1 [
70] is an improvement over the k
cat value calculated from experimental data.
It all depends on the choice of the optimization procedure. We chose to maximize the partial entropy production in the rate-limiting product-release step (the 4th catalytic step in the forward direction) [
2]. We noticed in 2017 [
70] how that choice led to the concomitant increase in the optimal net flux (from 14.4 to 20.77 s
-1), optimal catalytic constant (from 432 to 686 s
-1), optimal catalytic efficiency (from 7.86∙10
5 to 1.13∙10
6 M
-1s
-1), and optimal overall entropy production (from 9.9 to 14.2 s
-1). Within the restriction we used (fixed equilibrium constants for each catalytic step at their values calculated from the experimental data), there was a common 30% increase in the flux, efficiency, and dissipation.
Figure 4 above illustrates how a regular 30% increase between two points follows from the constant slope and perfect proportionality between enzyme kinetic parameters and its overall entropy production. We did not show the J dependence on dissipation because the dissipation function φ = J∙X, and force X was constant in the calculations we performed to construct that figure.
In the following subsection, we study the same efficiency-dissipation relationships when noise is introduced in the last catalytic step so that FORTRAN programs for calculating kinetic and thermodynamic parameters do not contain regular stepwise increases in catalytic constants k7 and k8.
Noise introduction in kinetic constants with selected restrictions
We asked what would happen after introducing random normal noise in forward (k
7) and backward (k
8) rate constants for the rate-limiting product-releasing step.
Figure 6. illustrates the advantage of using noise when looking for the combination of rate constants corresponding to higher enzyme efficiency. The highest efficiency of 1.6∙10
6 M
-1s
-1 is associated with the highest total dissipation in the RT units (20.3 s
-1) due to the perfect proportionality between the main enzyme performance parameter and the main physical parameter in irreversible thermodynamics. We used the same restrictions of constant overall force and constant equilibrium constants K
1 to K
4, which were required in deriving the partial entropy production theorem [
2,
30,
70]. Random normal noise was called once in the FORTRAN program as the Box-Muller transform (see Methods, eq (25)) with the shift +2 to ensure that only positive rate constants k
7 are the output. There was no need to call that function again for the multiplication with the observed k
8 value because we kept the no-change requirement for all equilibrium constants K
i (i = 1, 2, 3, 4) from our 2017 paper [
70].
In the next computational experiment, we introduce normal noise in the kinetic constant k
7 = 4000 s
-1 while taking representative initial values of rate constant k
8 = 25 s
−1, 32 s
−1, 40 s
−1, 100 s
−1, and 160 s
−1. Fixed equilibrium constants for each k
7-k
8 pair are then K
4 = 160, 125, 100, 40, and 25, respectively, each calculated using experimental value k
7 = 4000 s
−1. Experimental data for kinetic constants of the enzyme triosephosphate isomerase shown in
Table 1 are used for transitions other than transition 4. Enzyme efficiency k
cat/K
M as a function of dissipation/RT is shown in
Figure 7 for the forces X/RT equal to 0.2389, 0.4620, 0.6852, -0.6774, and -1.1474 corresponding to the equilibrium constants K
4=100, 125, 160, 40 and 25, respectively.
There is one additional condition besides those we used to create
Figure 6. We assumed a constant sum of free substrates and free products. It is a good approximation for the mass conservation of ligands only if the initial free enzyme concentration (50 nM) is much smaller than the concentrations of [S]+[P] for all points and all forces.
Figure 7 joins the results of five FORTRAN programs that include noise in the last forward catalytic step.
Careful examination of the case X
tot/RT = -1.1474 reveals slight curvature in the efficiency as a function of dissipation (magenta symbols). The ratios of rate constants k
3/k
7 and k
6/k
7 (from the K
M expression [
38]) are not constants in the output file because k
3 and k
6 are constants but k
7 changes in a random manner. The best linear fit slope increases when negative force values approach the thermodynamic equilibrium and then decreases when positive force is increased. Thus, we examined slope changes and goodness of linear fit changes for a wider span of force values ranging from -3 to +4 (
Figure 8). For that task, we constructed ten FORTRAN programs. The ratios of rate constants k
3/k
7 and k
6/k
7 (from the eq. (23) K
M expression) are not constants in the output files because k
3 and k
6 are constants, but k
7 changes randomly. The K
M exhibits small changes.
Figure 8 shows the output of these programs. It illustrates how the slope and the perfection of the seemingly straight-line proportionality increase with the approach to the thermodynamic equilibrium when net force and entropy production vanish.
Computational Optimizations of the TPI Catalytic Activity when noise is included
The question we studied in this subsection is how noise introduction affects various computational optimizations for TPI catalysis. In
Figure 9, variations of K
1 and K
4 were introduced by multiplication of K
4 = k
7/k
8 (see
Table 1) with the normal noise. The fixed force restriction X = X
tot/RT = 0.684 [
70] ensured concomitant variations in K
4 and K
1. There was no explicit requirement for the maximal entropy production. Still, after going randomly through the 1000 quasi-steady states, our FORTRAN program finds that the maximal overall dissipation corresponds to optimal enzyme efficiency (
Figure 9). The 7-fold efficiency improvement from 7.86∙10
5 to 5.585∙10
6 M
-1s
-1 follows after a 4-fold dissipation increase (see
Table 1).
The best combination of the backward rate constants k
2 and k
8, which resulted in even higher k
cat/K
M of 8.903∙10
6 M
-1s
-1, is k
2 = 74 s
-1 and k
8 = 2438 s
-1. The enzyme working in that state has 11 times higher catalytic activity (the highest point in
Figure 9) than the value of 7.86∙10
5 M
-1s
-1 calculated from the experimental data (
Table 1). Required changes in rate constants are two orders of magnitude changes in k
2 (decrease) and k
8 (increase). These rate changes describe the inhibition of substrate release from the ES complex and stimulation of product association with the free enzyme. The corresponding overall dissipation per RT of 21.3 s
-1 is approximately double the value calculated from the experimental data. Still, the dissipation needed to reach the maximal efficiency state is halved compared to maximal dissipation (
Figure 9).
Interestingly, the same dissipation value of 21 to 22 s-1 is connected with two very different catalytic efficiency values of 8.9∙106 and 1.94∙106 s-1, respectively. Thus, when specific restrictions are imposed, the nonlinear system may be able to jump between two quasi-steady states characterized by high and low efficiency and a minor change in dissipation. How to force the system to live in about a 10-fold higher efficiency state with only a two-fold higher price in terms of overall dissipation is outside the scope of this paper.
Optimal efficiency can be obtained for fixed force when other pairs of equilibrium constants are varied by introducing noise. We did not show corresponding efficiency (dissipation) dependence because optimal k
cat/K
M values for the dissipation maximum were considerably lower from the 8.9∙10
6 M
-1s
-1 value obtained after K
1-K
4 variations. A reader can verify that conclusion from the
Figure 10 coordinates (20.6, 1.95∙10
6) and (21.7, 1.9∙10
6) obtained after K
2-K
4 and K
3-K
4 variations. Still, the requirement that total entropy production is maximal and corresponding restrictions on equilibrium constants for chosen catalytic steps can produce higher catalytic efficiencies for fixed force than the maximal requirement for selected partial entropy productions (
Figure 10).
The primary purpose of
Figure 10 is to illustrate the relationships among different methods for obtaining higher than referent values for catalytic efficiency (
Table 1). That task led to the map of dissipation-efficiency points when the x-axis is for dissipation and the y-axis is for enzyme efficiency. The perfect efficiency-dissipation proportionality is the straight line fit to 1000 points after each kinetic constant is multiplied with the normal noise invoked only once in the corresponding FORTRAN program. It is the consequence of assuming fixed values for all equilibrium constants K
i, meaning that the overall force is also identical for all data points (their referent values can be found in
Table 1). We also used those assumptions in our previous publications [
2,
30,
70]. We did not consider noise in these publications. The first two highlighted points (9.9, 0.79∙10
6) and (14.2, 1.13∙10
6) are centered at the linear fit. They are the dissipation and efficiency values calculated from the experimental data and the modest improvement achieved after the requirement that partial entropy production P
4 in the rate-limiting product release step is maximal [
2,
70,
71].
We discussed above the results after introducing noise in the pairs of equilibrium constants K
1-K
4, K
2-K
4, and K
3-K
4. These are off-line points in
Figure 10, respectively: (39.4, 5.585∙10
6), (20.6, 1.95∙10
6), and (21.7, 1.9∙10
6). Thus,
Figure 10 clearly shows the advantage of noisy substrate and product association with free enzyme (the highest point). After considering many different optimization methods for entropy production (either ours or by other authors), the K
1-K
4 variations with constant force restriction resulted in best theoretical increase of the TPI catalytic efficiency above its observed value [
32].
Computational optimizations of TPI kinetics by some other authors [
72,
73] used the substrate and product concentrations similar to each other. It reversed the net flow in the direction of product → substrate due to negative flux and force and resulted in higher total entropy production values of several orders of magnitude. For instance, Šterk et al. [
72] used the constraint k
1*∙k
3∙k
5∙k
7 = K
+ = constant equal to the observed value. Such optimization required that the product of all kinetic constants in the forward direction and all kinetic constants in the backward direction k
2∙k
4∙k
6∙k
8* remain fixed when other parameters change. However, Šterk et al.[
72] used the steady state concentrations [S] = 31.45 and [P] = 8.55 μM. The corresponding force is then highly negative X
tot/RT = - 4.47. When multiplied with the high negative flux, it produced such a high dissipation that the optimal values (5685, 1.8∙10
6) could not be illustrated as the (x, y) point within the confines of
Figure 10. The authors found the maximum in overall entropy production, but it was about 570-fold higher than the calculated value from the experimental data. As expected, for the backward-directed enzyme turnover, the corresponding optimal efficiency for forward catalysis of 1.8∙10
6 M
-1s
-1 (the right-hand arrow pointing outside
Figure 10) is substantially smaller than our best results.
Noise introduction without restrictions other than all ki>0
When normal noise without shift is introduced in all rate constants k
i, some can vanish or become negative. To avoid such cases, we replaced negative with observed k
i values (see
Table 1).
Figure 11 illustrates that reasonable proportionality exists between efficiency and entropy production when there are no other restrictions on kinetic constants and on equilibrium constants for the TPI enzyme. The advantage of calling random numbers eight times (once for each of eight kinetic constants) is an extended range of possible steady states and forces. The highest efficiency state has 30-fold better efficiency and 160-fold higher dissipation compared to values calculated from experiments. The corresponding force for that state is X
tot/RT = 6.335.
The basic assumption we used in calculating entropy production values is that each of the 10000 computational steps probes a new quasi-steady state in which all parameters of interest can be calculated by using the Terrel Hill method [
34,
35]. We found the maximal efficiency value in the 1078th step. It corresponds to an unusually high information entropy of 1.181 and a low Michelis-Menten constant K
M = 0.000015. Interestingly, only the kinetic constants k
2, k
6 and k
7 significantly differed from their experimental values, all being much smaller, 56, 15, and 4 times, respectively. An increase in the k
1 value (from 400 to 1144 s
-1) may have resulted from increased substrate concentration or increased second-order rate constant for the association between the substrate and enzyme to form the ES complex. There was no change from experimental values for the kinetic constants. k
4, k
5, and k
8.
Enzyme turnover became slightly slower (kcat decreased from 432 to 348 s-1), but the division with considerably smaller KM (from 5.5∙10-4 to 1.474∙10-5) ensured surprisingly high efficiency. As is usually the case, the most illustrative representation is the profile of changes in the equilibrium constants or free energy changes. The equilibrium constant K1 increased about 160 times (from 0.057 to 9.07), and the K3 constant increased nearly 15 times (from 0.667 to 9.75). It led to a significant increase in the total equilibrium constant (from 1.98 to 564) despite a decrease in the K2 (from 0.333 to 0.174) and K4 (from 156.25 to 36.59).
There were 3580 points corresponding to the force X = X
tot/RT ≤ 0. Thus, for 35.8% of sets with random values for kinetic constants, the enzyme can still work in the reverse direction, converting products into substrates. Most k
i octuplets simulated the major physiological role of the TPI enzyme in converting DHAP to GAP. The best case of k
cat/K
M = 2.36∙10
7 M
-1s
-1 is also for the forward-directed net flux. However, we used the same forward catalytic efficiency definition for X > 0 and X ≤ 0. All the experimental data in the literature were extracted for the force X > 0 and flux J > 0 under the conditions when substrate concentration greatly exceeded the product concentration. The initial concentrations were [S] = 40 μM and [P] = 0.064 μM. Variations in k
1 and k
8 allowed the changes either in the second-order rate constants or in concentrations. The extreme case when X = -10.19 was obtained with k
1 = 3.4 s
-1 and k
8 = 27.0 s
-1. If the change in k
1 occurred only due to the change in [S], the substrate concentration would decrease almost 120 times. Therefore, although we included the points with negative force and flux in this figure and other simulations from the literature considered such cases [
72,
73], there is no experimental or physiological justification for retaining them.
Simulating dynamics using an agent-based modeling approach
Dynamics can be simulated using an agent-based modeling approach without solving differential equations. Agent-based model (ABM) is stochastic by nature. For instance, the stochastic noise inherent to the NetLogo computer language can be used to construct models for the stochastic interaction of an enzyme with its substrates, products, and inhibitors [
75,
76]. In the following example, noise is introduced in the TPI kinetics through random-float values (uniform noise) added to selected rate constants, not by Gaussian random number values (
Figure 12). Additional noise in rate constants is due to random encounters among ligands and [E]
free and among enzyme conformations [ES]↔[EZ] and [EZ]↔[EP], which is also specified with several different random-float values. Random changes occur in all computational steps (named „ticks“). Ticks can be in chosen time units. Agent-based programming requires dimensionless numbers as the input. However, when these numbers are specified as 40000 for substrates, 64 for products, and 50 for enzymes (for the TPI kinetics), they correspond to [S]
initial = 40 μM, [P]
initial =0.064 μM, and [E]
initial = 0.1 μM. The mass conservation of all ligand forms ([S], [P], [ES], [EZ], and [EP]) and all enzyme forms ([E]
free, [ES], [EZ], and [EP]) is an explicit requirement for each tick in all our NetLogo programs. Thus, [S]
initial + [P]
initial = [S] + [P] + [ES] + [EZ] + [EP] and [E]
initial = [E]
total = [E]
free + [ES] + [EZ] + [EP] because we left the system to itself and never added ligands or enzymes.
The intermediate equilibrium constants K
2 and K
3 never changed from their experimental values. Due to slight changes in the substrate and product concentrations the equilibrium constants K
1 and K
4 also underwent small changes. However, we kept the assumption that
where
is the second-order rate constant for the enzyme and substrate association, while const1 is the ratio determined from the observed
value, which does not change during the simulation. Similarly,
where
is the second-order rate constant for the enzyme and product association, while const2 is the ratio determined from experiments, which does not change during the simulation.
Since initial product concentration is small (64 nM), each stepwise increase in the product concentration is seen as a jump from one straight line fit to another in four steps „a“ to „d. “ It increased the product concentration to 67 nM. Thus, the proportionality between enzyme efficiency and entropy production (dissipation) remained almost perfect. Maximal efficiency values close to the 3∙106 M-1s-1 are about 4-fold higher compared to the value calculated from experiments. Similar 4-fold increase is for the corresponding dissipation.
When simulation time was extended to 2137 ticks, product concentration increased from 64 to 74 nM, while the driving force decreased from X
tot/RT = 0.68 to X
tot/RT = 0.54 with the same stepwise slope increase for the efficiency-dissipation dependence (
Figure 13). The best efficiency value of 3.3∙10
6 M
-1s
-1 corresponded to the dissipation/RT = 27.2 s
-1 . Free enzyme concentration dropped from 100 to 12 nM for this case of more extended simulation.
The insight from these figures would be that maximal catalytic efficiency remains approximately the same during the system relaxation toward the thermodynamic equilibrium (when chemical affinity, net flux, and total entropy production all reach their zero values). The slope of the efficiency-dissipation line keeps increasing toward an infinitely high value at the thermodynamic equilibrium when dissipation vanishes. Also, the perfection of the straight line approximation for the fit connecting all (x, y) values keeps increasing in discrete jumps (for each unit change in the product concentration) while chemical affinity decreases. The same time-development rule holds when equilibrium is spontaneously approached from the high positive or negative initial forces (see
Figure 8). Better efficiency to dissipation proportionality for positive forces stems from the
definition of catalytic efficiency, where both the catalytic constant and the Michaelis-Menten constant are defined for the forward direction
. It is easy to see that nonlinear is considerably better than linear fit for efficiency-dissipation proportionality in the case of higher negative X (see
Figure 8).
For longer simulations, the concentrations of enzyme conformations ES, EZ, and EP after each step (tick) go through the typical Michaelis-Menten kinetics: slow initial increase, faster, nearly constant increase, a broad maximum with minor changes, and prolonged decrease. That pattern repeats itself with the ES complex, after some delay with the EZ complex, and finally with the EP complex.
We next examined if a broader scope search for better enzyme performance is possible when Gaussian noise g
i (see Methods) is multiplied with each microscopic rate constant k
i (
Figure 14). The best catalytic efficiency of k
cat/K
M = 2.22∙10
7 M
-1s
-1 is indeed better than previous NetLogo simulations and similar to the best result we obtained after FORTRAN simulation for the TPI kinetics (
Figure 11).
Ketosteroid Isomerase (KSI) Case: What is Different when the Operating Range is Farther from Equilibrium?
Paul Talalay discovered in 1951 [
77] the
Pseudomonas testosteroni bacterium (presently named
Commamonas testosteroni [
78]) from soil beneath a rosebush on the Berkeley campus. The bacterium could grow in a medium containing testosterone as its only carbon and energy source. That was a clever and brave approach because, at that time, many steroid metabolites were known, but enzymic transformations of steroid hormones and metabolites were yet undiscovered. Paul Talalay and his collaborators purified highly active small bacterial enzyme ketosteroid isomerase from that bacterium and reported their findings from 1955 onward [
77]. The alternative name for the KSI enzyme is 3-oxo-Δ
5-steroid isomerase (EC:5.3.3.1).
Anna Radzicka and Richard Wolfenden reported typical high values for the
catalytic constant, catalytic efficiency, and catalytic proficiency of KSI as respectively 6.6∙10
4 s
-1, 3.0∙10
8 M
-1s
-1, and 1.8∙10
15 M
-1 [
12]. The catalytic proficiency is the catalytic efficiency rate enhancement (k
cat/K
M)/k
uncat when a nonenzymatic reaction rate constant k
uncat can be found for a corresponding spontaneous chemical reaction without the enzyme (1.7∙10
-7 s
-1 in our case). Thus, KSI is one of the fastest enzymes with extraordinary catalytic power. The formation of essential steroid hormones would take months to millions of years without enzymes such as KSI [
79]. The equilibrium constant K
eq = 2400 [
80] corresponds to the far from equilibrium conditions, high positive force, and the preference for the forward isomerization rate of 5-androstene-3,17-dione (a substrate for KSI) to its conjugate isomer 4-androstene-3,17-dione. Elucidating how the KSI reaction mechanism is connected to structure, kinetics, electrostatics, and thermodynamics was a highly challenging but worthy task through the last 50 years [
81,
82,
83,
84]. Hopefully, the rational design of KSI enzymes with augmented catalytic efficiency would benefit green chemistry goals for the pharmaceutical industry in manufacturing specialized steroid chemicals [
85].
Mammalian steroid isomerases have multifunctional activity and a more complex structure than bacterial KSI enzymes [
86]. Although crucial in all mammals, their structure-function connection has not been as extensively examined as in the case of the model enzyme KSI from bacteria. Thus, we shall use the best predicted KSI rate constants for bacterial KSI [
31] that agree well with those reported earlier by [
80,
87].
Our first task was a wide exploration of possible system states when noise is introduced into each of the eight rate constants for the 4-state kinetic scheme (
Figure 15). Our FORTRAN simulation kept the concentrations of substrates and products fixed at their initial values (
Table 2, second column: [S] = 10
-4 M, [P] = 5∙10
-5 M). Nevertheless, due to random changes in all rate constants, the force changed in the range 0.72 < X
tot/RT < 17.17. Repeated runs produced an identical output. The third best efficiency value from the (1.71∙10
4, 1.66∙10
9) point reveals that 5-fold higher efficiency can be achieved when corresponding entropy production is 10-fold smaller than their experimental values. That is a rare case when the choice of rate constants results in high catalytic activity despite low dissipation for the KSI enzyme.
As for the case of triosephosphate isomerase, the perfect efficiency-dissipation proportionality followed after the no-change requirement in the equilibrium constants for all catalytic steps (not shown). When noise is called only once, a nearly perfect linear fit survives for the efficiency-dissipation dependence, no matter how many rate constants ki are multiplied with the normal noise function (not shown). The consequence of fixed equilibrium constants Ki is constant overall force, too.
Regular dependence of enzyme efficiency on overall dissipation follows when noise is introduced only into one or two kinetic constants without fixed K
i requirements (
Figure 16 and
Figure 17). However, that dependence is very different if the overall force X
tot/RT is allowed to vary too (
Figure 16) and when overall force is kept at the constant initial value X
tot/RT = 8.426 (
Figure 17, see
Table 2).
Figure 17 confirms the observation from
Figure 9 that the maximum in overall entropy production exists when variations in K1 and K4 equilibrium constants are introduced and fixed overall force is maintained in all steps of FORTRAN simulation. Total entropy production is maximal in the point (1.3∙10
5, 4.7∙10
8)(
Figure 17). The corresponding optimal efficiency is about 50% higher than the observed value of. 3.02∙10
8 M
-1s
-1. Still, the point with the highest efficiency (1.96∙10
4, 8.15∙10
8) corresponds to 5.9 times smaller dissipation than the value 1.16∙10
5 s
-1 calculated from the experimental data. That is another rare case when randomly chosen equilibrium constants within imposed restriction (constant overall force) resulted in a high catalytic efficiency despite low dissipation for the KSI enzyme.
Agent-based modeling extended and confirmed the simulation results for the KSI kinetics using the FORTRAN software. The advantage of NetLogo simulation (
Figure 18) is that random changes in the concentrations of substrates, products, free enzymes, and enzyme complexes are allowed. Typical Michaelis-Menten kinetics for concentration changes, which we described for the NetLogo simulation of the TPI kinetics, is also seen for the KSI kinetics (not shown). Initial concentrations were [E]
free = 5 μM, [S]
free = 100 μM, [P]
free = 50 μM. Final foncentrations at the 6977th tick were [E]
free = 4 μM, [S]
free = 95 μM, [P]
free = 54 μM, [ES] = 0.3 μM, [EX] = 0.4 μM, [EP] = 0.3 μM. The mass conservation conditions [E
tot] = [E
free] + [ES] + [EX] + [EP] and [ligands] = [S]
free + [P]
free + [ES] + [EX] + [EP] were satisfied through all time jumps (ticks). Freedom to change equilibrium constants in each tick enabled the exploration of a wide range for overall force (1.1 < X
tot/RT <16.8), catalytic efficiency, and overall dissipation. The best pair of dissipation-efficiency values (4.6∙10
5, 2.6∙10
9) corresponded to approximately 4-fold higher dissipation and almost 10-fold higher efficiency in the comparison with values calculated from observed data (1.2∙10
5, 3.0∙10
8).
CA I, CAII, CAII-T200H chapter (also 4-state enzymes)
Carbon dioxide conversion into biomass is essential for the survival and spreading of life in all terrestrial environments. Carbon sequestration is also crucial for the survival of our carbon dioxide-producing civilization, which is unfortunately addicted to fossil fuels burning and breaking all life-supporting balances the biosphere has developed through eons. Nature developed multiple means and different organic structures for the fast conversion of carbon dioxide to bicarbonate – the first step toward carbon fixation. Carbonic anhydrases (CAs) are universal enzymes responsible for that process in all three life domains: Bacteria, Archaea, and Eukarya [
88]. With rare exceptions [
89], CAs are metalloenzymes containing the metal ion (most often zinc) in their central active-site cavity. From their discovery in red blood cells in 1932, the scientific interest in CAs continued to grow, as seen from the abundance of more than 900 solved CA structures deposited in the Protein Data Bank [
90].
The spontaneous reaction of CO
2 with water can produce bicarbonate
, but that reaction is too slow to support respiration [
91,
92] and other biological processes catalyzed by different CAs. Eight unrelated families of carbonic anhydrase (CA) enzymes represent different ways nature performed the feat of fast catalytic inter-conversion between carbon dioxide and carbonic oxide [
93], reaching the catalytic turnover of 1 μs
-1 or even higher [
94]. There is little or no sequence homology among CA families α, β, γ, δ, ζ, η, θ, and ι [
89,
95]. Molecular biologists concluded that convergent biological evolution performed the spectacular function-enhancing feat at least seven times because different CAs evolved to perform an identical function [
93,
96,
97].
Mammals possess 16 different CA isoenzymes from the alpha class family [
98]. All are metalloenzymes, with the Zn II hydride located at the enzyme center anchored by three histidines. CA isoforms are involved in a variety of physiological functions. Human CA isoforms are well-recognized drug targets for designing isoenzyme-specific inhibitors [
99,
100] to help fight glaucoma, epilepsy, obesity, cancer, and other diseases. Also, human CA II is one of the most efficient known enzymes. Its calculated catalytic efficiency from experimental data is 1.5∙10
8 M
-1s
-1 [
100]. Earlier efficiency calculations also positioned CA among “perfect” enzymes working close to the diffusion limit [
14,
56].
Genetic defects of specific CA isoforms can cause osteopetrosis, cerebral calcifications, retinal problems, hyperammonemia, hyperchlorhidrosis, neurodegenerative and other metabolic diseases [
101], which is a good enough reason to look for CA activators [
98] or other means for increasing the activity of these isoforms. Memory enhancement can be achieved through CA activation [
102]. It opens the possibility for targeted improvement of brain CA performance to enhance cognition and slow the aging process [
98,
103]. Some CA mutants can accelerate the proton transfer, the rate-limiting step for CA turnover [
94]. Another reason for increasing CA activity is the urgent need for green ways of industrial CO
2 sequestration [
104] we mentioned above.
This subsection deals with the theoretical possibilities for catalytic efficiency improvements of human CAs I, II, and the T200H variant of CAII with His200 replacing Thr200 [
105]. There may be better models than the four-state kinetic model for reversible Michaelis-Menten-type kinetics (
Figure 19). Still, it is based on the publication [
105] that contains all microscopic rate constants needed to calculate and compare the enzyme’s performance with associated dissipation. Referent (initial) state values can be found in
Table 3.
The FORTRAN simulation of noisy CA I kinetics did not change any of the initial concentrations (
Table 3), and it still found in the 246th step the dissipation-efficiency point (5.14∙10
5, 1.12∙10
8) with the 4.5-times higher catalytic efficiency from the calculated value based on the observed kinetic data. The corresponding overall force was positive (X
tot/RT = 5.0) and closer to the upper end of the force range (X
tot/RT = 8.3). However, the substantial efficiency increase (4.5-fold) was „paid for “ with the 18 times higher overall dissipation. Closer inspection of the performance parameters from the 246th computational step (concerning observed initial values) revealed the 6.3-fold increase in the turnover number and 2.8-fold increase in the overall force as the main reason for the improved efficiency.
Figure 20.
The FORTRAN simulation of the relationship between catalytic efficiency and overall dissipation for the carbonic anhydrase I, when each of eight rate constants k
i is multiplied with independently called normal noise function gi without shift (see Methods, eq (25)). Only positive gi values were allowed while the program ran through two thousand steps. There was no change in the substrate, product, and buffer concentration from their initial values (
Table 3).
Figure 20.
The FORTRAN simulation of the relationship between catalytic efficiency and overall dissipation for the carbonic anhydrase I, when each of eight rate constants k
i is multiplied with independently called normal noise function gi without shift (see Methods, eq (25)). Only positive gi values were allowed while the program ran through two thousand steps. There was no change in the substrate, product, and buffer concentration from their initial values (
Table 3).
The NetLogo simulation of noisy CA I kinetics (
Figure 21) slightly changed initial substrate and product concentrations (
Table 3). The 6-fold efficiency increase point, which we found halfway through the simulation, was „paid-for“ with the 22-fold dissipation increase. That quasi-steady state corresponded to about a 3-fold increase in the k
cat and overall force.
Human red cell isoenzyme CAII is superior to CAI when their catalytic efficiencies are compared [
105] (see
Table 3). Thus, simulations for the CAII will have the advantage of starting from a better initial state. Here, we show only the NetLogo simulation (
Figure 22). The CAII mutant T200H, constructed by Behravan et al. [
105], was an attempt to find the single amino acid substitution that would lead toward catalytic parameters of CAI. The NetLogo simulation (
Figure 23) indicates the evolutionary potential for improving the performance of CA-T200H as being indeed between CAI and CAII but closer to CAII (
Figure 24).
We compared enzyme performance and associated dissipation for three CAs isoenzymes: CAI, CAII, and the T200H mutant of CAII (
Figure 24). Values calculated from observed rate constants [
105] are confined near the origin of that figure, while the best-simulated values are expanded in the order CAII > CAII T200H > CAI. Improved catalytic efficiency is associated with increased dissipation in the same order.
PC1-β-Lactamase
When we maintain the same restrictions of unchanged initial values for the equilibrium constants (except for the changes in the substrate and product concentrations), identical normal noise introduction in all forward kinetic constants leads to only slight changes in the nearly perfect proportionality between catalytic efficiency and dissipation (
Figure 25). Besides noise, an additional reason for changes in k
1 and k
6 is the decrease in the free substrate concentration and an increase in the free product concentration during enzyme cycling scheme E+S ↔ ES ↔ EP ↔ E+P. It produces a slight decrease in the first equilibrium constant K
1 = k
1*∙[S]/k
2 and an increase in the third equilibrium constant K
3 = k
5/(k
6*∙[P]. Increased product concentration is the main reason for the gradual force decrease, from initial X
tot/RT = 11.4 to final X
tot/RT = 10.8 after 5168 ticks of the NetLogo simulation. At the 381st NetLogo simulation tick, we found the best efficiency value of k
cat/K
M = 4.18∙10
7 M
−1s
−1, which corresponds to forward rate constants k
1 = 1.35∙10
5 s
-1, k
3 = 717 s
-1, k
5 = 398 s
-1, catalytic constant k
cat = 252 s
-1, and dissipation/RT = 2823.6 s
-1. That is the same 4.1-fold increase for all of these parameters with respect to their observed values (see
Table 4).
The next goal is to look for limits to the evolvability of PC1 β-lactamase subject to variable rate constants k
1 and k
2 in the first catalytic step (association-dissociation between the free substrate and free enzyme: E+S ↔ ES).
Figure 26 illustrates how the multiplication of k
1 and k
2 with the, respectively, Box-Muller normal noise functions named g
1 and g
2 can find a quasi-steady state with 6.5 times higher catalytic efficiency and merely 1.2 times higher dissipation in comparison with those values calculated from experiments (
Table 4). That is a significantly better result than all previous optimizations [
30] based on the requirement of maximal partial entropy production in the proton transfer catalytic steps 2 (ES ↔ EP) and 3 (EP ↔ E+P). For instance, joint optimizations of both catalytic steps for maximal transitional entropy production in these steps find about 2-fold higher efficiency, which is "paid for" by the 183 times higher dissipation. The maximum total entropy production requirement combined with the obligatory K
+ = k
1∙k
3∙k
5 = const constraint leads to 333 times lower catalytic efficiency (
Figure 26, and [
30]).
No further gain in enzyme efficiency follows after normal noise is independently introduced in four or all six rate constants. The maximal k
cat/K
M was in a narrow range from 5.9∙10
7 to 6.2∙10
7 M
-1s
-1. Thus, we shall keep the best result from
Figure 26 (6.5∙10
7 M
-1s
-1) to compare evolutionary potential with other enzymes.
Lac1-β-Lactamase
For the case of the Lac1-β-lactamase kinetics, we explored several options for the independent noise introduction in each kinetic constant from chosen pairs (
Figure 29,
Figure 30 and
Figure 31). It turned out that the k
1, k
2 pair is the best choice because it led to the catalytic efficiency value of 1.25∙10
8 M
-1s
-1, which is also inside the range 10
8–10
10 M
-1s
-1 for diffusion-limited enzyme reactions [
63].
When normal noise is introduced only once in the forward rate constants k
1, k
3, and k
5 with the proviso that the equilibrium constants K
1* = k
1*/k
2, K
2 = k
3/k
4, and K
3* = k
5/k
6* do not change from their observed values, the perfect proportionality follows for all efficiency-dissipation pairs (
Figure 32). Due to higher initial product concentration (see
Table 4), the constraints K
1* = const1, K
2 = const2, K
3* = const3 are almost equivalent to the requirement that initial equilibrium constants K
1, K
2, K
3 never change during the NetLogo simulation for the Lac-1 β-lactamase kinetics. Since (k
cat/K
M)/dissipation expression depends only on equilibrium constants and the ratios of rate constants (see Appendix), there is no reason for the slope change in the efficiency-dissipation dependence (
Figure 32). The best efficiency value of 9.68∙10
7 M
-1s
-1 is close to the diffusion-limit range's lower end (10
8 M
-1s
-1). It was reached at the 2182nd tick for the X
tot/RT = 8.23 and 3.7-fold higher turnover k
cat = 7086 s
-1.
To summarize, we have seen nearly perfect kinetic-thermodynamic proportionality for the PC1 (
Figure 25), RTEM (
Figure 27), and Lac-1 (
Figure 32). Corresponding FORTRAN programs confirmed it for all three β-lactamases (
Figure 42). We also confirmed the efficiency-dissipation proportionality for the triosephosphate isomerase kinetics (
Figure 4,
Figure 6,
Figure 7,
Figure 10,
Figure 12 and
Figure 13). It will likely hold whenever the no-change condition is imposed for the equilibrium constants in all catalytic steps (see Appendix for more details). The capture-release initial step leads to different relationships when the no-change condition is imposed on all first and second-order rate constants except k
1 and k
2 (
Figure 26,
Figure 28 and
Figure 31 for β-lactamases). A fast enzyme efficiency increase can then occur for limited dissipation. The potential for the exponential-like efficiency increase is likely to be a general phenomenon for all Michaelis-Menten enzymes after lowering the activation barrier for the E+S→ES transition and increasing the activation barrier for the reverse ES→E+S transition.
Figure 33 illustrates the relationship between evolutionary distance and overall entropy production for PC1, RTEM, and Lac-1 lactamase. We found blue points and corresponding fit lines (black) from the simulation of experimental data [
30,
106]. The dissipations associated with the red points (and red fit line) are from the best catalytic efficiency points in
Figure 26,
Figure 28 and
Figure 31. The dissipation increased in an almost linear manner for more evolved β-lactamases. Noise introduction and searching for the highest enzyme efficiency confirmed the proportionality between the time passage (evolutionary distance) and overall entropy production. It is to be expected if we can consider the cumulative entropy production as a surrogate for time passage.
Glucose isomerase
Glucose isomerase (GI abbreviation) fulfills nutritional requirements mainly in bacteria [
124]. It is also known as xylose isomerase because GI reversibly isomerizes D-glucose and D-xylose to D-fructose and D-xylulose, respectively. Glucose to-fructose conversion is rather inefficient but critical for the commercial production of high-fructose corn syrup (HFCS) due to the specificity (the absence of non-metabolizable or toxic side products) and mild ambient conditions [
125]. Together with other industrial applications through decades, such as bioethanol production [
126], the GI maintained a high market share of the food industry among other industrial enzymes despite its inherently low activity [
124,
127]. That is why there are so many research attempts to use molecular engineering to improve GI performance for different applications [
125,
126,
128,
129,
130]. Besides academic interest, that is also the reason for theoretical research devoted to improving GI performance when kinetic and thermodynamic limits are taken into account.
Considering previous examples for other enzymes, the best option is to initiate the research with all the microscopic rate constants inferred from the observed data. Unfortunately, the best such example [
27] is for the GI preparation from
Streptomyces murinus with very low measured activity. Nevertheless, the principles we employed in this section to significantly improve the GI performance may be applicable for predicting activity gains of the most promising GI variants for green industry applications. We also wanted to test our hypothesis about catalytic efficiency proportionality to dissipation by using the example of an inefficient enzyme working close to the thermodynamic equilibrium. The drawback is the restriction to the two-state model (
Figure 1a), that is, the reversible Briggs-Haldane mechanism used in early and recent proposals for the kinetic mechanism of immobilized GI [
27,
131,
132,
133,
134]. In the pseudo-steady state, the solution is the Michaelis-Menten equation and corresponding performance parameters k
cat, K
M, and k
cat/K
M. In our notation (
Figure 1a) k
3 = k
cat and K
M = (k
2 + k
3)/k
1s , where k
1s is the second order rate constant. Unusual experimental conditions used to determine these kinetic parameters include an elevated temperature (65
◦C) besides the immobilization of enzymes. These conditions are not responsible for observed low k
cat and k
cat/K
M values. If anything, they slightly increase the performance parameters.
Table 7.
Initial values of microscopic rate constants from experimental and estimated data [
27,
28] for GI to which we added our calculations of other initial kinetic and thermodynamic parameters in the case of
Streptomyces murinus catalyzed conversion of glucose (substrate) to fructose (product) at 65
◦C.
Table 7.
Initial values of microscopic rate constants from experimental and estimated data [
27,
28] for GI to which we added our calculations of other initial kinetic and thermodynamic parameters in the case of
Streptomyces murinus catalyzed conversion of glucose (substrate) to fructose (product) at 65
◦C.
Rate constants |
Observed values [27] |
Calculated values [28] |
k1* |
3.8 M−1min−1
|
0.063 M−1s−1
|
k2
|
1.23 min−1
|
0.021 s−1
|
k3
|
1.75 min−1
|
0.029 s−1
|
k4* |
4.9 M−1min−1
|
0.082 M−1s−1
|
Other relevant parameters |
|
Initial values (this paper) |
[S] |
|
2.0 M |
[P] |
|
0.2 M |
[E] |
|
0.01 M |
k1
|
|
0.126 s−1
|
k4
|
|
0.0164 s−1
|
kcat
|
|
0.029 s-1
|
KM
|
|
0.794 M |
kcat/KM
|
|
0.0365 M−1s−1
|
Keqtot
|
|
10.61 |
Xtot/RT |
|
2.31 |
|
|
Initial value (this paper) |
P |
|
0.0392 s-1
|
Figure 39 NetLogo simulation with noise independently introduced in all rate constants illustrates the absence of a strong proportionality relationship between catalytic efficiency and entropy production for an inefficient enzyme, such as glucose isomerase. The best point with the coordinates (0.21, 0.21) was found at the 1715th tick. It is associated with the X
tot/RT = 4.7, about 2.5-fold higher k
1, 9-fold smaller k
2, 24-fold higher equilibrium constant K
1, and approximately 11-fold higher partial entropy production
P1. Thus, the association and dissociation of substrate with enzyme should be shifted toward complex ES formation for the significant 5.8-fold increase in enzyme efficiency and 2.3-fold increase in the turnover number.
Imposing some constraints on the system can recover the efficiency-dissipation relationship. For instance, rate constants k1 and k2 can be independently multiplied with the normal noise without changes in the rate constants k3 and k4 other than those caused by increased product concentration. The correlation R2 then jumps to 0.881 for kcat/KM dependence on dissipation. However, NetLogo simulation goes through restricted search space and finds somewhat lower values for the best efficiency (not shown).
We also performed the FORTRAN simulations to verify that different software and ways for noise introduction can still produce an approximately linear response in the catalytic efficiency to the dissipation expressed as forcing XJ combination of external force and internal current (
Figure 40). Normal noise was called only once and used to multiply all four rate constants. The best catalytic efficiency of 0.18 M
-1s
-1 was comparable to the best result for the NetLogo simulations with noise. It was also considerably better than the 0.0215 M
-1s
-1 catalytic efficiency, easily calculated from the Dobovišek et al. optimization [
28]. Incidentally, Dobovišek's result was obtained for the positive force X
tot/RT = 0.51 when the net flow is in the forward direction, and our two-state expressions for k
cat and K
M (see Methods) are appropriate to use for the calculation of k
cat/K
M. It emerged due to the special nature of the quasi-steady state these authors obtained after imposing the no-change constraint for the product of forward rate constants:
= const.
Normal noise with shift +1 (see Methods) was called four times in the next FORTRAN simulation so that each rate constant was multiplied with its own Box-Muller transform (
Figure 41). The best catalytic efficiency of 0.226 M
-1s
-1 result was found for lower overall dissipation in the RT units (0.12 s
-1) compared to the previous NetLogo simulation (
Figure 39).
Dissipation, Evolution, and Catalytic Power of Enzymes
The evolution of the Universe can be described as the universal evolution. It created time, space, and myriad beautiful objects such as galaxies, stars, planets, and living beings [
135]. The invisible but not less important product of universal evolution is increased entropy and entropy production. Through time passage, the product of absolute temperature and entropy production (the dissipation) kept rising from the mysterious initial singularity, which must have had very low entropy, extremely high temperature, and zero entropy production. A new phase of universal evolution started with the appearance of objects that can be associated with the huge jump in entropy production. The first class of such objects are black holes, astronomical objects originating after the explosive death of some massive stars. The second class of objects originated in an aqueous environment endowed with rich chemistry as the first living cells. The same volume of some bacterial cells, mitochondria, or chloroplast produces many orders of magnitude higher dissipation than an equivalent average volume of a sun-like star despite the star's much higher temperature [
136,
137]. The specific variety of complex life and mineralogy we enjoy here on Earth is not likely to exist anywhere else in the Universe [
138]. Thus, we should protect it, study it, and, if possible, understand it as a natural consequence of universal evolution.
It is hard to imagine life without housekeeping enzymes. The simplest and the most successful description of how many such enzymes work is generalized Michaelis-Menten kinetics [
36,
139,
140,
141]. The increased complexity of life through eons required means for increasing the catalytic efficiency of such enzymes. According to our knowledge, the scientific literature has never seriously considered how catalytic efficiency is related to an enzyme's entropy production. Banerjee and Bhattacharyya's finding [
18] that the more efficient enzyme involves higher total dissipation is in accord with the results presented in this paper. The finding is based only on three pairs of dissipation-efficiency values for a single enzyme (β-galactosidase). Still, it is gratifying that their different method for calculating overall entropy production produced the same result (2553 s
-1) as Terrel Hill's approach [
34] we used throughout this paper (see the last row from
Table 5). How changes in dissipation can lead to an increased catalytic efficiency was not the main interest of these authors.
Martyushev and Seleznev [
142] anticipated a fruitful connection between optimal kinetics parameters and entropy production for strongly nonequilibrium processes. However, it is surprising that the relationship between overall dissipation and frequently measured specificity constant k
cat/K
M was never thoroughly examined. These two parameters connect laboratory biochemistry with the fundamental thermodynamics of nonequilibrium processes. Wofenden and coworkers found with other authors that nothing makes a sharper distinction between life and non-life from the massive jump in the catalytic power, which enzymes show when the speed and specificity of the reaction they catalyze is compared to an equivalent reaction in the presence of inorganic catalysts [
10,
12,
13,
14,
15,
16,
143,
144,
145].
Great strides have been made in uncovering the 3D structure of proteins with enzymatic activity [
146]. Still, dynamic changes essential for understanding their catalytic activity are difficult to trace structurally [
147]. Structural studies did not help as much as we hoped in answering how enzymes work [
4,
148]. Initial expectations that entropy changes are the most important contribution to the stabilization of the transition state [
149,
150] were not confirmed [
151], nor the expectations that determining an enzyme's mechanism of action is enough to understand its catalytic power [
152].
Despite all the complexity of the enzymes we mentioned here, they are simpler than organisms for investigating evolution [
153]. That "cleaner" opportunity brings us back to our approach to the mystery of how enzymes work. Evolution in physics is firmly connected to entropy production. At the same time, the marvelous biochemistry of enzymes is tied to the evolutionary enhancements of enzymes' catalytic rate (up to 10
26-fold, according to Edwards et al. [
16]). Corresponding catalytic proficiency for alkylsulfatase is an astronomical number: (k
cat/K
M)/k
uncat = 10
29 M
-1. The mystery lies in the vast increase of enzyme efficiency k
cat/K
M during pre-biological and biological evolution. Since an increase in entropy production speeds up the physical evolution of any nonequilibrium system (it undergoes faster relaxation from the initial far-from-equilibrium state), we can assume the connection with the evolution of catalytic efficiency.
The initial impetus for the evolution must be the relaxation tendency, that is, the appearance of fluxes in the system, which spontaneously dissipate external force gradients. All products of forces and corresponding fluxes together are nothing else but overall dissipation. Thus, the dissipation is a dynamic phenomenon including external driving forces, fluxes through the system, and output fluxes toward the system's exterior. A system in a quasi-steady nonequilibrium state is an almost perfect converter. It converts free energy influx into internal dissipative fluxes, and external entropy increase, which is usually detected as the heat flux. A small percentage of free energy influx can be converted into transient free energy storage within the system. For biological systems in nearly stationary, far-from-equilibrium states, a hierarchy of internal free-energy conversions can lead to dissipative adaptation [
154]. A slight change in the system's properties increases its ability to absorb free energy from the environment. At the same time, the system must produce a slight increase in its dissipative output to remain in a stable quasi-steady state. After long enough free-energy accumulation, some systems acquire the essential biological property of self-replication.
A large outflow of entropy from the growing cells is the major thermodynamic process [
155]. Only a small portion of available free energy is used by cells for synthetic and mechanistic goals. For instance, the free energy converted into chemical bonds is a minor contribution compared to the free energy change from catabolism. Still, an almost perfect correlation exists between the total heat released and the amount of dry mass grown or the total amount of oxygen consumed during the aerobic growth of a yeast culture [
156]. Theoretical studies also concluded that a higher maximal growth rate would be achieved by replicating a system capable of producing more heat [
154,
157]. Thus, higher entropy production can be an advantage during the evolution of organisms. As a rule, total entropy production reaches its maximum value before a slow decrease when microorganisms are fully supplied with free-energy sources and engaged in vigorous growth during their short-term evolution in batch experiments. This pattern is recapitulated in the life of every individual organism. Metabolic heat production per surface area reaches the maximal value early, with a subsequent decline over the lifetime [
158].
The metabolic heat production is due to enzymes. Using the microcalorimetry method, Sica et al. found in 1987 [
159] the proportionality between enzyme activity k
cat and observed heat flow. Thermal power was directly proportional to the reaction rate for dihydrofolate reductase (EC 1.5.1.3). Todd and Gomez [
160] extended that observation to representative enzymes from each EC classification (a total of 11 different enzymes), assuming the validity of the Michaelis-Menten equation [
21]. Isothermal titration microcalorimetry can be used for direct, nondestructive, and precise reaction rate measurement. Todd and Gomez [
160] found a reasonably good agreement between kinetic parameters k
cat, K
M, and k
cat/K
M assayed colorimetrically and with other methods. Riedel et al. [
161] confirmed the agreement of calorimetric and kinetic parameters k
cat and K
M for catalase, urease, alkaline phosphatase, and triosephosphate isomerase. These authors also studied the feedback effect of released heat in enzyme-catalyzed reactions. Surprisingly, some enzymes asymmetrically release enough heat to increase their diffusion in the presence of substrates.
It is the chemotactic behavior. The enzyme preferentially diffuses towards higher substrate gradients, thus increasing its activity due to increased dissipation. Therefore, nature found a way to channel dissipated and "useless "energy released by enzyme catalysis towards helpful purposes, such as directed motion. For instance, directional water ejection of hexokinase classifies it as a biological pump despite being a cytoplasmic enzyme [
162,
163]. Hexokinase-2 (HK-II) in humans has been described as a gatekeeper between life and death [
2] because it integrates glycolysis with oxidative phosphorylation in healthy cells, contributing to stronger respiratory control. HK-II chemotaxis, due to the entropically favored expulsion of water outside the enzyme pocket, is probably instrumental in its preferential access to ATP generated by mitochondrion.
Assuming that KM does not change for a chosen enzyme, the observed proportionality between the enzyme's entropy production and the turnover number kcat implies a linear increase in catalytic efficiency kcat/KM with dissipation. When noise is present in rate constants, approximate KM constancy will still hold for no changes in k3/k1 ratio (2-state equation (17)), or k5/k1 and k5/k4 ratio (3-state equation (20)), or k3/k1, k3/k5, and k5/k7 ratio (4-state equation (23)) when equilibrium constants do not change.
Through this paper, we could have presented both k
cat and k
cat/K
M proportionality with overall entropy production for each enzyme. One example of the proportionality between k
cat and dissipation is for three β-lactamases (
Figure 42). Besides k
cat to dissipation/RT proportionality (in the units of inverse seconds), that figure also illustrates nearly linear connection between the evolutionary distances of PC1 (1.19), RTEM (1.44), and Lac1 (1.60) lactamase and either k
cat or overall dissipation (see
Figure 33, [
30], and [
106] for the evolutionary distances we put in parenthesis). We obtained the same result after comparing experimental results for kinetic and thermodynamic parameters of the A-class β-lactamases and after looking for the maximal partial dissipation in the rate-limiting steps [
2,
106] In these and other publications [
30], we stressed that the optimization for the turnover numbers should be based on the physical principle of maximum transitional entropy production, not on the uncritical acceptance of the maximal catalytic efficiency or maximal catalytic constant as the selection or optimization criterion.
There was no need in the present study to make a priory assumption of either physical or biological principle reigning supreme. We only required some mechanism for reasonable variations in the microscopic rate constants. A crowded cellular milieu and unavoidable errors in translation and transcription offer several such means for noise introduction in kinetic parameters. Stochastic fluctuations are always present and are relevant for applying the Michaelis-Menten type kinetics inside cells [
164]. Our simulations are, admittedly, a crude and artificial way of considering the noise. Better methods for dealing with physical and biological noise sources are undoubtedly possible. However, we were primarily interested in whether different means of noise introduction can uncover regular relationships between the most important thermodynamic and kinetic parameters for highly active enzymes that work arbitrarily far from equilibrium. Making use of thermal and non-thermal noise through stochastic fluctuations and dynamic disorder [
165,
166] may well have been beneficial during biological evolution [
167,
168,
169,
170].
Standard evolutionary theory [
171] has a simple answer to the question of how new variations can arise: random mutations and natural selection, ensuring the adaptation of the organisms to their environment. Thus, a particular noise class (chance genetic changes) is adapted to provide a better fit among organisms and environments. That view has been extended recently by taking into account the physicochemical evolutionary driving forces [
172], including the maximization of dissipation [
1].
Figure 42.
FORTRAN simulations for the k
cat dependence on overall dissipation in the case of three β-lactamases (PC1, RTEM, and Lac1). Each forward rate constant was multiplied with the identical normal noise function while corresponding backward rate constants were determined from the no-change requirement to equilibrium constants.
Table 4 parameters were used for each enzyme. Concentrations were not allowed to change from their initial (observed) values. As for other figures, the P
tot label at the x-axis is the dissipation/RT in inverse seconds. The figure illustrates the proportional increase or decrease of the turnover number with dissipation from observed (calculated) points for PC1 (689, 61), RTEM (6757, 975) and Lac-1 (14526, 1905) (see
Table 4 and Juretić et al., 2019 [
106]). The highest points have coordinates (3035, 268) for PC1, (3∙10
4, 4303) for RTEM, and (6.4∙10
4, 8394) for Lac-1. As for catalytic efficiencies (
Figure 33), both observed, and the highest points (dissipation, k
cat) are nearly proportional to the evolutionary distance from the putative common ancestor in the order PC1 (1.19) < RTEM (1.44) < Lac1 (1.60).
Figure 42.
FORTRAN simulations for the k
cat dependence on overall dissipation in the case of three β-lactamases (PC1, RTEM, and Lac1). Each forward rate constant was multiplied with the identical normal noise function while corresponding backward rate constants were determined from the no-change requirement to equilibrium constants.
Table 4 parameters were used for each enzyme. Concentrations were not allowed to change from their initial (observed) values. As for other figures, the P
tot label at the x-axis is the dissipation/RT in inverse seconds. The figure illustrates the proportional increase or decrease of the turnover number with dissipation from observed (calculated) points for PC1 (689, 61), RTEM (6757, 975) and Lac-1 (14526, 1905) (see
Table 4 and Juretić et al., 2019 [
106]). The highest points have coordinates (3035, 268) for PC1, (3∙10
4, 4303) for RTEM, and (6.4∙10
4, 8394) for Lac-1. As for catalytic efficiencies (
Figure 33), both observed, and the highest points (dissipation, k
cat) are nearly proportional to the evolutionary distance from the putative common ancestor in the order PC1 (1.19) < RTEM (1.44) < Lac1 (1.60).
Computational Improvements of the Catalytic Power for Specific Enzymes
The catalytic power of enzymes is measured as k
cat or k
cat/K
M. Some authors did not recommend using k
cat/K
M as an index for comparing the catalytic effectiveness of enzymes [
173]. The majority consensus is that k
cat/K
M is the appropriate measure for the specificity of noncooperative Michaelis-Menten enzymes [
24,
174]. In rare cases when all microscopic rate constants have been determined [
175], k
cat and k
cat/K
M can be connected to partial and total entropy production when an enzyme reversibly cycles through all of its functionally important conformations (this work and [
2]). Moreover, after variations of rate constants around their observed values, we can analyze optimal rate constants k
i and dissipations associated with the highest performance parameters k
cat and k
cat/K
M. What are, if any, common features of the states with the highest enzyme efficiency, and how do the thermodynamic and kinetic parameters of these states differ from the same values calculated or inferred from the experimental data?
Table 9 helps deal with that question. Our choice in this paper was to examine the best k
cat/K
M values for corresponding k
i, partial, and total dissipation.
Table 9 gives the partial entropy production in the first forward catalytic step because it exhibited the highest increase regarding the observed value. That is the consequence of increased forward rate constant k
1 and decreased backward constant k
2 for the substrate-to-enzyme association and dissociation.
In the case of triosephosphate isomerase (TPI), there was a 1454-fold increase of the partial entropy production
P1 for the E+S↔ES transition, which became the 42% instead of 6% contribution to the total dissipation (NetLogo result from
Figure 14). For other enzymes and software simulations in the presence of noise,
P1 increased one to two orders of magnitude, and its percentage also increased for the best enzyme efficiency results.
The single exception is the carbonic anhydrase. For CAI, CAII, and the T200H CAII mutant, an absolute increase in
P1 was not accompanied by an increase in its percentage. A possible reason is a different kinetic scheme for the CA enzyme (
Figure 19) and an inadequacy of the standard k
cat and k
cat/K
M expression (equations (22)-(24) ) for that scheme. Krishnamurthy et al. [
92] (
Table 1 from [
92]) compared all known CA enzymes for their
and
values for the catalytic hydration of CO
2 and the dehydration of bicarbonate:
That is the first half-reaction. The buffer is considered the second substrate in the overall two-substrate ping-pong reaction, which recovers free enzymes.
The best efficiency-fold improvement is seen for β-galactosidase, which also reaches the highest efficiency-to-dissipation fold ratio (
Table 8). However, ketosteroid isomerase has the best evolutionary potential in our simulations. That can be connected to the two proton transfer reactions catalyzed by KSI [
176,
177,
178] and a powerful electric field [
152]. Electric field catalysis needs a strong and correctly oriented field. The measured field of 1.44∙10
10 V/m is enough to account for 72% of the total acceleration rate [
152]. The transient appearance of billion volts per meter electric field strength in the interior of active proton-shuffling enzymes frequently speeds up catalysis [
2]. The isomerization of 5-androstene-3,17-dione in solution through the same mechanism utilized by KSI is slow. That is why KSI catalytic proficiency is so high. As mentioned in the main text, Radzicka and Wolfenden [
12] estimated it as 1.8∙10
15 M
-1 based on k
uncat = 6∙10
-7 s
-1.
Interestingly, the efficiency-fold improvement (
Table 8) is similar for the best (KSI) and worst enzyme (glucose isomerase). The k
cat = 0.029 s
-1 (experimental) and k
cat = 0.031 s
-1 or 0.068 s
-1 (optimal) values for the k
3 in
Table 9 GI results are two orders of magnitude smaller than the turnover numbers 2 s
-1 and 11 s
-1 reported in the literature [
128,
179]. The Converti et al. [
27] data we used to initiate simulations pertain to weakly active GI working close to thermodynamic equilibrium. Nevertheless, our method for the theoretical increase in the catalytic activity is robust enough to ensure its close to 10-fold increase (from 0.0365 to 0.2262 M
-1s
-1,
Figure 41).
Possible Benefits of Considering Unanswered Questions
Most enzymes did not use their potential to evolve higher catalytic efficiencies due to the absence of selection pressure to maximize it for individual enzymes [
19]. When metabolic demand existed, the superstars of enzyme evolution developed, often named perfect enzymes [
175]. Our simulations suggested the theoretical possibility of increasing the k
cat/K
M of either moderately efficient or perfect enzymes. In practice, more than one amino acid substitution is needed to improve the performance parameters. Several orders of magnitude improvement typically requires at least 5 to 10 beneficial mutations [
180].
Living far from equilibrium is an essential asymmetry of present-day life [
181,
182]. Total entropy production is a convenient measure of how far the system is removed from thermodynamic equilibrium. A certain distance from equilibrium must be maintained to increase the beneficial catalytic efficiency after increased dissipation by whatever means for enzymes we studied in this paper. The plausible inference is that some abiotic driving forces, such as proton gradients in alkaline hydrothermal vents, must have operated to maintain far-from-equilibrium situations and high entropy production during life emergence on Earth. According to that assumption, bioenergetics and vectorial biochemistry are older than the genetic code and the first universal common ancestor [
2]. It enabled enzyme-less and cell-less synthesis of amino acids, sugars, nucleotides, and lipids. Nonlinearity and far-from-equilibrium conditions are two requirements for driving proto-metabolism toward autocatalysis and self-organization. Accelerated accumulation of organic molecules followed in the presence of long-lived abiotic protonmotive force to jump-start the development of life [
183]. The efficiency of organic synthesis with proto-enzymes was surely low compared to present-day enzymatic catalysis. However, such self-reinforcing reactions increased the efficiency of dissipating available free-energy gradients. The present-day connection between dissipation and catalytic efficiency we studied in this paper is thus likely to reflect the linkage between the higher dissipation potential and accelerated synthesis of ever more complex organic compounds, which was already present at the origin of life. Entropy production increases faster due to the enzyme's activity, albeit in the microscopic world.
Within biology, we cannot find the answer to why dissipation was crucial for the emergence of life, as it is essential for the present-day catalytic efficiency of uni-uni enzymes. Can entropy production have an autocatalytic role too? Namely, did increased entropy production promote the selection of the organic structures capable of increasing entropy production? That question has yet to be answered in biology or physics of nonequilibrium processes. The evolution of all systems in the universe may be coupled to decreasing their free energy in the least possible time [
184]. Thus, the living systems and biological macromolecules can be regarded as manifestations of physical principles about dissipation intensity rather than ends in themselves [
185].
We mainly dealt with the academic interest in answering the why and how of life emergence during universal thermodynamic evolution. However, there is also a practical goal of enhancing the desired activity of natural enzymes or competing with nature in a rational design of artificial enzymes with better catalytic performance. These research fields are still in their infancy. Microwave irradiation can enhance enzyme activity and entropy production under chemiosmotic conditions [
186]. Also, the enhancement in these biochemical and physical parameters can result from distal mutations that do not change individual equilibrium constants for each catalytic step or the overall equilibrium constant of the reaction [
187]. De novo enzyme design for green chemistry and medical goals has a huge potential [
188,
189,
190,
191,
192,
193]. It has been recently explored by combining computational methods and directed evolution experiments [
194]. Still, something needs to be added to our insights about enzymatic catalysis. Artificial enzymes are generally inferior in catalytic efficiency compared to their natural counterparts [
190]. While role or reorganization energy is recognized in rational protein design [
190], that is not the case with the catalytic efficiency to dissipation proportionality for uni-uni enzymes we described in this paper. After all other means are employed to identify possible beneficial mutations for increasing the catalytic efficiency with a given substrate, the computer-aided enzyme design can be extended with an additional selection for higher overall entropy production
. In principle, mutations can be predicted based on their contribution to total entropy production, not only their contribution to the transition state stabilization and reorganization energy.