Roles of Local Non-equilibrium Free Energy in the Description of Biomolecules

When a system is in equilibrium, external perturbations yield a time series of non-equilibrium distributions, and recent experimental techniques give access to the non-equilibrium data that may contain critical information. Jinwoo and Tanaka (L. Jinwoo and H. Tanaka, Sci. Rep. 2015, 5, 7832) have provided mathematical proof that such a process's non-equilibrium free energy profile over a system's substates has Jarzynski's work as content, which spontaneously dissipates while molecules perform their tasks. Here we numerically verify this fact and give a practical example where we analyze a computer simulation of RNA translocation by a ring-shaped ATPase motor. By interpreting the cyclic process of substrate translocation as a series of quenching, relaxation, and second quenching, the theory gives how much individual sub-states of the ATPase motor have been energized until the end of the process. It turns out that the efficiency of RNA translocation is $48\sim 60\%$ for most molecules, but $12\%$ of molecules achieve $80\sim 100\%$ efficiency, which is consistent with the literature. This theory would be a valuable tool for extracting quantitative information about molecular non-equilibrium behavior from experimental observations.


I. INTRODUCTION
Counting is the basis of thermodynamics, which describes the behavior of biomolecules [1]. Free energy is (negative logarithm of) the number of configurations with more weights to lower energy conformations [2]. If one restricts to a specific sub-state, it gives the concept of conformational free energy, (negative logarithm of) the weighted number of conformations that belong to that sub-state. Thus when the second law of thermodynamics states that molecules tend to minimize their free energy, it means that molecules would be in (sub)states that include more accessible configurations with lower energies.
We may transform a system's (conformational) free energy by supplying energy through external control, such as in the case of single-molecule experiments [3][4][5][6]. Since molecules behave stochastically, we build an ensemble by repeating experiments to apply thermodynamics in understanding molecules' behavior [2,7,8]. For example, Jarzynski's work fluctuation theorem enables one to convert fluctuating work into the difference of free energies between the end states of external control [9].
On the one hand, Jarzynski's work fluctuation theorem can be very informative since free energy provides comprehensive information about the state of a system [6,[10][11][12]. On the other hand, it might be insufficient to provide detailed information on molecular nonequilibrium behavior on the level of sub-states [13][14][15]. For example, when a molecular machine hydrolyzes ATP, hydrolysis energy will probably increase the molecules' energy, allowing it to overcome the free energy barrier. However, since the process proceeds in an entire non-equilibrium situation and free energy provides only macro-information, a detailed description of this process is elusive. * jinwoolee@kw.ac.kr Jinwoo and Tanaka have recently revealed that those trajectories that reach the individual sub-state of a system contain essential details for the thermodynamics of each sub-state [16]. Significantly, they show that each substate's local non-equilibrium free energy has Jarzynski's work as content, allowing us to see how the introduced energy by external control affects molecular substates.
In this paper, we numerically verify Jinwoo and Tanaka's local version of Jarzynski's work fluctuation theorem. Mainly, we elucidate local non-equilibrium free energy, which contains detailed information about biomolecules performing their functions while allocating energy between sub-states and dissipating that energy.
As an application, we analyze a simulation carried by Ma and Schulten for RNA translocation by a ring-shaped ATPase motor [17]. There, we reveal Jarzynski's work required for reaching each substate at the end of the translocation process.
We organize the paper as follows: Section 2 introduces the theoretical framework while clarifying the terms used. In Section 3, firstly, we provide simulation results, confirming that local non-equilibrium free energy and an ensemble average of Jarzynski's work coincide at each time and substate. Secondly, we provide an application showing that the theory enables one to extract critical quantities from experiments or simulations. Section 4 discusses the implication of the results.
Boltzmann constant, and T is the heat bath temperature. External control λ t at time t of the system determines the system's free energy F (λ t ), (negative logarithm of) the number of configurations with more weights to lower energy configurations: where E(z; λ t ) is the internal energy of the system configuration z. Eq. (1) indicates that the lower the free energy, the richer the accessible configurations with lower energies. We note that external control λ t directly changes the internal energy E and the set of accessible configurations, Ω. If one considers free energy restricted to a sub-state Γ, it gives the concept of conformational free energy G(Γ), (negative logarithm of) the weighted number of configurations restricted to that sub-state: Eq.
(2) tells that the lower the conformational free energy of sub-state Γ, the richer the accessible microstates with lower energies. Thus with λ t fixed, after exploring the conformational free energy landscape for an extended period, the smaller the conformational free energy of substate Γ, the higher the probability that a system stays in Γ: We consider a situation where we pre-determinedly varying external control λ t as time t varies from 0 to τ . A system's microstate z t at time t would form a trajectory, {z t } 0≤t≤τ . If one accumulates along trajectory {z t } 0≤t≤τ the increment of internal energy E at each z t due to external control λ t defines Jarzynski's work: On the other hand, if one accumulates along trajectory {z t } 0≤t≤τ the increment of internal energy at each z t due to thermal fluctuation defines heat: where • indicates 'stochastic multiplication' in the Stratonovich sense. Combining Eq. (4) and Eq. (5) gives the first law of thermodynamics along trajectory {z t } 0≤t≤τ [30]: where ∆E = E(z τ ; λ τ ) − E(z 0 ; λ 0 ) indicates the increment of internal energy along the path.

