A novel mechanism for energy activation in biomolecules
1 A novel mechanism for energy activation in biomolecules
Shanshan Wu and Ao Ma* Department of Bioengineering The University of Illinois at Chicago 851 South Morgan Street Chicago, IL 60607 *correspondence should be addressed to: Ao Ma Email: [email protected] Tel: (312) 996-7225 Abstract
An activated process consists of energy activation and barrier crossing; the former is a prerequisite for the latter. Barrier crossing has been studied extensively, but energy activation has been overlooked due to a lack of means to gauge its progress. We define reaction stability π ! as the probability that reactive trajectories pass a vicinity in phase space; it enabled us to analyze energy activation of a biomolecular isomerization. This process follows a mechanism fundamentally different from presumed mechanisms in standard reaction rate theoriesβit features accumulation of high kinetic energy in reaction coordinates, achieved by precise synergy between them coordinated by momentum space. Activated process is essential to all molecular systems, ranging from chemical reactions of small molecules to conformational dynamics and enzymatic catalysis of proteins. For proteins, all the functionally important processes are activated processes because it gives them well-defined rates, which is critical for the functional roles of proteins in the cellular context because proper timing is key for functioning. The defining feature of an activated process is that the system needs to cross an activation barrier that is much higher than the thermal energy π " π , where π " is the Boltzmann constant and π is temperature. There are two critical questions concerning the mechanism of an activated process: (1) how the reaction coordinates acquire sufficient energy to cross the activation barrier, and (2) how the reaction coordinates cross the barrier once they acquired adequate energy. The first concerns energy activation and the second concerns barrier crossing [1]. The former is a prerequisite for and must precede the latter. Barrier crossing has been intensively studied both theoretically, through efforts in refining the reaction rate theories [1-3], and numerically, with molecular dynamics simulations. The complexity of numerical studies has grown with the increase in computational power: from early studies of gas phase reactions [4], to later studies of reactions in solution [5], to biomolecular dynamics in the past two decades with the help of transition path sampling method [6, 7] to harvest unbiased reactive trajectories and committor [8] to parameterize the progress of barrier crossing. In contrast, energy activation only received attention at the early stage of the development of reaction rate theories, with three lines of ideas. In transition state theory ( TST ), the problem of energy activation is bypassed by the quasi-equilibrium assumption between transition state and reactant. In Lindmanβs mechanism for unimolecular reactions in gas phase, it was assumed that collisions between a reactant molecule and buffer gas provide the reactant with extra energy, which is quickly equilibrated among all the coordinates via fast intra-molecular vibrational redistribution [1]. The subsequent barrier crossing of the high-energy reactant is under this condition of equipartition of energy over the entire molecule. In Kramers theory, it was assumed that the reaction coordinate acquires energy from thermal fluctuations of the solvents, modeled as a random force [9]. Therefore, the reaction coordinate can climb up the activation barrier during rare fluctuations in which the random force increases the total energy of the reaction coordinate. The overwhelming interest in activated processes in proteins over the past few decades has presented new challenges as well as opportunities. Since protein reactions always occur in solutions, TST and Kramers theory were widely adopted. However, a complex molecule like a protein differs fundamentally from a simple molecule, on which Lindmanβs and Kramersβ ideas were based. The large size of a complex molecule means it has sufficient degrees of freedom that it can provide the reaction coordinates with adequate energy for activation on its own. In contrast, simple molecules require an external energy source, be it buffer gas or solvents, for activation. More importantly, in complex molecules reaction coordinates and the thermal bath are connected by chemical bonds and the interactions between them involve a wide range of spatial scales. Bonded interactions impose complex and strong constraints on the motions of bath coordinates. This is in stark contrast with buffer gas molecules that essentially undergo free motions and solvent molecules that undergo a combination of libration and diffusion. The existence of interactions with different spatial scales presents a much more complex situation than the linear coupling between reaction coordinates and thermal bath assumed in standard reaction rate theories. These unique features of complex molecules make it possible that their energy activation may follow a mechanism fundamentally different from mechanisms suggested by Lindmanβs and Kramersβ ideas. This could be the reason behind the challenges encountered in applications of TST and Kramers theory to activated processes in proteins [10, 11]. One instance that hints at this issue is an unexpected observation from our previous study of a prototype of biomolecular isomerization dynamicsβthe πΆ β πΆ transition of an alanine dipeptide in vacuum [12]. Alanine dipeptide is the smallest molecule in which the non-reaction coordinates constitute a large enough thermal bath for activation, namely, it is the smallest example of complex molecules. Previous studies have identified dihedrals π and π ( (Fig. 1) as the essential reaction coordinates for this processβthey are sufficient for determining the value of committor π " [7, 13]. Committor is the probability that a dynamic trajectory initiated from a configuration space point, with momenta drawn from Boltzmann distribution, to be a reactive trajectory. It is the reaction probability in configuration space and provides a rigorous parameterization of the progress of barrier crossing. Our previous study showed that the activation barrier is located on the path of π , and π ( helps π to cross this barrier by directly transferring kinetic energy to π [12, 14]. The unexpected observation concerns failed attempts of π to cross the activation barrier. As shown in Fig. 2, π already reached the critical value that marks the onset of a successful transition at t = 1.03 ps. However, instead of moving forward and crossing the activation barrier, π reversed its direction and receded back to the reactant basin. A closer examination revealed that the critical difference between the failed and the succeeded barrier crossing is the position of π ( . The incorrect position of π ( makes the force from π ( to π (Fig. 2), defined as ΞπΉ ) ! β+ =β β« , " -./0β2,+,) ! ! ππ ( ( π/π 1β3 is the system potential energy), strongly repulsive, which drives π away from the activation barrier and reverses its direction of motion.. Even though the force that pushes π up the activation barrier is orders of magnitude stronger than what π normally experiences during a successful transition, ΞπΉ ) ! β+ increases much faster until it exceeds the force that facilitates π βs barrier crossing and reverses the direction of the total force acting on π . These results show that the motions of π and π ( need to be precisely coordinated: it appears that π ( acts as a gating mechanism on the motion of π . The observation above suggests that some critical events occurred in the region where π " = 0 ubiquitously, the region that π " cannot parameterize. This indicates a region where the momentum space is critical for activation [15], as π " provides rigorous parameterization of activation in configuration space. To understand the difference between a failed attempt for barrier crossing and a successful one, we need a rigorous parameterization of this region. The first step is to find a proper reference point. The natural choice is a phase space point Ξ on an existing reactive trajectory. If we perturb the momenta of Ξ slightly, we obtain a phase space point Ξ in its close vicinity. The probability that a dynamic trajectory launched from Ξ is a reactive trajectory is less than 1. This probability reflects the likelihood that reactive trajectories pass the vicinity around Ξ . Based on this idea, we can define a new parameter π ! , which we call the reaction stability. For a given phase space point Ξ on an existing reactive trajectory, we can generate an ensemble of phase space points πΈ(Ξ ) in its vicinity. A point Ξ in πΈ(Ξ ) is generated from a small random perturbation π (e.g. 20%) to the momenta of Ξ . From Ξ we can launch a dynamic trajectory and check if it is reactive. In this way, we obtain the probability of trajectories launched from πΈ(Ξ ) to be reactive, which is the value of π ! for Ξ . A large value of π ! suggests that Ξ lies in a phase space region that has dense populations of reactive trajectories, so that a perturbation to Ξ will more likely land on a phase space point on another reactive trajectory. Along a given reactive trajectory, the system in general moves from region with low density towards region with high density of reactive trajectories. Therefore, π ! provides a rigorous parameterization of the progress of activation in phase space. The region with π ! β (0, 1] and π " = 0 precedes the region with π " β (0,1] ; the former corresponds to the energy activation phase and the latter corresponds to the barrier crossing phase. To analyze the mechanism for energy activation, we use the energy flow theory we recently developed. Within this theory, we define both potential ( PEFs ) and kinetic (
KEFs ) energy flows during an activated process. The PEF through a coordinate π is its work: βπ (π‘ ( , π‘ ) = β C ππ(πβ)ππ ππ $ (3 " )% $ (3 ! ) (1). According to Eq. (1), Ξπ (π‘ ( , π‘ ) is the change in the potential energy of the system due to the motion of π along a dynamic trajectory in the time interval [ π‘ ( , π‘ ] . It is a projection of the change in the total potential energy onto the motion of π . Therefore, it is a measure of the cost of the motion of π in terms of potential energy. Accordingly, the change in the total potential energy of the system can be decomposed into PEFs through different coordinates: Ξπ(π‘ ( , π‘ ) = π(π‘ ) βπ(π‘ ( ) = β β Ξπ (π‘ ( , π‘ ) :6;( , where the summation is over all coordinates of the system. A major finding from our previous PEF analysis was that reaction coordinates are the coordinates with high PEFs during barrier crossing [12]. The KEF through a coordinate π is [14]: π < πΎ = ππΎππΜ ππΜ + KππΎππ L %0β % ,<0β ππ = ππ‘ Mπ πΜ + KππΎππ L %0β % ,<0β πΜ O (3), where πΎ is the system kinetic energy, πβ = (π ( , π , β― , π , π , β― , π : ) is the system position vector in internal coordinates with π removed, and π£β = (πΜ ( , πΜ , β― , πΜ : ) is the velocity vector. Since π < πΎ is the change in the system kinetic energy caused by changes in (π , πΜ ) alone, which fully describes the motion of π , it rigorously defines the KEF through π . Similarly, we have: ππΎ = β π < πΎ . To gain mechanistic insights, we need to look at how the PEFs and KEFs of individual coordinates change with the progress of activation. We first project the PEF or KEF onto a projector π(Ξ) that parameterizes the progress of activation, then average over the ensemble of reactive trajectories: β©πΏπ΄ (π β )βͺ = β« πΞπ(Ξ)πΏπ΄ (π(Ξ) β π(Ξ) + ππ)πΏ(π(Ξ) β π β )β« πΞπ(Ξ)πΏ(π(Ξ) β π β ) β©Ξπ΄ (π ( β π )βͺ = C β©πΏπ΄ (π)βͺ A " A ! (4) Here, is the probability of finding the system in an infinitesimal volume around a point Ξ in phase space in the transition path ensemble; πΏ(π₯) is the Dirac Ξ΄-function; πΏπ΄ (π(Ξ) β π(Ξ) + ππ) is the change in π΄ in a differential interval ; β©Ξπ΄ (π ( β π )βͺ is the change in π΄ in a finite interval [π ( , π ] , Ξπ΄ can be either Ξπ or Ξ < πΎ [12, 14]. For the barrier crossing phase, the optimal projector is π = π " ; for the energy activation phase, the optimal projector is π = π ! . Our results show that the duration of the energy activation phase is generally about 10 times longer than the duration of the barrier crossing phase, suggesting the former is a more complex and important process than the latter. Figure 3 shows the PEFs and KEFs through all the coordinates in the system. As expected, only the two dominant reaction coordinates, π and π ( , experience significant PEFs during energy activation. Moreover, the PEFs through π ( during energy activation and barrier crossing are of opposite signs. While π ( receives energy through PEFs during barrier crossing, there is a barrier of ~ β mol -1 on its path of motion during energy activation. Importantly, the energy flows in π cannot start until π ( reaches the top of this barrier at π ! β 0.4 . This explains why π ( can act as a gating mechanism on π : the necessary Ο ( Ξ ) d Ξ d Ξ ΞΎ ( Ξ ), ΞΎ ( Ξ ) + d ΞΎ β‘β£ ) condition for barrier crossing of π is its energy flow, which cannot start before π ( crossed its own barrier and became βreadyβ to help π , as shown in ref. [14]. A surprising observation is that kinetic energy plays the dominant role in energy activation for π . After π ( crossed the barrier on its path of motion, kinetic energy starts to accumulate in π until it reaches ~10 kJ β mol -1 , which is about the same as the total amount of potential energy consumed during its barrier crossing. This means that the energy cost for the barrier crossing of π is mainly paid by the kinetic energy it gathered during a long energy activation process. In addition, the momentum space is critical in this process, i.e. the proper alignment of momenta of all the coordinates in the system is critical for successfully building up kinetic energy in π . Moreover, a significant portion (~60%) of kinetic energy accumulated in π is due to the loss of kinetic energy in π ( and π . This also explains why the motions of reaction coordinates are highly coordinated as their momenta have to align properly to ensure successful transfer of kinetic energy from π ( and π to π . The mechanism of energy activation discussed above, the persistent accumulation of high amount of kinetic energy in the reaction coordinates, is fundamentally different from the physical picture for energy activation suggested by Lindmanβs and Kramers theories, which are the foundation of our understanding of activated dynamics. It suggests that energy activation in complex molecules like proteins follows fundamentally new mechanisms. Since alanine dipeptide is the smallest complex molecule, real proteins likely employ energy activation mechanisms that are more sophisticated and effective. This might be the reason behind the extraordinary efficiency of enzymes: proper energy activation is the foundation for protein functions and efficient energy activation leads to efficient functionality. Figure 1 : Two representative structures of an alanine dipeptide for the πΆ (solid color) and πΆ (semi-transparent color) states. The part of the molecule that does not change in the πΆ βπΆ transition completely overlaps between the two representations. The essential reaction coordinates π and π ( are indicated. ΟΞΈ C C Figure 2: (Upper) Time evolution of π (green) and π ( (red) along a reactive trajectory that includes a failed attempt of barrier crossing marked by an arrow. Blue dashed line: time evolution of π " . The horizontal dashed line marks the critical value of π that marks the onset of barrier crossing. The two vertical dashed lines mark the region of the successful barrier crossing. (Lower) Time evolution of the total force acting on π ( ΞπΉ + ; green) and the force from π ( to π ( ΞπΉ ) ! β+ ; red). Note the huge difference in the magnitudes of ΞπΉ + and ΞπΉ ) ! β+ between the failed attempt around t = 0.4 ps and the successful transition during π‘ β [1, 1.3] ps. -2-101 t (ps) -400-2000200 π , π ! ( r a d ) Ξ πΉ ! , Ξ πΉ " ! β ! ( kJ β m o l - β r ad - ) Failed Attempt
Success. Transit. Figure 3 : The PEFs ( β¨Ξπ β© π = π, π ( , π ; dashed lines) and KEFs ( β¨Ξ < πΎ β© π = π, π ( , π ; solid lines) during energy activation phase (left panel) and barrier crossing phase (right panel). The blue arrow marks the top of the barrier on the path of π ( during energy activation. The gray dashed lines are the PEFs and KEFs of the other coordinates (57 in total) in the system. p S -10-50510 ΞΈ Ο Ο p B β¨ Ξ π ! β© , β¨ Ξ " πΎ ! β© ( kJ β m o l - * Energy Activation
Barrier Crossing References
1. Berne, B.J., M. Borkovec, and J.E. Straub,
Classical and Modern Methods in Reaction-Rate Theory.
Journal of Physical Chemistry, 1988. (13): p. 3711-3725. 2. Pechukas, P., Statistical Approximations in Collision Theory , in
Dynamics of Molecular Collisions Part B , W.H. Miller, Editor. 1976, Plenum: New York. p. 269. 3. Berezhkovskii, A. and A. Szabo,
One-dimensional reaction coordinates for diffusive activated rate processes in many dimensions.
Journal of Chemical Physics, 2005. (1). 4. Karplus, M., R.N. Porter, and R.D. Sharma,
Exchange reactions with activation energy. I. Simple barrier potential for (H, H2).
J. Chem. Phys., 1965. (9): p. 3259-3287. 5. Montgomery, J.A., D. Chandler, and B.J. Berne, Trajectory Analysis of a Kinetic-Theory for Isomerization Dynamics in Condensed Phases.
Journal of Chemical Physics, 1979. (9): p. 4056-4066. 6. Bolhuis, P.G., Chandler, D., Dellago, C. and Geissler, P. L., Transition path sampling: Throwing ropes over rough mountain passes, in the dark.
Ann. Rev. Phys. Chem., 2002. : p. 291-318. 7. Bolhuis, P.G., C. Dellago, and D. Chandler, Reaction coordinates of biomolecular isomerization.
Proceedings of the National Academy of Sciences of the United States of America, 2000. (11): p. 5877-5882. 8. Ryter, D., On the Eigenfunctions of the Fokker-Planck Operator and of Its Adjoint.
Physica A, 1987. (1-3): p. 103-121. 9. Kramers, H.A.,
Brownian motion in a field of force and the diffusion model of chemical reactions.
Physica, 1940.
VII (4): p. 284-304. 10. Schwartz, S.D. and V.L. Schramm,
Enzymatic transition states and dynamic motion in barrier crossing.
Nat Chem Biol, 2009. (8): p. 551-8. 11. Garcia-Viloca, M., et al., How enzymes work: Analysis by modern rate theory and computer simulations.
Science, 2004. (5655): p. 186-195. 12. Li, W. and A. Ma,
Reaction mechanism and reaction coordinates from the viewpoint of energy flow.
J Chem Phys, 2016. (11): p. 114103. 13. Ma, A. and A.R. Dinner,
Automatic method for identifying reaction coordinates in complex systems.
Journal of Physical Chemistry B, 2005. (14): p. 6769-6779. 14. Li, H. and A. Ma,
Kinetic energy flows in activated dynamics of biomolecules. http://arxiv.org/abs/2007.06733 (submitted to J. Chem. Phys.), 2020. 15. Antoniou, D. and S.D. Schwartz,
Phase Space Bottlenecks in Enzymatic Reactions.
J Phys Chem B, 2016. (3): p. 433-9.(3): p. 433-9.