Theoretical Study on the Kinetics of the Rubisco Carboxylase Reaction by a Model Based on Quantum Chemistry and Absolute Reaction Rate Theory.

The rate of the Rubisco carboxylase reaction is evaluated by statistical mechanics and hybrid density functional theory (DFT). The Rubisco molecular model given by Kannappan et al. was modified and used in the present calculation. The activation energies of CO2 addition reaction, H2O addition reaction, C2-C3 bond scission, and C2 protonation are estimated. We calculated the turnover number (TON) for each of the four reaction steps based on a revised absolute reaction rate theory, which became applicable to soft matter reactions. The molecular parameters used in TON calculations were obtained by DFT calculations. The TON of the total Rubisco reaction was finally evaluated using rate equations. The calculation in a vacuum gave the total TON to be around 5 × 10-5, which was much lower than the experimental value. The DFT calculation in water solvent gave the total TON to be around 0.1, which agreed reasonably well with experimentally reported values (∼2.71). The rate-limiting process was the scission reaction. The present calculation showed that both the phosphate groups in the substrate accelerate each reaction step. The present calculation showed that a more comprehensive molecular model including enolization and quantum chemical methods is necessary to make a more precise reaction model including the irreversibility of some reactions.

Introduction

The global climate change is due to anthropogenic carbon dioxide (CO2) emissions.[1] The Sustainability Science Consortium was established in August 2010 in Japan for the purpose of better cooperation among various fields such as physics, chemistry, engineering, biology, economics, politics, and quality control. Since the reduction of fossil fuel burning is becoming an urgent issue, various types of alternatives, renewable energy resources have been proposed,[2] such as solar energy, biomass, wind power, and wave power. Among other research studies aiming to mitigate global warming, the mechanism of photosynthesis is extensively studied.[3−7] Rubisco plays a key role in photosynthesis. At the stage of year 2008, crystallographic, mutagenesis, kinetic, and computational studies on Rubisco over 3 decades have revealed much about its catalytic mechanism and the role played by several active-site residues. However, key questions remain unanswered. Specific details of the carboxylase and oxygenase mechanisms, required to underpin the rational re-engineering of Rubisco, were still speculative. Kannappan et al. address critical gaps in knowledge with a definitive comprehensive computational investigation of the mechanism of carboxylase activity at the Rubisco active site. Density functional theory calculations (B3LYP/6-31G(d,p)) were performed on active-site fragment models of a size up to 77 atoms, not previously possible computationally, and gave a reasonable energetics of the total enzyme reactions on Rubisco.[6] Then, a more precise and large-scale model of the active center of Rubisco is reported.[7a−7c] The mechanism should be compared with the experimental results. Some experimental research to analyze the reaction rate of each reaction step is in progress.[8] We presented the turnover number (TON) of the total reaction as well as of each reaction step. It will be a future study to assign the reported experimental isotopic study[8] to the theoretically obtained elemental reaction steps.

A theoretical study of the reaction between solid catalytic surfaces and substrate gasses in a vacuum chamber is presented prior to that of the reaction in protein enzymes. In the first stage of solid surface studies, the conventional theory of absolute reaction rates[9] is used. The theory was developed more than 80 years ago. Quantum chemistry, however, was not usually used[10] because of the limited performance of computers. Therefore, the numerical accuracy of the calculations was limited.[10] More recently, one of the current authors (S.O.) reported the systematical simulation of CO2 and H2 adsorption on the surface of Cu and Zn microclusters in a vacuum chamber,[11] based on semiempirical molecular orbital calculation PM6 and the theory of absolute reaction rates, as well as the technique of statistical mechanics, providing the results that agreed reasonably well with experimentally reported values. The theory of absolute reaction rates, proposed in ref (11), gives a more rigorous introduction of the Langmuir adsorption isotherm, than the conventional one,[9] whereas the partition function of the adsorbed substrate in the transition state was approximately obtained as the partition function of the substrate–surface complex divided by the partition function of the isolated surface. This approximation is not reasonable in the case of soft matter because an enzyme is reformed largely with and without a substrate.

As a natural extension of our previous studies,[10−14] we present here a theoretical study. The purpose of the present paper is to report the physical model to give the kinetic rate of the Rubisco carboxylase reaction and the turnover number (TON) for each step of the reactions shown in Figure (addition of CO2, hydration, scission, and protonation). It has not been reported so far even though this is a critical value for the future catalyst design learning from this splendid natural system. We report the TON of the total Rubisco CO2 addition reaction using rate equations.

