Proton tunnelling in hydrogen bonds and its implications in an induced-fit model of enzyme catalysis

(i) Motion of the proton between two atoms/molecules.

The locations between which the proton tunnels back and forth can be regarded as the sites of an interaction network in such a way that each site j = {1, 2} is associated with the bonding orbital σ_Xj − H (figure 3a). Then the proton is expected to move from one site to another in accordance with the Hamiltonian

where ${\hat{n}}_j=a_j^{†}{a}_j$ is the proton number operator at site j, a^†_j and a_j are, respectively, proton creation and annihilation operators that obey the fermion anticommutation relations. On-site energy W_j can be taken to be the total potential felt by a proton at the jth site. J₁₂ stands for the orbital interactions that drive proton tunnelling. V₁₂ is introduced to penalize the case in which there is one proton at each site. λ₁₂ is a constant responsible for the total intermolecular interactions between X₁ and X₂.

• Download figure
• Open in new tab
• Download PowerPoint

Each of the coefficients in equation (2.1) exhibits a different functional dependence on geometric parameters such as the length of the single X_j − H covalent bond (r_j), the separation between X₁ and X₂ (R₁₂) and the angle (ϕ₁₂≡∠HX₁X₂). The key coefficient here is the hopping constant J₁₂. Its functional dependence on the geometric parameters can be written in a similar way to the coupling constant of diabatic state models [43,44] as:

where J₀, b₀ and R₀ depend on the chemical identities of X₁, X₂ and their environment. This functional dependence guarantees the directionality that is a critical property of H-bonds. J₁₂ becomes J₀ e^{−b₀(R₁₂−R₀)} for linear H-bonds and decreases slightly by the increase of ϕ₁₂. When ϕ₁₂ reaches to the value of π/2, it vanishes and the proton can not tunnel between X₁ and X₂ anymore even if the separation R₁₂ is sufficient for the formation of a quantum H-bond. However, it is important to note that our results in this study do not require the exact form given in (2.2). The only requirement is that J₁₂ should be non-zero at small bond angles, and should vanish when the bond angle becomes larger.

To obtain a pseudo-spin Hamiltonian by preserving the anti-commutation relations, we apply the Jordan–Wigner transformation to a_j, a^†_j and ${\hat{n}}_j$ in (2.1) with the convention for Pauli z operator that σ^(j)_z = |0_j〉〈0_j| − |1_j〉〈1_j|. Then, we arrive at a two-spin Heisenberg XXZ model under one homogeneous and one inhomogeneous magnetic field,

with parameters J_x =  − J₁₂/2, J_z = V₁₂/4, $\tilde{\lambda }={\lambda }_{12}+(2W_1+2W_2+V_{12})/4$, B =  − (W₁ + W₂ + V₁₂)/4, b =  − (W₁ − W₂)/4. The eigensystem of this pseudo-spin Hamiltonian can be written in an increasing order of eigenvalues as: where ${\eta }_{\pm }=(b\pm \sqrt{b^2+J_x^2})/J_x$. The ground state |e₁〉 is just a two-qubit representation of the state α|X₁ − H · · · X₂〉 + β|X₁ · sH − X₂〉 that describes a generic H-bond. Unlike the one-qubit representation of the same state that is used in the spin-boson and the diabatic states models found in the literature, this two-qubit representation can investigate quantum correlations such as entanglement and discord generated in H-bonds. This is an advantage of our quantum informational perspective that stands to yield new insights since there is a high amount of correlations in biomolecular interactions and quantum physics is better suited to describe all correlations contained in a complex system in a rigorous way [45]. Also, the ground state energy e₁ goes to the expected energy of a classical bond, i.e. W₁ + λ, when J_x vanishes. Moreover, this energy seems to be suitable for energy decomposition analysis as it equals W₁ + λ + J₁₂ for symmetric H-bonds.

Besides these properties of the ground state, the last two excited states respectively represent the state of X₁ and X₂ atoms of the ionic intermediate structures S^‡₊ and S^‡₋ formed in water-mediated tautomerization.