B. Microscopic Reversibility
Since microscopic reversibility forms the basis of work fluctuation theorems, we discuss it briefly [31][32][33]. To this end, we consider the time-reversed process the same as filming the forward process and playing it backward. For control λ t that varies from time 0 to τ , let timereversed control λ t ≡ λ τ −t , where ≡ denotes 'defined by'. As time t goes from 0 to τ , time-reversed control λ t strictly follows λ t in a time-reversed manner. For trajectory {z t } 0≤t≤τ , let time-reversed trajectory {z t } 0≤t≤τ ≡ {z τ −t } 0≤t≤τ , which again follows the forward trajectory in a time-reversed manner as time t goes from 0 to τ .
The microscopic reversibility condition says that a lot of spontaneous heat flow from the heat bath to the system along a trajectory is less probable or a rare event. In detail, the conditional probability of time-forward path {z t } 0≤t≤τ conditioned on the initial configuration z 0 decreases exponentially as heat flow from the heat bath to the system along the path increases [34]: Here the forward probability is compared to the conditional probability of the time-reversed trajectory {z t } 0≤t≤τ conditioned on the initial conformation z 0 for the backward process. This comparison to the timereversed probability p makes the condition Eq. (7) invariant under time reversal. If one reversed time, the left-hand side of Eq. (7) would flip, and heat Q absorbed by the system from the heat bath would be released upon time reversal, flipping the right-hand side of Eq. (7), too. Thus, condition Eq. (7) states precisely the same for the time-reversed trajectory, so its name derives.

C. Work Fluctuation Theorems
Let us consider a single molecule (e.g., RNA) pulling experiment using an optical tweezer. We assume the system is in equilibrium at external control λ 0 [4,6,25,35]. As external control λ t pulls RNA during 0 ≤ t ≤ τ , it would give energy to RNA by W as in Eq. (4), changing the free energy landscape of the molecule (see Eq. (1) and Eq. (2), which depend on energy E). Repeating the experiment (which would hypothetically yield a nonequilibrium probability of sub-states Γ at each time t) and taking the following average of Jarzynski's work W would give the free energy difference between the two end states λ 0 and λ τ [33,36]: where · indicates the average over the repeated experiments, and F is free energy as defined in Eq. (1). Jarzynski's work fluctuation theorem, Eq. (8), takes into account all the paths generated while the forward process repeats. On the other hand, Jinwoo and Tanaka considers only those paths that reach a sub-state Γ at final time τ . In detail, they have shown that the following average of Jarzynski's work of paths that reach Γ at time τ gives local non-equilibrium free energy F of sub-state Γ as follows [16]: where · Γ,τ on the left-hand side indicates the average over the conditioned paths that reach Γ at time τ . On the right-hand side of Eq. (9) appears the difference between the local non-equilibrium free energy F of sub-state Γ at time τ and free energy F at external control λ 0 . One may interpret Eq. (9) as that the local non-equilibrium free energy of each sub-state at time τ has the ensemble of Jarzynski's work of each path that reaches that sub-state at time τ as content. Here, the local non-equilibrium free energy F of Γ at τ is composed of conformational free energy G of sub-state Γ as in Eq.
(2) and stochastic entropy, −k B ln p, as follows: where p(Γ, t) is the non-equilibrium probability of RNA being in a sub-state Γ at time t. We note that rewriting Eq. (10) gives Here the Boltzmann factor in the nominator is the same as the equilibrium distribution, Eq. (3). In the denominator appears local non-equilibrium free energy F of substate Γ at time t with Jarzynski's work as content, as in Eq. (9). Thus Eq. (11) represents that a sub-state Γ at time t endowed with energy through Jarzynski's work gets a higher probability of being realized, overcoming energy barriers due to the Boltzmann factor.