Figure 1

Energy diagram relative to Structure O (Rubisco + substrate + CO2 (gas)) with the kinetic scheme of the Rubisco reaction in water solvent. NI to NIX are the population in Structures I to IX. TON is the turnover number. VT.CO, VT.H, VT.SUBSTRATE/2, and VT.H are the thermal velocity of CO2, H2O, the half-separated substrate, and hydrogen atom, respectively. Eq. means the thermal equilibrium. δ is the length of the plateau of a transition state.

In general, the chemical reaction rate, k, is given bywhere EA is the activation energy, T is the temperature, KB is the Boltzmann constant, and A is a constant called the frequency factor. The quantum chemical calculation gives the value of EA. The theory of absolute reaction rates,[9] which is based on statistical mechanics, can be employed to derive a mathematical formula for the frequency factor. Each of the reaction steps on Rubisco, except the CO2 capture step, can be calculated by a slightly modified version of Eyring’s original theory of absolute reaction rate,[9] to be shown later.

On the other hand, the rate of the CO2 capture reaction must be calculated by considerably modifying Eyring’s original theory of absolute reaction rate,[9] not only because Rubisco has a surface site on which a single CO2 molecule can adsorb,[11] but also because Rubisco is soft matter. A statistical mechanical calculation[15] gives a mathematical formula. It gives the surface coverage, θ (total number of stably adsorbed CO2 molecules divided by the total number of sites), as a function of the temperature T, and the CO2 gas pressure p. In the present paper, we calculate the quantitative value of θ for the CO2 molecule in the transition state. We developed a more rigorous absolute reaction rate theory than that of the previous study.[11] That is, we calculate the partition function of the transition state as a single unit, without separating the enzyme and the substrate. Thus, the theory is reasonably applied to soft matter calculations. We used ab initio quantum chemical calculations to determine the qualitative values of molecular parameters to be used in the statistical mechanical calculations, such as chemical bonding energies, frequencies of molecular normal mode vibrations, and activation energies. When the value of θ at the CO2 adsorption for the transition state is obtained, the absolute reaction rate can be calculated in a straightforward way, through the theory of absolute reaction rates.[9]

The current model consists of only the nearest neighbor parts acting in the Rubisco carboxylation reaction. We made only a slight modification to Kannappan’s minimal molecular model[6] by adding a second phosphate base because the real substrate of the real Rubisco reaction has two phosphate bases. Kannappan et al. suggested that the dominant steps in the Rubisco carboxylase reaction are irreversible. The present study also showed that reasonable TON values are obtained when assuming that only a forward reaction occurs and no reverse reaction occurs in each step, in accordance with Kannappan’s proposal. The results of the calculation showed that both the phosphate bases accelerate all of the steps of the reaction.

Model and Details of Quantum Chemical Computation

The molecular model of Rubisco and quantum chemical calculations are almost the same as the preceding study given by Kannappan,[6] except for a small addition. The present section shows the details of our molecular model of Rubisco and quantum chemical calculations, clarifying on what basis the kinetic model is constructed.

We introduced a Rubisco and substrate (RS) structure, which is obtained by removing CO2 from Structure I. Structure O consists of CO2 and RS. Kannappan et al. assumed 4 series of reactions and 10 structures (O, I, II, III, IV, V, VI, VII, VIII, and IX) in their model of the total reaction process of Rubisco (see Table ).[6] Following their model, these reactions are the addition of a carbon dioxide molecule to the substrate (reactant O, transition state II, and product III), hydration of the substrate (reactant III, transition state IV, and product V), scission of a C2–C3 bond in the substrate (reactant V, transition state VI, and product VII), and protonation of the substrate (reactant VII, transition state VIII, and product IX). The sequence of reactions is shown in Scheme . Each of the nine structures is a loosely coupled cluster consisting of Mg2+ with fragments of surrounding residuals and the substrate (d-ribulose-1,5-bisphosphate, RuBP).

Table 1

Reactants, Transition States, and Products of Rubisco Reactions