In the following, each of the locations allowed for the protons is regarded as a pseudo-spin, explained above. The formation of a H-bond between two atoms/molecules X_j and X_j^′ in the biological scenario under consideration is then described as the coupling of the corresponding pseudo-spins by the Hamiltonian H^(jj′)_HB given in (2.3).

(ii) Tautomeric transitions due to the motion of protons.

Here, we will represent S, P, and all the possible transition states between them by different eigenstates of a single Hamiltonian. This may seem like an oversimplification, since each of these different molecules is actually corresponding to the ground state of a separate Hamiltonian with many electronic levels. However, it is just a projection of all the extrema of the potential energy surfaces of the nucleotide onto a single energy spectrum, representing a unification of all the tautomeric forms along possible reaction coordinates from S to P into a hypothetical generic molecule. This type of modelling, which assigns a global time-independent Hamiltonian for a transition from an initial to a final state is quite common in quantum thermodynamics.

In this respect, we collect all the other degrees of freedom apart from the two pseudo-spins in the molecular structure into a single macromolecular configuration C and label its state as |ζ〉_C. That is to say, the whole state |N〉 representing the nucleotide's tautomeric form N = {S, P, S^‡, S^‡₋, S^‡₊, …} is taken to be the product of its configuration state |ζ_N〉_C and its pseudo-spin state |ψ_N〉₁₂.

Let us first consider the states of S and P. To prevent a direct transition between these forms, both ϕ₁₂ in S and ϕ₂₁ in P are taken to be π/2. Then the hopping constant J₁₂ substituted in H⁽¹²⁾_HB that describes the pseudo-spin interaction in either form vanishes, i.e. |ψ_S〉₁₂ = |10〉₁₂ and |ψ_P〉₁₂ = |01〉₁₂. In other words, the H-bonding interaction between X₁ and X₂ does not have a quantum character. When one of the forms is converted to the other, the energy of this classical interaction also changes. However, no significant difference in the molecular structure is expected due to this change. Thus, we can assign the same configuration state |G〉_C to both forms: |S〉 = |G〉_C⊗|10〉₁₂ and |P〉 = |G〉_C⊗|01〉₁₂.

Similarly, neither the protonation of X₂ nor the deprotonation of X₁ is likely to alter the rest of the molecular structure. Hence, the ionic transition states S^‡₋ and S^‡₊ can, respectively, be described by |G〉_C⊗|00〉₁₂ and |G〉_C⊗|11〉₁₂. Conversely, a more energetic configuration, |E〉_C, should be assigned to the neutral transition state S^‡ as the proton tunnelling between X₁ and X₂ requires a bond angle less than π/2 as is given in figure 1. Then we end up with |S^‡〉 = |E〉_C⊗|e₁〉₁₂.

Note that the ground and excited states of the molecular configuration, |G〉_C and |E〉_C, which we can visualize them like macromolecular logic qubits |0〉_L and |1〉_L, need not include any quantum degree of freedom.

At this point, it is possible to describe tautomerization processes as transitions between energy levels of a single Hamiltonian constructed as:

with an eigensystem: where ϵ −  = e₁(r^′, R^′, ϕ^′ < π/2) and ϵ +  = e₂(r^′, R^′, ϕ^′ < π/2). Also, ${\tilde{\text{S}}}^{‡}$, ${\tilde{\text{S}}}_{-}^{‡}$ and ${\tilde{\text{S}}}_{+}^{‡}$ are the first possible excitations from S^‡, S^‡₋ and S^‡₊, respectively.

(iii) Open quantum system dynamics.

The configuration of atoms/groups and the pseudo-spins are expected to be coupled to different phonon environments in different ways. For example, the configuration states are likely to exchange energy with a heat bath B in a reversible manner