III. RESULTS
In order to verify the local version of Jarzynski's work fluctuation theorem Eq. (9), we carry out Brownian motion simulations by solving the over-damped Langevin equation [37,38]: which tells that a system's configuration is subject to a force in the direction of the lower energy during random walks caused by thermal fluctuations. Thermal fluctuation ξ increases with temperature and should be uncorrelated at different times, i.e., ξ(t) ξ(t ) = 2k B T ζδ(t − t ).
Here · indicates the average over all realizations and δ is the Dirac delta function. We set friction coefficient ζ = 1 and k B T = 2 and consider one-dimensional domain 0 ≤ z ≤ L with L = 10 with reflecting boundaries.
We partition the domain into 20 bins, forming sub-states {Γ i : i = 1, 2, · · · , 20}. The initial probability distribution p(z, 0) is uniform on the domain. We consider two different external controls. The first one is single quenching, and the second is complex scheduling, including continuous changing, quenching, and relaxation.

A. A quenched process
At time t = 0, we abruptly turn on bi-stable potential E(z) depicted in Figure 1c and solve Eq. (12), taking 4k time steps with step size δt = 0.01 for a single trajectory during 0 ≤ t ≤ 40. We repeat the experiment 200k times. The initial quenching would cause the initial equilibrium distribution into a non-equilibrium one and subsequently generate time series of non-equilibrium probability distributions p(z, t).
Panels (a1) and (a2) in Figure 1 show the time series of non-equilibrium free energy profiles. Initially, stochastic entropy (− ln p) is constant due to the uniformly distributed initial condition. Thus the conformational free energy G in Eq. (2), which resembles E(z) depicted in Figure 1c in our case, determines the shape of nonequilibrium free energy F in Eq. (10). Since a particle is subject to a force in the direction of the lower energy, particles tend to move towards local minima, z = L/4 and z = 3L/4 in the early rapid stage within 50 steps (see Panel (a1) in Figure 1). If stochastic entropy (− ln p) and conformational free energy G are in balance, nonequilibrium free energy F becomes locally flat. Then the process of balancing the non-equilibrium free energy on a larger scale between regions z < 2L/4 and z > 2L/4 (see Panel (a2) in Figure 1). This large-scale process proceeds very slowly during the remaining 3.5k steps because a single particle that experiences energy E(z) takes time in overcoming barrier region L/4 < z < 2L/4 in Figure  1c.
The color-coded plus symbols in Panels (b1) and (b2) of Figure 1 represent the ensemble average of Jarzynski's work on the left-hand side of Eq. (9) multiplied by e −βF (λ0) to match F. Plus symbols are overlaid with the profiles of the non-equilibrium free energy for comparison. When quenching occurs, the particles gain different energies depending on their position. Since quenching occurs only once, gained work of individual particles is constant throughout the process. Therefore, it is worth questioning how the ensemble average of this work tends to dissipate. It is because the particles occupying a specific sub-state at any time come from different places and are mixed. If one considers a sub-state near L/4 in the early stage, particles mainly located in region 0 ≤ z ≤ 2L/4 would gather in that sub-state and are mixed so that the ensemble average of Jarzynski's work over those particles would be equalized. This mixing will proceed over time, similar to how the non-equilibrium free energies find balance.
In Panel (a3) in Figure 1, the dashed blue line repre- sents the time series of the average over sub-states of the non-equilibrium free energy shown in Panels (a1) and (a2). The dotted red line shows the time series of the mean over sub-states of the work content represented in Panels (b1) and (b2). Both agree well and tend to dissipate.