reaction	reactant	transition state	product
O to III	O	II	III
III to O	III	II	O
III to V	III	IV	V
V to III	V	IV	III
V to VII	V	VI	VII
VII to V	VII	VI	V
VII to IX	VII	VIII	IX
IX to VII	IX	VIII	VII

Scheme 1

Sequence of Reactions Catalyzed by Rubisco, in Which Only Substrates that Have Two Phosphate Bases Are Shown

We adopted two phosphate models because we assumed that both electron-withdrawing phosphate groups of RuBP play an equally important role in carboxylation, hydration, and C2–C3 bond-scission processes (see Section S.IV for details). A fragment of residual HIS327 was added to balance the total structure respecting the previously reported theoretical study.[16]

In this paper, quantum chemical calculations were carried out to analyze the reactions. All geometry optimization procedures and vibrational frequency analyses confirming the minima and transition states were performed using the Gaussian 16 program package[17] with hybrid density functional theory (DFT) calculations.[18] The exchange–correlation term was considered in the B3LYP functional.[19−21] The 6-31G(d,p) basis sets[22,23] were adopted for all atoms.

The structure optimization was applied only to the additional parts (i.e., a second phosphate group and residual HIS327) by freezing the original structure[6] in a vacuum. Then, a full structural relaxation was carried out for whole models in the vacuum. Each of the structures optimized in a vacuum was reoptimized in diethyl ether solvent and in water solvent, wherein the solvent effect was incorporated using the solvation model based on the density (SMD) method supplied by Gaussian 16.

The energies for Structures I to IX, having two phosphate bases, are shown in Figure . Optimized Structures for I to IX are shown in Figures and 3.

Figure 2

Optimized Structures I to VI. Each structure is the nearest neighbor part of the reaction site of a Rubisco, including two phosphate bases. In all figures, atoms of Mg and CO2 are highlighted by the ball and stick style, H2O, RuBP, CABP, and GAP are shown in the licorice style, and the residuals are presented in the wireframe style for clarity.

Figure 3

Optimized Structures VII to IX.

Structure I consists of the active site of the Rubisco molecule, one substrate molecule, and loosely coupled CO2. One RS structure consists of the active site of the Rubisco molecule and one substrate molecule. One RS structure and one gaseous CO2 molecule (Structure O) precede Kannappan’s Structures I to IX in the sequence of reactions catalyzed by Rubisco. The details are shown in Section . The relative value of energy for Structures I to IX and the activation energies of these three reactions are summarized in Table .

Table 2

Energies of Rubisco Reactions

structure	relative energy (kcal/mol)	activation energy (kcal/mol)
	(a) In a Vacuum
CO2 + RUB.SUB (O)	0.00
I	–0.03
II	8.94	8.97
III	5.82
IV	14.56	8.74
V	4.98
VI	17.52	12.53
VII	12.40
VIII	16.02	3.62
IX	–15.01
(b) In Diethyl Ether
CO2 + RUB.SUB (O)	0.00
I	1.18
II	9.21	8.02
III	4.35
IV	10.58	6.23
V	2.80
VI	16.65	13.84
VII	10.96
VIII	15.52	4.56
IX	–16.16
(c) In Water
CO2 + RUB.SUB (O)	0.00
I	–4.84
II	0.33	5.18
III	–3.75
IV	5.60	9.35
V	–5.38
VI	10.11	15.50
VII	7.95
VIII	12.01	4.05
IX	–22.65

Computational Model and Theory of the Kinetic Model to Give the Turnover Number

General Theory; Absolute Reaction Rate for Enzyme Reactions

Let us call the population of transition state (TS) NTS, the length of the plateau of the transition state δ, and the thermal velocity of the moving part of TS vT. The reaction rate RP, which is the number of products produced in unit time, is given by[9]The number of enzyme molecules concerning this reaction defined as NENZYME, and the turnover number TON is given by the definitionBy inserting eq to eq , we getHere, we assume[11]The value of and vT can be derived using statistical mechanics. Details are described in the Supporting Information (eqs S.I.4-8, S.I.4-10, S.II.5-3, and S.III.6-4).

Application to Each Step

Adsorption Step; Addition of CO2 to the C2 Atom of the Substrate

Definition of the Reactant, the Precursor State, the Transition State, and the Products of the Reaction