where b^†_k and b_k are phonon creation and annihilation operators associated with the kth oscillator mode of the bath. These bath operators can be related to the rotational vibration modes such as bending or libration modes that change the orientation of X₁ and X₂. Such rotational vibrations may occur due to the collisions with solvent molecules or as a result of the intramolecular nucleotide dynamics. As X₁ and X₂ are covalently bonded to the rest of the molecule, these vibrations should be dependent on the orientation of the whole other atoms/groups close to them.

However, dominant environmental effect on the proton should be originated from the X_j − H stretch vibrations that are not expected to affect rest of the molecule in a significant way. Charge fluctuations in the surrounding molecules may drive these oscillations and we can incorporate them into our model by coupling the position of the proton linearly to the equilibrium positions of phonons through

where ${\tilde{b}}_{k,j}^{†}$ and ${\tilde{b}}_{k,j}$ are phonon creation and annihilation operators of a second heat bath $\tilde{B}$ and they are associated with the kth oscillator mode at the jth proton location. $H_{\tilde{B}S}$ guarantees that the interaction of the pseudo-spins with this second heat bath $\tilde{B}$ destroys the quantum correlations between them and has no further effect on the configuration.

As the magnitudes of E_g and E_e are significantly higher than the other components in the eigenenergies given in equation (2.6), the evolution of configuration and pseudo-spins in the nucleotide can be assumed to be separated in such a way that the former is not affected by the latter, but not vice versa. In this respect, we can solve a master equation for the configuration and use this solution in the evolution of the pseudo-spins. Here, the evolution of either the pseudo-spin or the configuration state is described using the Markovian master equation in the Lindblad form [46], which is one of the key elements of the theory of quantum thermodynamics [47],

where H_LS is the Lamb shift Hamiltonian providing an environment-induced unitary contribution to the dynamics: whereas $\mathcal{D}$ is the dissipator responsible for the irreversible dynamics:

While solving this master equation for the pseudo-spin state, ω and A_j(ω) are, respectively, taken to be the Bohr frequencies and the Schrödinger picture eigenoperators of the pseudo-spin Hamiltonian H_HB in (2.3), i.e. ω = (e_j^′ − e_j^′′)/ℏ and $A_j(\omega )=\sum _{j^{\mathrm{\prime }}{j}^{\mathrm{\prime }\mathrm{\prime }}}|e_{j^{\mathrm{\prime }\mathrm{\prime }}}⟩⟨e_{j^{\mathrm{\prime }\mathrm{\prime }}}|A_j|e_{j^{\mathrm{\prime }}}⟩⟨e_{j^{\mathrm{\prime }}}|$ where A_j are the Pauli z operators in accordance with the equation (2.8). Also, the coefficients S_jj^′(ω) and γ_jj^′(ω) are defined as the imaginary part and one half of the real part of the one-sided Fourier transformation of the thermal correlation function of heat bath $\tilde{B}$.

Parameters of the pseudo-spin Hamiltonian H_HB used in this solution should change depending on the configuration state in line with the equation (2.5). To do so, if the configuration is in the state |G〉_C, ω and A_j(ω) are evaluated using the eigensystem of H_HB that is calculated with a vanishing hopping constant J_x. This makes the dynamics of the pseudo-spins fully separable from each other and corresponds to a pure two-qubit dephasing process. Conversely, a non-vanishing J_x is taken into account when the configuration is in the state |E〉_C. This couples the dynamics of the pseudo-spins to each other and leads to an evolution that brings the pseudo-spin states to a detailed balance only between |e₁〉₁₂ and |e₂〉₁₂ in the stationary state.

On the other hand, the evolution of configuration state should be independent from the pseudo-spin state. Hence, the master equation (2.9) is solved for the configuration state using the Bohr frequencies and eigenoperators of an effective Hamiltonian H_C ≈ E_g|G〉_C〈G| + E_e|E〉_C〈E|. However, the eigenoperators A_j(ω) cannot be directly constructed from H_C and H_BC. The interaction Hamiltonian (2.7) is first decomposed into the form $\sum _j{A}_j^{(C)}\otimes B_j^{(B)}$ with A^(C)_j and B^(B)_j are Hermitian operators on the Hilbert spaces of the configuration C and the bath B. Moreover, the coefficients S_jj^′(ω) and γ_jj^′(ω) are defined for the heat bath B this time. This solution describes a thermalization process.