B. A process including continuous changing, relaxation, and quenching
The scheduling of external control in the second process is as follows. At time t = 0, we abruptly turn on bi-stable potential E(z) depicted in Figure 1c. Then gradually diminishes its amplitude by multiplying a timedependent stiffness constant κ(t) that varies from 1 to 0 for 0 ≤ t < 20. At t = 20, we apply the second quenching by setting κ = 1/2 and keep the value during the remaining steps so that the system would relax. Figure 2c represents the resulting E(z; λ t ) = κ(t)E(z) for 0 ≤ t ≤ 20. We solve Eq. (12) in the same way as in the first process, taking 4k time steps with step size δt = 0.01 for a single trajectory during 0 ≤ t ≤ 40 and repeating the experiment 200k times. Figure 2a shows the time series of the non-equilibrium free energy profiles for 0.5 < t ≤ 40, omitting the initial relaxation stage that exhibits the same pattern as the case of the first process. In Panels (b1) and (b2) of Figure 2, the color-coded plus symbols (overlaid with F for comparison) represent the ensemble average of Jarzynski's work (on the left-hand side of Eq. (9) multiplied by e βF (λ0) ) for 1 < t ≤ 21. The selected interval includes the second quenching and thus exhibits patterns not seen in the first process. Data for final times in Panel (b1) of Figure 2 continues to the data for initial times in Panel (b2) of Figure 2. We see that F and work-content agree well during continuous changing, abrupt quenching, and the final relaxation. Panel (a2) of Figure 2 shows the averages over sub-states of F and work-content, and we see a convex pattern for 5 < t < 20, which is worth mentioning. It is because the relaxation speed of particles during the large-scale relaxation stage cannot keep up with the change speed of the external control. Due to this, particles near z = L/4 lose energy, and particles near z = 3L/4 gain energy but more quickly than the former, increasing the average F and balancing the two regions forcefully. Even after the forceful balance, the particles in the right-well keep gaining energy so that they have excess non-equilibrium free energy, which would be relaxed subsequently, decreasing the average F (see the final times of Panel (b1) of Figure 2).

A. Mechanism of RNA Translocation
Here we analyze Ma and Schulten's simulation data for RNA-translocation by Rho hexameric helicase [17]. Rho hexameric helicase [39], a ring-shaped motor, translocates the substrate (RNA) during the rotary reaction from the state R to S [17]. The state R represents a stable state taken from the crystal structure [39] with the ATP-mimics (ADP·BeF3) at the six active sites being replaced by the equivalent ATP molecules (at four sites) and ADP+Pi (at one site), leaving the remaining one site empty. The state S is identical to R except 60 • clockwise rotation. The transition R → S occurs through subunitsubunit interface conformational changes, during which RNA is propelled.
The transition R → S of the motor, which is not spontaneous, contains a dwell phase of a relative long dura-tion and a motor-action phase [40]. During the dwellphase, the motor is energetically charged through three step processes: (1) ATP binding in the empty site, (2) release of existing (say "old") ATP hydrolysis product (ADP+P i ), and (3) hydrolysis of ATP into "new" product (ADP+P i ) [17]. The subprocess (2) is considered rate limiting [41,42]. Thus the subprocesses (1,3) take place at the very beginning of the transition, Rho is charged energetically, and waits until the subprocess (2) occurs. This state defines I (see Figure 3a). Then, subprocess (2) initiates the motor-action phase during which Rho translocates RNA through subunit-subunit interface conformational changes, which completes the cycle.
We interpret this cyclic process as a series of quenching, relaxation, and second quenching. In detail, we prepare the initial state as the state that subprocesses (1,3) have just taken place. Then, we take the first quenching that realizes subprocess (2), translocating RNA during relaxation. Then, we take the second quenching that realizes subprocess (1,3), which would turn the system into the initial state after another relaxation. We do not involve this final relaxation and compare the initial probability distribution and the final non-equilibrium caused by the second quenching.