Structure O (see Figure ) consists of the reaction center of Rubisco, the substrate (enediolate), and free gaseous CO2. Structure I (see Scheme and Figure ) is considered the precursor state of CO2 addition, which consists of the reaction center of Rubisco, the substrate (enediolate), and loosely captured CO2. The chemical reaction from Structure I to Structure III, which is the CO2 addition reaction, occurs through a transition state (i.e., Structure II in Figure ). Structure III (see Scheme and Figure ) consists of the reaction center of Rubisco and the substrate firmly coupled with CO2 (CKABP). We call the number of this product III as NIII.

No transition state between Structure O (see Figure ) and precursor Structure I was found in the present calculations. The transition Structure II consists of the reaction center of Rubisco, the substrate (enediolate), and loosely captured CO2 drifting within a plateau of length δ with a thermal velocity vT.CO. We call the number of the transition states NII. The structure consisting of the reaction center of Rubisco and the substrate (enediolate), is called the RS complex, hereafter, whose number is called NRS. Structure O consists of the RS and the gaseous CO2 molecule. We assume that there are NCO number of CO2 gaseous molecules in volume V.

Canonical Ensemble for the Reaction Rate of CO2 Addition to the C2 Atom of the Substrate

First, we describe the reaction from Structure O to Structure III. In absolute reaction rate theory, the reactants (Structure O = RS + gaseous CO2 molecule) and the transition states II are assumed to be in thermal equilibrium. There are NRS number of RS and NCO number of gaseous molecules, before the addition reaction occurs. During the addition reaction, NII number of transition states exists in equilibrium. Then, in the equilibrium, (NRS – NII) number of RS remain uncombined with CO2, and (NCO – NII) CO2 gaseous molecules remain uncombined with an RS. We assume that the volume of gaseous CO2 remains V in the equilibrium. The free energy (ΔH – T·ΔS) of the total system can be given as a function of NII. The value of NII fluctuates by thermal noise. The RS, the gaseous CO2, and the transition states II are assumed to be in thermal equilibrium. Thus, they form the canonical ensemble with the average value, ⟨NII⟩. According to eq , the turnover number for CO2 addition to the C2 atom of the substrate is given in eq , and the detailed derivation is in eq S.I.4-12 of the Supporting Information.where PCO is the partial pressure of the CO2 gas, δ is given by eq , KB is the Boltzmann constant, T is the temperature, mCO is the mass of a CO2 molecule, h is the Plank constant, and NVM.II is the total number of the vibrational normal modes of a transition state II. One normal mode out of NVM.II modes has an imaginary number frequency, reflecting the transition state. ν is the real number vibration frequency of the iII-th vibrational normal mode of the transition state II. NVM.RS. is the total number of the vibrational normal modes of each RS molecule. All of the vibrational normal modes of an RS molecule have a real number of vibrational frequencies. HACT.ADD is the activation energy of the CO2 addition, which is the energy level of II measured from the energy level of a gaseous CO2 molecule and an RS (i.e, Structure O). The values of ν, ν, and HACT.ADD are calculated by quantum chemical calculations.

j(T) is the rotational and vibrational partition function of a CO2 molecule,[15] as defined belowwhere ICO is the moment of inertia of a CO2 molecule, NVM.CO is the total number of the vibrational normal modes of a CO2 molecule, and ν is the frequency of the i-th vibrational normal mode of a CO2 molecule.

Second, we describe the thermal equilibrium between Structure O (RS and CO2 molecule) and Structure I (the precursor state). Here, we include Structure I. Structure O and the precursor state (Structure I) are assumed to be in thermal equilibrium. In the present quantum chemical calculations, no transition state for the reaction from Structure O to Structure I was found. We assume that there were originally NRS number of RS and NCO number of gaseous molecules, before the capture of gaseous CO2 molecules, that is, before the formation of Structure I. After the equilibrium state is achieved, there are NI precursor states in Structure I. Thus, in the equilibrium, (NRS – NI) number of RS remains uncombined with precursor Structure I and (NCO – NI) number of gaseous CO2 molecules remain uncombined with RS. In this stage, it is assumed that the value of N has no influence on the population balance. We assume that the volume of gaseous CO2 remains V in the equilibrium. The free energy of the total system is given as a function of NI. The value of NI fluctuates by thermal noise. The RS, the gaseous CO2 molecules, and precursor Structure I form the canonical ensemble with the average value, ⟨NI⟩. The complete derivation for ⟨NI⟩ is given in eq S.III.6-4where NVM.I is the total number of the vibrational normal modes of precursor Structure I, and ν is the vibrational frequency of the i-th normal mode in Structure I. Note that normal vibrational frequencies of precursor Structure I are all real numbers. HCAPI is the heat generated when one CO2 gaseous molecule is captured by RS and one precursor state (Structure I) is formed. The j(T) is given by eq .