(iv) Enzyme's conformational changes during the induced-fit mechanism.

Two allowed proton locations in the binding site of E are also regarded as pseudo-spins, e.g. |ψ_E(t₀)〉₃₄ = |01〉₃₄. Furthermore, we take into account the conformation of this site in a similar way and denote it by C^′: |ζ_E(t₀)〉_C^′ = |G^′〉_C^′. However, unlike the nucleotide's conformation C, C^′ is assumed to change in a classical and gradual way because of the intermolecular interaction so that its pseudo-spins are continuously tilted. Thus, the bond angles ϕ₁₃ and ϕ₄₂ become smaller and the resultant decrease in the energy of intermolecular H-bonds strengthens the binding interaction.

In accordance with the biological scenario under consideration (see §2a), the time-scale τ₁ of the binding interaction due to the induced-fit mechanism is supposed to be much larger than τ₂, the time-scale over which E undergoes a conformational change that converts C^′ back to the initial state |G^′〉_C^′.

As the interaction Hamiltonian changes slowly in time relative to the conformational change that follows, the instantaneous eigenstates of the binding interaction are assumed to evolve independently. Hence, the enzyme-substrate complex (ES) is considered to be in the corresponding ground state at every instant of time until it reaches the maximum complementarity between E and S. On the other hand, as E rapidly returns back to its initial state after this point due to the relation τ₂≪τ₁, the reverse conformational change is supposed to be driven by a sudden post-selection measurement on the pseudo-spins of E which accompanies the loss of energy from its configuration C^′ to the heat bath B.

(i) Quantum correlated H-bonds.

We assume without any loss of generality that the quantum correlations of each three H-bonds in figure 2b is sustained until E undergoes a conformational change that converts it back to the initial state. In other words, τ₂ is assumed to be small compared to the decoherence time τ_D enforced by the heat bath $\tilde{B}$. This is quite reasonable as the H-bonded atoms/groups are partially isolated from their environment until the detachment of the ES complex.

As explained in §§2a,b(iv), the second conformation motion of interest drives C^′ from |E^′〉 to |G^′〉 and is governed by the interaction Hamiltonian H_BC^′. During the course of this relatively fast motion, the pseudo-spins of E should detach from the pseudo-spins of S and return back to the initial product state |01〉₃₄. We describe this conformation-induced process by a sudden post-selection measurement M = I₁₂⊗|01〉₃₄〈01| and find that some of the four-qubit entanglement of the reaction intermediate ES can be transferred to the pseudo-spins of S when M converts the pseudo-spins of E back to their initial state:

where $b^{\mathrm{\prime }}=b/\sqrt{b^2+d^2}$ and $d^{\mathrm{\prime }}=\sqrt{1-{b^{\mathrm{\prime }}}^2}$. We can analyse the open system dynamics of the nucleotide using the master equation approach introduced in §2b(iii), starting at t = t₄. The stationary solution becomes

This means that the enzyme converts the nucleotide into the product state with a probability of p(P) = b^′2P_g at equilibrium. Also, neither p(S^‡) nor $p({\tilde{S}}^{‡})$ changes when compared to equilibrium probabilities in the uncatalysed case.

b^′ depends on the strength of the H-bonds formed in the enzyme-substrate interaction, as does the efficiency of the catalysis. If the three H-bonds shown in figure 2b generate a maximal W-type four-qubit entanglement, i.e. $a=b=c=d=e=f=1/\sqrt{6}$, the nucleotide reaches an equilibrium such that p(S) = p(P) = P_g/2 that equals 0.5 for our chosen model parameters.