B. Interpretation
We introduce hypothetical external control λ t that turns on and off pre-determinedly the interaction energy between Rho hexameric helicase, ATP, and ADP+P i . For an initial state, we set λ t (for t < 0) as the state where the subprocesses (1,3) have just taken place from the state R. In detail, regarding subprocess (1), the interaction energy between ATP and an active site of Rho is turned on. For subprocess (2), the interaction energy between the existing ATP hydrolysis product (ADP+P i ) and Rho is kept turning on, describing the state where the "old" product is not released. Concerning subprocess (3), We change the interaction between an ATP and Rho to the interaction between ADP+P i and Rho, describing the hydrolysis of ATP into a "new" product (ADP+P i ). Let us shorten this initial state as A, i.e., λ t = A for t < 0.
Then, the equilibrium probability distribution p 0 of this state reads: where G is the conformational free energy of sub-state Γ i at state A, and F is the free energy at state A. Ma and Schulten [17] defined reaction coordinates Γ as a collection of positions of key residues at the six subunit-subunit interfaces, which contributes significantly to the relative motion of subunits. The transition R → S is then partitioned into 51 groups in terms of Γ so that we have Γ i = 0, 1, · · · , 50, which we adopt as the sub-states. They calculated G(Γ i ; A) − F (A), which we digitized to obtain p 0 (Γ i ; A). We normalized so that 50 0 p 0 (Γ i ; A) = 1. In Figures 3b and 3c, blue curves represent the conformation free energy and the initial probability distribution given λ t = A.
At time t = 0, we abruptly turn off the interaction energy concerning subprocess (2) so that the "old" ATP hydrolysis product (ADP+P i ) is released. Let us denote this quenched state by B so that λ t = B for 0 ≤ t < τ . This quenching would initiate a motor-action phase, yielding the time-series of non-equilibrium probability distributions p: where F(Γ i , t) in the denominator is the non-equilibrium free energy of sub-state Γ i at time t. We do not change λ t until Rho translocates RNA during relaxation. Let τ be the time long enough to complete the process. We note that τ should be fixed during the repeat of this (hypothetical) experiment. After the relaxation, the final equilibrium probability distribution p 1 would be: Ma and Schulten [17] calculated G(Γ i ; B) − F (B), which we digitized to obtain p 1 (Γ i ; B). We normalized so that 50 0 p 1 (Γ i ; B) = 1. In Figures 3b and 3c, orange curves represent the conformation free energy and the final probability distribution given λ t = B. At time t = τ , we perform the second quenching by setting λ τ = A, the same state as the initial one. Then the final equilibrium distribution turns into a non-equilibrium one: We apply the local work fluctuation theorem to this process λ t for 0 ≤ t ≤ τ to obtain We can obtain the right-hand side of Eq. (17) by comparing the initial probability distribution Eq. (13) and the final non-equilibrium probability distribution Eq. (16) as follows: which is represented in Figure 3d.