The coverage θ is derived in eq S.III.6-6 in the Supporting Information, and the result is

According to absolute reaction rate theory, the reactants (precursor Structure I) and the transition states (Structure II) are assumed to be in thermal equilibrium. There are NI number of precursor structures, before the I to III reaction occurs. In the equilibrium state after the I to III reaction starts, there are N transition states. Then, in the equilibrium, (NI – NII) number of precursor Structure I remains without changing to Structure III. The free energy of the total system is given as a function of N. The value of N fluctuates by thermal noise. Precursor Structure I and the transition states II are assumed to be in thermal equilibrium. Thus, they form the canonical ensemble with the average value, ⟨NII⟩. According to eq , the turnover number of this reaction from Structures I to III is given in eq S.II.5-11 in the Supporting Informationwhere NVM.TS is the number of vibrational normal modes of the transition state (Structure II). One vibrational frequency at the transition state is an imaginary number. ν is another real number vibrational frequency belonging to the transition state. HACT is the activation energy of the reaction, which is the energy level of II relative to I. NVM.R is the total number of the normal vibration modes of the reactant I. ν is the vibrational frequency (all are real) of the reactant (Structure I).

In the case of the present study, the equationcan hold, when the following two conditions are satisfied

Equation states that the reaction rate of the direct gaseous CO2 molecule addition to form a product III is nearly equal to the reaction rate of transition from I to III, wherein state I is in equilibrium with gaseous CO2 molecules and the CO2 coverage on I is θI (see Figure ).

Hydration at the C3 Center

Following the CO2 addition reaction, the next reaction from Structure III to Structure V (see Scheme and Figure ) occurs, having the transition state Structure IV. This is a reaction of hydration at the C3 carbon of the substrate. In absolute reaction rate theory, the reactants (Structure III) and the transition states (Structure IV) are assumed to be in thermal equilibrium. There are NIII number of reactant III structures, before the reaction, and then this reaction starts. There are NIV number of transition states in equilibrium. In this equilibrium, (NIII – NIV) number of reactant III structures remain without changing into Structure V. The free energy of the total system is given as a function of NIV. The value of NIV fluctuates by thermal noise. The reactant III and the transition states IV are assumed to be in thermal equilibrium. Thus, they form the canonical ensemble with the average value, ⟨NIV⟩. According to eq , the turnover number of this reaction from Structures III to V is given bywhere vT.H is the thermal velocity of a H2O molecule. In eq S.II.5-11 in the Supporting Information, it is shown thatwhere mH is the mass of a H2O molecule, NVM.TS is the total number of the vibrational normal mode of the transition state IV, NVM.R is the total number of the vibrational normal mode of the reactant III, and HACT is the activation energy of the reaction, which is the energy level of IV relative to III (see around TON(III→V) in Figure ).

C2–C3 Bond Scission

After the hydration of the C3 carbon, the next reaction continues from Structure V to Structure VII (see Scheme and Figures and 3). The transition state of this reaction is with Structure VI (see Scheme and Figure ). This is a reaction of the C2–C3 bond scission. Structure VII is a loosely combined aggregate including two fragment products from the divided substrate molecule. In an absolute reaction rate theory, the reactants (Structure V) and the transition states (Structure VI) are assumed to be in thermal equilibrium. There are NV number of reactant V before this reaction occurs. In the equilibrium state after this reaction starts, there are NVI number of transition states. Then, in the equilibrium, (NV – NVI) number of reactant V structures remain without changing into Structure VII. The free energy of the total system is given as a function of NVI. The value of NVI fluctuates by thermal noise. The reactant V and the transition state VI are assumed to be in thermal equilibrium. Thus, they form the canonical ensemble with the average value, ⟨NVI⟩. According to eq , the turnover number of this reaction from Structures V to VII is given bywhere vT.SUBSTRATE/2 is the thermal velocity of one of the fragments of the divided substrate molecule. In eq S.II.5-11 in the Supporting Information, it is shown thatwhere mSUBSTRATE/2 is the mass of one of the fragments of the divided substrate molecule, NVM.TS is the total number of the vibrational normal mode of the transition state VI, NVM.R is the total number of the vibrational normal mode of the reactant V, and HACT is the activation energy, which is the energy level of VI relative to V (see around TON(V→VII) in Figure ).