We can compare this result with another obtained in water-mediated tautomerization. When the amounts of entanglement generated through each water-nucleotide H-bond are chosen to be maximal, the nucleotide reaches an equilibrium with p(S) = p(P) = 0.25. Here, when the overall four-qubit entanglement is chosen to be maximal, the pairwise entanglement generated through each enzyme-nucleotide H-bond is not maximal and the equilibrium probabilities corresponding to both tautomers are twice as high as 0.25. Moreover, the ionic transition states S^‡ ± do not emerge unlike the water-mediated case in which p(S^‡ ± ) = 0.25. In the meantime, the equilibrium probability of the neutral transition state is still kept at the order of 10⁻³⁰ at physiological temperatures.

In this context, not only the conformational change induced in the enzyme increases the likelihood of the formation of quantum H-bonds with the nucleotide, but also the subsequent conformational change that brings the enzyme back to its initial state makes quantum H-bonds a more efficient resource for the tautomerization process.

(ii) Classically correlated H-bonds.

The assumption of τ₂ < τ_D is not necessary to derive our insight into the possible role of proton tunnelling in an induced-fit model of enzyme catalysis. Once the quantum correlations are created by proton tunnelling events through H-bonds, the enzyme can use them as a resource to convert S to P by the same mechanism under consideration even if these correlations rapidly decay to classical correlations during decoherence.

In the case of τ₂ > τ_D, the final state of the pseudo-spins of ES at the end of the binding stage should be a classically correlated state rather than a four-qubit entangled state:

After the pseudo-spins of E are detached from the pseudo-spins of S and returned back to their initial state by the post-selection measurement induced by the second conformational change, the nucleotide ends up in a mixture of S and P as below:

where b^′ and d^′ are the re-normalized b and d, as before. If we feed master equation (2.9) with this incoherent superposition state instead of the coherent superposition state |ψ(t₄)〉_N given in (3.8), the steady-state solution ρ^eq_N|cHB obtained is the same density matrix with ρ^eq_N|qHB in (3.9). Hence, the equilibrium probabilities remain the same in both cases. Our conclusion still holds even if the quantum correlations of H-bonds found in the reaction intermediate ES are converted to classical correlations by strong decoherence just at the beginning of the catalytic stage.

(iii) Activation energy.

A state transformation from ρ to ρ^′ is said to be thermodynamically favourable only if the free energy F goes down, i.e. ΔF = F[ρ^′] − F[ρ] ≤ 0. Biological systems can carry out thermodynamically unfavourable state transformations by coupling them to favourable ones. However, even thermodynamically favourable biochemical state transformations generally involve the formation of transition states with higher free energies than initial states. Enzymes lower these activation barriers but are not expected to remove them completely. The difference in the free energy between the transition and initial (or final) states guarantees that this catalysed transformation is sufficiently slow for the organism inside which it takes place.

It is worth remembering that the tautomerization process of interest is not assumed to be the ultimate catalytic function of the multifunctional enzyme, but rather supports a secondary catalytic function required for the ultimate one. In this respect, an interesting question to ask is whether all enzyme-catalysed reactions have a rate-limiting step between the initial and final states. In other words, can an enzyme provide the formation of a transition state with a free energy lower than the free energy of the initial (or final) state in specific cases? The simple model investigated in this study points to a surprising possible answer to this question in what follows.

The free energy of the nucleotide can be calculated by:

where R is the universal gas constant, 1.987 cal (mol K)⁻¹. As the states of S, P and S^‡ are pure states, i.e. as tr[ρ²] = 1 for these tautomers, their von Neumann entropies vanishes. These states are also eigenstates of H_N which means that average energies taken over them are nothing but the corresponding eigenenergies. Hence, we end up with We consider the free energy of the transition state ES^‡ that appears after the binding stage of enzyme catalysis and study the case in which quantum correlations that are generated through the H-bonds of ES is preserved until the start of the catalytic stage, i.e. τ₂ < τ_D. |ES^‡〉 is separable (see equation (3.8)) and the reduced state of the nucleotide equals |ES^‡〉_N = |G〉_C⊗(b^′|01〉 + d^′|10〉)₁₂ where $d^{\mathrm{\prime }}=\sqrt{1-{b^{\mathrm{\prime }}}^2}$. The expected value of the energy over this state is tr[|ES^‡〉〈ES^‡|_NH_N] = b^′2ϵ₂ + (1 − b^′2)ϵ₁. As it is a pure state, its von Neumann entropy vanishes and its free energy becomes which corresponds to a free energy just in-between the free energy of S and P, a value much lower than expected in biochemical reactions (figure 4).