C. Analysis
Let us follow what happens to the ensemble of rho hexameric helicases (with RNAs) during external control λ t varies from 0 to τ . The initial probability distribution with λ t = A (blue curve in Figure 3c) tells that most instances would be near Γ i = 5. If the molecules' "old" ATP hydrolysis products are released, or λ t is set to B, the conformational free energy profile would change from the blue curve in Figure 3b to the orange one so that the ensemble of conformations near Γ i = 5 would have excess non-equilibrium free energy. As time flows, the excess free energy would become flat locally, covering regions 0 < Γ i < 28. Subsequently, global balancing would proceed very slowly between regions Γ i < 28 and Γ i > 28. This stage corresponds to individual molecules' overcoming the free energy barrier near Γ i = 28. As the time approaches τ , the non-equilibrium free energy would become flat for all Γ i , and most instances would be near Γ i = 41 as the orange curve in Figure 3c indicates. Then, we set λ t = A; then again, they would have excess free energy, which is shown in Figure 3d.
According to the local work fluctuation theorem Eq. (17), the excess free energy also tells that how much the molecules have been energized during the process, especially by the two quenching steps. We note that λ t = A is the state where an ATP is already hydrolyzed before the process begins at time t = 0. During the intermediate 0 < t < τ , the release of the "old" ATP product causes a motor action. Only at time t = τ , when we set λ τ to A, the hydrolysis of new ATP happens. So the net effect of the two quenching upon a single instance of the system is the hydrolysis of one ATP molecule.
Analysis of the simulation result for RNA translocation by Rho hexameric helicase: (a) The state R represents a stable state taken from the crystal structure [39] with the ATP-mimics at the six active sites being replaced by the equivalent ATP molecules (at four sites) and ADP+Pi (at one site), leaving the remaining one site empty. The state S is identical to R except for 60 • clockwise rotation of subunit-subunit interface conformations (which are schematically represented by different colors), which is achieved by conformational changes of subunit-subunit interfaces along the reaction coordinate. This rotary reaction R → S is energetically prepared during a dwell phase: (1) ATP binding at the empty site and (2) ATP hydrolysis in another already occupied (by ATP) site [17]. This state defines state I. The motor reaction phase is then initiated by the release of ATP-bound product (ADP+Pi) so that the stored energy during the dwell phase causes the rotary reaction. At time τ , it is the most probable that the system is in state Γ i ≈ 41 and its non-equilibrium free energy value, F(Γ i = 41, τ ) − F (A), is approximately 12k B T . Provided that ATP hydrolysis energy is 20 ∼ 25k B T , depending on the environmental condition [43], the efficiency of the motor action corresponds to 48 ∼ 60%. It is interesting that the non-equilibrium free energy contains work values for another states at time τ . For example, sub-state Γ i = 45, the value of ∆F(Γ i = 45, τ ) is about 20k B T , corresponding to 80 ∼ 100% efficiency. The final probability value for this state is p 1 (Γ i = 45, τ ) = 0.12, indicating that 12% of the molecules in the ensemble are achieving this efficiency. We note that some literature reports near 100% efficiency of molecular motors [44].

V. CONCLUSIONS
We have verified the local version of the work fluctuation theorem linking Jarzynski's work and local nonequilibrium free energy. In various conditions of external controls, including quenching, relaxation, and continuous changing, the non-equilibrium free energy and Jarzyn-ski's work agree very well for each substate at any time.
As an application, we analyzed a simulation result of RNA translocation by Rho hexameric helicase. By treating one cycle of a rotary reaction mechanism as a series of quenching, relaxation, and second quenching, we could obtain Jarynski's work content of each substate of Rho at the end of the translocation process. This theory would enable one to extract non-equilibrium work content from modern single-molecule experiments and simulations.

VI. APPENDIX
A. Proof of Eq. (9) We consider the local non-equilibrium free energy f for micro-state z: f (z, t) = E(z; λ t ) + 1 β ln p(z, t), Let {Γ i } be sub-states. We calculate the following average of local non-equilibrium free energy for microstates where we used Eq. (2), Eq. (10), and Eq. (19). We will use this fluctuation theorem linking f and F to prove the local work fluctuation theorem. Let us consider external control λ t from 0 to τ . We assume that the initial probability for the reversed process is the same as the final probability of the forward process, i.e., p (z 0 , 0) = p(z, τ ). By the property of conditional probability, we have p [{z t } 0≤t≤τ ] = p [{z t } 0≤t≤τ |z 0 ]p(z 0 , 0) and similarly, p [{z t } 0≤t≤τ ] = p [{z t } 0≤t≤τ |z 0 ] p (z 0 , 0).