C2 Protonation

After the C2–C3 bond-scission reaction, the final reaction from Structure VII to Structure IX (see Scheme and Figure ) occurs. The transition state of this reaction is realized in Structure VIII (Figure ). This is C2 protonation. In absolute reaction rate theory, the reactants (Structure VII) and the transition states (Structure VIII) are assumed to be in thermal equilibrium. We assume that there are NVII number of reactants VII, before this reaction occurs. In the equilibrium state after this reaction starts, there are NVIII number of transition states. Then, in the equilibrium, (NVII – NVIII) number of reactant VII remains without changing to Structure IX. The free energy of the total system is given as a function of NVIII. The value of NVIII fluctuates by thermal noise. The reactant VII and the transition state VIII are assumed to be in thermal equilibrium. Thus, they form the canonical ensemble with the average value, ⟨NVIII⟩. According to eq , the turnover number of this reaction from Structures VII to IX is given by

where vT.PROTON is the thermal velocity of one proton. In eq S.II.5-11 in the Supporting Information, it is shownwhere mPROTON is the mass of one proton, NVM.TS is the total number of the vibrational normal mode of the transition state VIII, NVM.R is the total number of the vibrational normal mode of the reactant VII, and HACT is the activation energy, which is the energy level of VIII relative to VII (see around TON(VIII→IX) in Figure ).

Rate Equations

Now, we have prepared every necessary piece, and we are ready to discuss the feature of reactions shown above using rate equations for the θ values. Here, θ1 was given by eq . θ3, θ5, and θ7 are the values of θ (each population divided by NRUB.SUB) for Structures III, V, and VII, respectively.

The stationary state conditions are

By inserting eq into eq , we obtainwhereThe TON(I→III), TON(III→V), TON(V→VII), and TON(VII→IX) are TONs for reactions I to III, III to V, V to VII, and VII to IX, respectively. The total TON(TOTAL) of the stationary state is given by

Results and Discussion

Results Given by Rate Equations

We calculated the value of TON(TOTAL) considering only the forward reactions. The rationalization will be discussed later. That is, we assumed

By utilizing these rate equations that we derived in this study, we evaluated the kinetics of the reactions for three representative solvents; in a vacuum, in diethyl ether, and in water with CO2 pressures 4 × 10–4 and 8 × 10–4 atm as a function of temperature; the values of TONs thus given by eq are shown in Figures –6. For example, the TON(TOTAL) at 300 K in 4 × 10–4 atm CO2 gas in the water solvent is 0.113 reflecting the double phosphate Rubisco model in the present work. On the contrary, the TON of reactions III to V of the single phosphate Rubisco reaction at 300 K in a vacuum is 6.38 × 10–5. This shows that the double phosphate model in water solvent seems to be more reasonable than the single phosphate model in a vacuum. Notice that the experimentally reported[2] value of the TON for the total Rubisco reaction is around 1. However, it is a matter of course that further studies are necessary, such as more extended molecular dynamics calculations to examine whether water molecules are distributed in the active center of Rubisco as the solvent.

Figure 4

Results for the absolute reaction calculation in a vacuum (double phosphate model).

Figure 6

Results for the absolute reaction calculation in water solvent (double phosphate model).

We calculated the TON of the total CO2 addition reaction and the TON of each reaction step in the total reaction. In Figure , the calculated values of the TON in a vacuum at a CO2 partial pressure of 4.0 × 10–4 atm, which is the partial pressure of CO2 in the present atmosphere are shown. The calculation of the TON at a CO2 partial pressure of 8.0 × 10–4 atm is also shown to illustrate the dependence of the TON on the CO2 partial pressure. θ1, θ3, θ5, and θ7, are the population of Structures I, III, V, and VII divided by the population of total Rubisco, respectively. TONs calculated in diethyl ether solvent and in water solvent are shown in Figures and 6.

Figure 5