• Download figure
• Open in new tab
• Download PowerPoint

In the case of ${\tau }_2\nless {\tau }_D$, |ES^‡〉 is separable again (see equation (3.11)) but the reduced state of the nucleotide equals |ES^‡〉〈ES^‡|_N = |G〉〈G|_C⊗((1 − b^′2)|10〉〈10|₁₂ + b^′2|01〉〈01|₁₂). As tr[|ES^‡〉〈ES^‡|_NH_N] is also equal to b^′2ϵ₂ + (1 − b^′2)ϵ₁, the expected value of the energy between S and P along the reaction coordinate is not affected by decoherence. However, unlike the previous case, the reduced state of the nucleotide is not a pure state, which means that its von Neumann entropy does not vanish: S[|ES^‡〉〈ES^‡|_N] =  − b^′2log₂b^′2 − (1 − b^′2)log₂(1 − b^′2) > 0. The free energy then drops further by the loss of quantumness of correlations in the H-bonds as below:

However, this entropy-based decrease is not likely to be significant as RT < 1 kcal mol⁻¹ at physiological temperatures and S[|ES^‡〉〈ES^‡|_N] ≤ 1. Thus, F[|ES^‡〉〈ES^‡|_N] is expected to be in-between F[|S〉〈S|] and F[|P〉〈P|] once again.

We can numerically compare the activation energies in the catalysed and uncatalysed reactions using the model parameters that we chose above (figure 4). Let us fix F[|S〉〈S|] to zero for simplicity. To do so, both of the conformation and pseudo-spin terms of ϵ₁ are equated to zero, i.e. E_g = 0 and W₁ + λ₁₂ = 0. Furthermore, we assume without any loss of generality that ((W₁ − W₂)² + 4J²₁₂)^1/2 ≈ W₂ − W₁. Then we end up with F[|P〉〈P|] ≈ 14 kcal mol⁻¹ and F[|S^‡〉〈S^‡|] ≈ 42 kcal mol⁻¹. When we take T = 37.5°C and b^′2 = 1/2 (or equivalently $b=d=1/\sqrt{6}$) as before, F[|ES^‡〉〈ES^‡|_N] becomes 7 kcal mol⁻¹ if τ₂≤τ_D, and 6.38 kcal mol⁻¹ otherwise. When b^′2 is raised to 3/4, these values change to 10.5 and 10 kcal mol⁻¹, respectively. On the other hand, a decrease in b^′2 to 1/4 also lowers F[|ES^‡〉〈ES^‡|_N] to 3.5 for τ₂ < τ_D and 3 kcal mol⁻¹ for τ₂ > τ_D. Hence the rate enhancement the enzyme brings about is in the range of 22 to 27 orders of magnitude for 1/4 ≤ b^′2 ≤ 3/4 and T = 37.5°C, based on the first-order rate equation k = (k_BT/h) e^{(−ΔF‡)/RT}, where k_B is the Boltzmann constant and h is the Planck constant. The lower limit of this rate enhancement can be decreased down to 5 × 10¹⁹ by increasing the upper limit of b^′2 up to unity. Both limits can be pulled down by squeezing the energy levels of the Hamiltonian given in (2.5), e.g. when 1/4 ≤ b^′2 < 1, T = 37.5°C, E_e − E_g = 34 kcal mol⁻¹ and ϵ +  − ϵ −  = 19 kcal mol⁻¹, the enzyme leads to a rate enhancement in the range of 10–20 orders of magnitude. Note that although this enhancement is more reasonable when compared with the experimentally estimated upper limits, it does not change our main proposition that the free energy of the transition state is smaller than the free energy of the product in the enzyme catalysed reaction.