Results for the absolute reaction calculation in diethyl ether solvent (double phosphate model).

The calculation in water solvent gives the reasonable value of the TON, which is consistent with the highly polar environment of the Rubisco reaction center. The rate-limiting process is the reaction from Structure V to Structure VII, which is the scission reaction.

About Irreversible Reactions

Kannappan et al. suggest that hydration at the C3 carbon (i.e., reactions III to V) is an irreversible reaction and only a forward reaction occurs, followed by no reverse reaction.[1] This suggestion may be natural since, in the product, the distance between a H and an O atom (originally in water molecule) is as far as 2.394 Å. The C2–C3 bond scission (i.e., reactions V to VII) is also an irreversible reaction.[1] It should be noted that in the transition state, the distance C2–C3 is 2.432 Å and in the product, the distance C2–C3 is as far as 4.142 Å. The C2 protonation (i.e., reactions VII to IX) is also irreversible because Structure IX is energetically extremely stable.

Based on absolute reaction rate theory, we also show that a reasonable total TON value is obtained only when assuming that addition of CO2 (i.e., reactions I to III) is irreversible. To the best of our best knowledge, genuine absolute reaction rate theory neglects the detailed collision geometries and cannot explain the irreversibility of reactions I to III, III to V, and V to VII. Detailed studies will be necessary to give a complete reaction rate theory for irreversibility in the future; here, we introduced a preliminary step of the qualitative model as follows.

The activation energy of the reverse reaction of the C2–C3 bond scission (i.e., reactions VII to V) is 2.16 kcal/mol. Therefore, the force that the particle feels, F, isBy the diffusion constant for the particle in the protein D, the drift velocity vD of the particle driven by the force F is given by Einstein’s relation as

In the condition that the drift velocity by the force F and the diffusion balances, the following equation can holdHere, n is the density of the particle, and it can be naturally considered asThen, we obtainTherefore, we obtainWhen we assume T = 300 K, this value isTherefore, taking into account the distances shown above, we put this in the x.Now, we can evaluate the density of the particle that we must considerAlthough the discussion given here may be preliminary, we can conclude that the density of the particle near the transition state is e1094 times smaller than that of the particle near the product. Therefore, the possibility of the reverse reaction to occur may be almost negligible.

Thus, it is argued that hydration and scission are irreversible because of the large distance between the atoms in the resulting molecular complexes. However, models (refs (7)), which include more of the enzyme’s detailed structure, suggest much shorter distances. Therefore, there remains a large question of whether the conclusions drawn here about irreversibility cannot be extrapolated to the enzyme reaction because the model ignores the enolization step and the molecular model of the enzyme is oversimplified. A more complete molecular model including enolization is necessary to discuss the irreversibility.

Origin of Errors in the Calculated Value of the Turnover Number

The experimentally obtained value of the turnover number of the Rubisco reaction of C3 plants is reported[24] to be around 2.71 at 25 °C. The calculated value of the turnover number given by the kinetic model of the present study is 1.13 × 10–1 at 27 °C. Therefore, the kinetic scheme calculation of the present paper estimates the total reaction rate to be 24 times smaller than that of the experimentally obtained value. There are several possible origins of this error as shown below.

The kinetic scheme of the present paper is based on the conventional transition state theory. In conventional transition state theory, several assumptions are involved. The first assumption is that the transition state is a point of no return and if the molecule once passed over the transition state to the product, it never returns to the reactant. If this assumption fails, the density of the transition state cannot be calculated assuming simple thermal equilibrium between the reactant and the transition state. In such a case, the calculated reaction rate is usually adjusted using the transmission coefficient κ. The typical value[25] of κ is found to be around 1/2. This means that the theory overestimates the reaction rate by factor 2, which is not the dominant origin of the error of the present kinetic calculation.

Ballistic scattering by the potential energy surface near the transition state often occurs at high temperatures, and tunneling often occurs in a small system at low temperatures. Rubisco is a large molecule and the reaction occurs near room temperature. So correction for ballistic scattering or the tunneling effect is not expected to be significant in the present kinetic calculation. Our calculations show that the transition state has only one imaginary numbered vibrational frequency, and all of the activation energies are positive numbers and the reactions are irreversible. Still, it is not clear at present whether the reaction coordinate is one dimension and completely separated from other freedom of movement because the system is soft and it is expected that some molecular vibrations deviate from the harmonic oscillator largely. A detailed quantum molecular dynamics study will be necessary to clarify this problem.

In the kinetic scheme of the present paper, the original transition state theory is modified so as to calculate soft matter, and the length of the plateau of the transition state δ is used explicitly in the equation giving the turnover numbers. The value of δ = 1.0 × 10–10 m is used in the present calculation, but the change in the absolute value of the turnover number was less than 0.001% even when the value of δ = 1.0 × 10–10 m or δ = 1.0 × 10–9 m is used in the calculation. Therefore, the error in the value of δ is not dominant in the error of the turnover value calculated in the kinetic scheme of the present paper.

In the kinetic scheme of the present paper, the rate-limiting process (except for the enolization step) was assigned to be the scission reaction, whose activation energy was calculated to be 15.49 kcal/mol. The experimentally reported value of the activation energy of the total Rubisco reaction of the C3 plants’ Rubisco reaction is[24] 12.9 ± 0.4 kcal/mol. By calculating at 300 K, we will find that the present kinetic calculation underestimates the turnover number by about 70 times. By calculating 70/24, the present calculation estimates the turnover number as only 2.9 times larger than the experimental values, if the activation energy of the scission reaction is corrected according to the experimentally reported value. At present, it is not clear that 2.9 times overestimation can be corrected somewhat by factor κ.

The activation energy of the reaction calculated by the present quantum mechanical calculation is 15.49 kcal/mol and is larger than the reported experimental value[24] 12.9 kcal/mol by 2.39 kcal/mol. The present calculation shows the stabilization energy of 3PGA to be −22.65 kcal/mol. However, it should be noted that the measured enthalpy of the overall carboxylase reaction is[26] found to be around −5 kcal/mol. This shows that the enolization step, which was not taken into account in the present kinetic model, should be included in the calculation and/or a more complete molecular model.[27]

It should be noted that the rate of reactions at room temperature is very sensitive to the value of the activation energy when compared with those at high-temperature reactions. The current model consists of only the nearest neighbor parts acting in the Rubisco carboxylation reaction. This oversimplification is possibly one of the significant reasons for the energy overestimation of 2.39 kcal/mol.

We adopted the B3LYP method in the quantum chemical calculation, and still, it is known that the B3LYP method overestimates the hydrogen bonds, and it is a question of how well the chosen basis set deals with Mg2+.

In the present state, it is clear that the Rubisco reaction at room temperature is an extremely delicate chemical reaction, and a sufficient number of active-site residues must be taken into account, and a more complete (e.g., QM/MM) description of the enzyme environment is essentially necessary to make a reasonable molecular model, and a more precise quantum mechanical method than B3LYP is necessary. It is expected that the kinetic scheme of the present paper will be a good tool with which to explore alternative molecular mechanisms for this reaction including issues regarding protonation states and the origin of the hydrating water molecule.

Conclusions

The absolute rate of the Rubisco carboxylase reaction is calculated using ab initio quantum chemical calculations and absolute reaction rate theory. We adopted a slightly modified Kannappan’s molecular model[6] of Rubisco. The present paper reports the exact absolute reaction rate model owing to the previously studied models.[12,13] Special attention is paid on how to make a kinetic model for soft matter enzyme reactions precisely. Kannappan et al. suggested that carboxylase reactions are irreversible.[6] The present paper also supports this. Our calculations show that assuming a reaction in water solvent is more reasonable than in a vacuum. Present calculations showed that both electron-withdrawing phosphate groups of RuBP accelerate carboxylation, hydration, and C2–C3 bond-scission processes. We proposed a simple model that explains why the second phosphate group, which is far from the C2–C3 bond, also accelerates the scission reaction (see the Supporting Information). In the present paper, a simple model was proposed to qualitatively explain the irreversibility. The calculated value of the turnover number is about 24 times smaller than the experimentally obtained value. The main cause of this error is estimated to be the overestimation of the activation energy by 2.59 kcal/mol. It is clearly shown that the Rubisco reaction is very delicate and a sufficient number of active-site residues must be taken into account. A more complete (e.g., QM/MM) description of the enzyme environment is essentially necessary and the enolization step must be taken into account to give a qualitative value of the turnover number and to discuss the irreversibility of some reactions precisely.