Nonequilibrium steady state of biochemical cycle kinetics under non-isothermal conditions

Nonequilibrium steady state of isothermal biochemical cycle kinetics has been extensively studied, but much less investigated under non-isothermal conditions. However, once the heat exchange between subsystems is rather slow, the isothermal assumption of the whole system meets great challenge, which is indeed the case inside many kinds of living organisms. Here we generalize the nonequilibrium steady-state theory of isothermal biochemical cycle kinetics, in the master-equation models, to the situation in which the temperatures of subsystems can be far from uniform. We first obtain a new thermodynamic relation between the chemical reaction rates and thermodynamic potentials under such a non-isothermal circumstances, which immediately implies simply applying the isothermal transition-state rate formula for each chemical reaction in terms of only the reactants' temperature, is not thermodynamically consistent. Therefore, we mathematically derive several revised reaction-rate formulas which not only obey the new thermodynamic relation but also approximate the exact reaction rate better than the rate formula under isothermal condition. The new thermodynamic relation also predicts that in the transporter system with different temperatures inside and outside the membrane, the net flux of the transported molecules can possibly even go against the temperature gradient in the absence of the chemical driving force.

The temperature gradient in living cells has been experimentally observed and shown to have a nonnegligible impact on many important biological functions [13][14][15]. This phenomenon has recently attracted significant interest from both cell biologists and physical chemists [16][17][18]. The isothermal assumption is only valid when the timescale of the biochemical kinetics is much slower than the timescale of heat exchange within the system. If the system under investigation is divided into several subsystems, and the heat exchange among them is sufficiently slow compared to the biochemical kinetics being investigated, the temperatures, which may be non-uniform, can only be well defined for each subsystem [13,19]. For example, in living cells [20,21], a transporter protein crossing the cell membrane interacts with both the intracellular and extracellular components of the cell, which may be at different temperatures [16]. A mesoscopic chemical kinetic system such as this is driven not only by chemical potential differences but also by the temperature gradient across the cell membrane [17].
Under isothermal or non-isothermal circumstances, mesoscopic biochemical cycle kinetics can be modeled using the master equation model. The theory of stochastic thermodynamics in terms of the master equation model has already been well developed, and the entropy production rate is expressed in terms of the transition rates between discrete states [4,11,12]. After taking the ensemble average, the mesoscopic stochastic thermodynamics should be consistent with the macroscopic nonequilibrium thermodynamics framework [22,23]. Such consistency under isothermal conditions has been illustrated previously [1,4,[24][25][26].
The problem of transport processes across a thermal gradient has been considered before in [16,17]. These references discussed the transport of ions across a plasma membrane, when the temperature inside the cell is different from the temperature in the extracellular medium. However, the temperature at the membrane interface where the molecular transporter occupies is assumed to be independent of the reaction coordinates. This assumption is quite uncertain and evades the circumstances when the assumption is not valid.
In our paper, we remove this assumption, and conduct a detailed and comprehensive theoretical analysis of the nonequilibrium steady state (NESS) of biochemical cycle kinetics under non-isothermal conditions. Starting from a four-state model of a transporter protein across the cell membrane, we derive a new thermodynamic relationship between the reaction rates in the master equation model and the thermodynamic potentials of the discrete chemical states involved based on the consistency between mesoscopic stochastic thermodynamics and macroscopic nonequilibrium thermodynamics. Although the temperature difference across the cell membrane is shown to be quite small (i.e., approximately 1 K [27]), we demonstrate that the entropy production can vary by 3% in real cells, revealing the significant contribution of non-uniform temperature to thermodynamics.
This thermodynamic relationship is not satisfied if one simply applies the isothermal transition-state rate formula for each reaction considering only the reactant' temperature. Instead, we mathematically derive several revised reaction rate formulas that not only obey the new thermodynamic relationship but also approximate the exact reaction rate better than isothermal calculations. Since we do not assume the independence of temperature upon the reaction coordinates, the temperature at the transition state of each reaction appears in our rate formulas, otherwise it will not appear as those in [16]. As long as one reaction is only assigned by one temperature, the transition rates between reactants and products have to follow the isothermal transition-rate formula under the same temperature reservoir, and the temperature at the transition state is irrelevant, just as the setup in [28,29]. However, in our setup, the reactants and products are assigned with different temperatures. So a new theory is in highly demand.
Additionally, we decompose the total thermodynamic driving force into its thermal and chemical components and predict that the net flux of the molecules transported by the molecular transporter can go against the temperature gradient in the absence of a chemical driving force, which does not violate the Second Law of Thermodynamics. Fluctuation theorems are also derived in the four-state model under non-isothermal conditions.
We then generalize the obtained results to a more complex six-state transporter model and general master equation models based on the existing cycle-decomposition theory of master equation models [30]. We derive the same new thermodynamic relation for each kinetic cycle of the system and validate its consistency with the revised reaction rate formulas. Multiple fluctuation theorems are also derived for each kinetic cycle.

A four-state transporter model under non-isothermal conditions
We start from a four-state model of transporter proteins across the cell membrane (figure 1). The nonequilibrium analysis of this model under isothermal conditions can be traced back to Hillʼs pioneering work [1,2]. E and E* are defined as the two conformations of the transporter protein that face the inside and outside of the cell membrane, respectively. The intracellular temperature is T 1 and the extracellular temperature is T 2 . When a counterclockwise cycle of the reaction diagram in figure 1 is completed, the transporter converts one molecule of M i on the inside of the cell to one molecule of M o on the outside.
Thermodynamic equilibrium means both thermal equilibrium T 1 =T 2 and chemical equilibrium (i.e., the chemical potentials inside and outside the membrane are equal to each other M Here, the chemical potential should be replaced with the electrochemical potential gradient, including the membrane potential effect, if the transported molecule is not electrically neutral.

Entropy production during the completion of a counterclockwise cycle and a new thermodynamic relationship
We consider the master-equation model of the transporter protein across the membrane, where k ij is the first-order or pseudo-first-order reaction rate from state i to state j, and p t i ( ) is the probability of state i at time t.
Therefore, in the steady state, during the completion of a counterclockwise cycle in figure 1, the log ratio of the forward and backward path probability, i.e. the entropy production, is The entropy balance equation in nonequilibrium thermodynamics tells that e p =ΔS−ΔS e , where ΔS is the change in the entropy of the transporter system (the molecular transporter as well as the molecules of M i and M o inside and outside the membrane) in figure 1 after the counterclockwise cycle, and ΔS e is the change in the entropy of the medium due to heat dissipation.
The resulting change in the entropy of the transporter system during this cycle is the entropy difference between molecules M o and M i ; i.e., and the change in the entropy of the medium due to heat dissipation is S e To perform a more detailed analysis of heat dissipation Q 1 and Q 2 , we consider the barrier-crossing picture for each reaction (figure 2). Transition state C lies at the saddle point of the potential energy surface along reaction coordinate x, while the reactant A and product B are located around the two local minima. In figure 2, we choose enthalpy as the potential energy, as is typically used for chemical reactions.
Here we assume that the reaction is under overdamped condition, i.e. the kinetic energy is neglectable, hence during a single reaction A B  , the change of internal energy before the reaction escapes the potential well of chemical state A in figure 2 is just H C −H A , which is the heat absorbed from the heat bath on the left-hand side of the transition state with temperature T 1 according to the First Law of Thermodynamics. Similarly, the heat released to the other heat bath on the right-hand side of the transition state with temperature T 2 is H C −H B . Here, H i denotes the average enthalpy of the chemical state i, where i=A, B, C.
Therefore, during the completion of a counterclockwise cycle in figure 1, the heat released to the two heat baths with temperatures T 1 where Q ĩ is the heat released to the heat bath with temperature T i during the regeneration process performed by the external device. Hence, the minimal entropy change of the medium due to heat dissipation during the regenerating process becomes which exactly represents the meaning of entropy production at steady state: all energy is dissipated as heat. Under steady-state conditions, the ensemble-averaged entropy production rate of the four-state model is [9,30] J k e epr , 10 in which J c is the cycle current, i.e. the averaged occurrence of the counterclockwise cycle per unit time.

The transmembrane flux of transported molecules
According to equation (4), the entropy production during a counterclockwise cycle in figure 1 can be rewritten as Figure 2. The enthalpy surface H(x) with two local minima A and B, where x is the reaction coordinate. Escape from A to B occurs at rate k AB . We assume that the temperature of the region on the left-hand side of X c is T 1 , and that the temperature on the right-hand side is T 2 . The temperature varies in the region around X c , which is called the connecting region.
Such a counterintuitive phenomenon cannot occur in a single reaction. Consider the single chemical reaction A B  in figure 2, where A is at temperature T 1 , and B is at temperature T 2 . The thermodynamic relationship of this single chemical reaction satisfies (see appendix A) since H AB * −H A is always positive. Thus, the net flux of molecules in non-cyclic kinetics always follows the temperature gradient in the absence of chemical driving forces.

The influence of small temperature differences on entropy production
It is thought that natural cell membranes, which have a thickness of ∼10 nm, cannot maintain a large temperature difference between their surfaces. Thus, temperature difference is likely to be small for a molecular transporter across a membrane and is estimated to be 1 K [27]. Such a small difference may have a nearly negligible influence on the dynamic/kinetic properties of the system. However, when considering thermodynamics, the results can be rather different. Figure 3 presents measured intracellular and extracellular concentration data for Na + and K + ions and membrane potentials [32][33][34][35].
Molecular dynamics simulations have already suggested that the activation energy for a transporter across the membrane is very high [33][34][35][36][37], so the enthalpy difference H H H figure 3 is assumed to be k T 20 48 kJ mol B » / , estimated from [37]. The entropy production can vary up to 3%, even if the temperature difference is only 1 K, provided the difference in the activation energies of the transmembrane kinetics for the molecular transporter bound or unbound with a small transported molecule is relatively high. When the temperature difference is approximately 3 K, the variation in entropy production can reach 5%-10%.

Fluctuation theorems
The stochastic net number of occurrences ν(τ) of the counterclockwise cycle in figure 1 up to time τ has been shown to satisfy the detailed fluctuation theorem  for any real λ. Therefore, for any physical quantity associated with the counterclockwise cycle (i.e., W(τ)=C·ν(τ) for any constant C), we can write for any w in which w/C is an integer, and e e . 1 9 For instance, in the four-state transporter model in figure 1, there are at least three thermodynamic quantities that we are interested in: e p , F chem , and F thermo . The fluctuation theorems described by equations (18) and (19) are applicable for each thermodynamic value.

A more general form of the new thermodynamic relationship
In the high-dimensional case, Barezhkovskii et al [38] showed that a one-dimensional reaction coordinate exists for any system whose transition rate is described by Langerʼs multidimensional generalization of the onedimensional Kramers' theory of diffusive barrier crossing. We can derive a more general form of the new thermodynamic relation (equation (5)) by regarding the temperature as a function of the one-dimensional reaction coordinate without further assumption or simplification.
During the completion of a counterclockwise cycle in figure 1, the entropy production can be represented as (see appendix B for detailed derivation) Notice that for the isothermal case with only one heat bath, the first term on the right-hand side of equation (21) vanishes, and hence, equation (21) reduces to (6). When the effect of the connecting region on the reaction rates is nearly negligible, equation (21) is reduced to (5).
3. Thermodynamically consistent revised reaction rate formulas 3.1. Simple application of the isothermal transition-state rate formula The well-known transition-state rate formula is derived under isothermal conditions [39,40]; that is, D ‡ ‡ and ΔS ‡ are the free energy, enthalpy, and entropy of activation, respectively. Applying the isothermal transition-state rate formula for each reaction in terms of the reactants' temperature, the following relationship can be derived (see appendix C): . S ij ‡ is the entropy of the transition state along the reaction coordinate from state i to state j.
There is no theoretical support for the consistent vanishing of the Δ term; thus, it is contradictory to the new thermodynamic relation (equation (5)).
The term Δ in equation (23) emerges because the transition states of the forward and backward reactions are attributed to different temperatures and are, thus, not reasonable. The transition state should have its own temperature T*, and a revised reaction rate formula is needed for the non-isothermal condition. Only when the temperature ν for each reaction is independent of reaction coordinate, the isothermal transition-state rate formula can apply. And only in this case, the local detailed balance condition for each reaction holds: in which j is the grand potential, also called Landau potential. However, in this paper, we drop this assumption, thus we need to seek for revised reaction-rate formulas.

The revised reaction rate formulas
The reaction rate from reactant A to product B crossing through transition state C as shown in figure 2 can be defined as the reciprocal of the mean first passage time from one local minimum x A of the enthalpy surface to the other local minimum at x B [41].
Consider the overdamped Langevin dynamics with the potential energy surface H(x) [39]: where the temperature T is a function of x. The stochastic integral should be interpreted as an Ito integral. The justification for this choice is due to the convergence theorem from the underdamping case to the overdamping case because this issue does not rise for the original second-order Langevin equation [42]. When the frictional coefficient is not independent of the position, the temperature-dependent component remains in the Ito form, while the frictional coefficient is in the anti-Ito form [43], also because of the local fluctuation-dissipation relation.
It is worth mentioning that the more general form of the new thermodynamic relation can be reproduced in the overdamped Brownian dynamics (see appendix E), in which the entropy production is defined as the log ratio of forward and backward path probabilities in the framework of stochastic thermodynamics.

One-dimensional case
In the one-dimensional case, we assume that near x A , x B , and x C ( figure 2), When the thermal fluctuation is rather small, we use the following approximation (denoted as k AB 1 ) for the transition rate from state A to state B (see appendix F): If the connecting region in figure 2 is narrow (see appendix F), equation (27) can be simplified to the following approximation (denoted as k AB 2 )(see appendix F): The conditions under which the revised rate formulas equations (27) and(28) approximate the exact reaction rate are summarized in table 1, and numerical validation is presented in figures 4(a), (b). Each approximation has its own range of applicability, and the most appropriate approximation for each region is always better than Kramers' approximation (figures 4(a), (b)).
Under isothermal conditions, the revised reaction rates in equations (27) and (28) are equivalent and the same as the established Kramerʼs rate formula. Whether these rate formulas are thermodynamically consistent will be verified using the more general high-dimensional versions below.

High-dimensional case
Without loss of generality, we denote the one-dimensional coordinate existing in the high-dimensional case as x and the other coordinates as y. Note that xäR 1 , and yäR n−1 . x R y y , : ) }is a set of parallel planes perpendicular to the x coordinate. We assume that, for fixed x, (x, 0) is the minimum of the potential energy H(x, y) on the plane {(x, y) : yäR n−1 }. We further assume that, around the reactant region x=x A , and that around the transition state region is the eigenvalue along the x direction. When the thermal fluctuation is rather small, we can use the following approximation for the transition rate from state A to state B (see appendix G): which is the high-dimensional version of equation (27).  , and that used in the high-dimensional case is Then, the entropy production after the completion of a counterclockwise cycle in figure 1 is (see appendix H) which is exactly the general form of the thermodynamic relationship (equation (21)), indicating that equation (30) is thermodynamically consistent. If we assume that the connecting region around the transition state x C is narrow and let T * =T(x C ), equation (30) can be simplified to which is the high-dimensional version of equation (28) and satisfies equation (29). However, even under the local Gaussian approximation, the rate formula equation (32) is not thermodynamically consistent with equation (5) if the dimension is not one (see appendix H). Therefore, we need a new approximation for equation (30).
Under the local Gaussian approximation, we derive another approximation of equation (30) (see appendix G): Notice that in calculating H C , the reaction coordinate x is not used, and H A involves one more coordinate than H C [40].
We rewrite the revised transition rate formula equation (33) in terms of the thermodynamic quantities which is exactly the new thermodynamic relation (equation (5)), showing that equation (33) is thermodynamically consistent. The three approximations k k , AB AB 1 2 , and k AB 3 are equivalent under isothermal conditions, which is exactly Kramers' rate formula. The first revised rate formula k AB 1 (equation (30)) holds even when the connecting region is large and is consistent with the thermodynamic relation equation (21). The second revised rate formula k AB 2 (equation (32)) is a modified version of Kramers' rate formula (equation (29)) but is not consistent with the thermodynamic relation equation (21) or (5). When the connecting region is narrow, the third revised rate formula k AB 3 (equation (33)) achieves a better approximation than the isothermal Kramers' rate formula (figures 4(c), (d)) and is consistent with the new thermodynamic relationship (equation (5)). The conditions under which the three revised formulas hold and their thermodynamic consistency are summarized in table 2.    for any real λ. Therefore, for any physical quantity associated with cycle i, where i=a, b, c, denoted as W(τ)=C·ν i (τ) for any constant C, we can write

A more complicated co-transporter model
for any w in which w/C is an integer, and e e . 4 6

General master equation model
Consider a general master equation model for N states {1, 2, L, N}, where k ij is the reaction rate from state i to state j, and p t i ( ) is the probability of state i at time t. Suppose that the temperature of state i is T i and that the enthalpy of the transition state along reaction i j  ) , we assume that after the completion of this cycle, there are n i molecules of chemical species S i at temperature T S i being converted into m j molecules of chemical species P j at temperature T P j . Then, the new thermodynamic relation for the cycle, which is similar to equation (5) is the cycle affinity. This is consistent with the revised rate formulas equations (28) and (33). We can also obtain a more general form of the new thermodynamic relation similar to equation (21) k which is consistent with the revised rate formula equation (30). Under steady-state conditions, the dynamics of the master equation can be broken down into cycles. The ensemble-averaged entropy production rate is [9,30]  for any real λ. Therefore, for any physical quantity associated with cycle c denoted as W(τ)=C·ν c (τ) for any constant C, we can write for any w in which w/C is an integer, and e e . 5 4

Discussion and conclusion
In summary, we have introduced temperature differences as an additional driving force of the chemical potential difference in biochemical cycle kinetics. Under non-isothermal conditions, the thermodynamic relationship between the reaction rates and thermodynamic potentials should be modified. Our approach is based on the consistency of the macroscopic nonequilibrium thermodynamics and mesoscopic stochastic thermodynamics, and assumptions of the separation of time scales between the internal dynamics (e.g., the dynamics of the molecular transporter), the dynamics of its surrounding environment, and the separation of time scales among the dynamics within each subsystem and those across their boundaries. We have derived several revised rate formulas for a single chemical reaction under non-isothermal conditions, most of which are consistent with the new thermodynamic relationship and approximate the exact reaction rate better than the isothermal transitionstate rate formula. The thermodynamic analysis presented here also suggests that the transmembrane flux of molecules in the cyclic transporter model can even go against the temperature gradient in the absence of a chemical driving force.
In [18,28,29,46], they also consider mesoscopic systems in contact with multiple temperature baths. But we have very different local equilibrium assumptions: in their approach, the different temperatures is assigned to different reactions (i.e., parallel reaction paths with different temperatures for each chemical reaction) and local equilibrium (detailed balance) for each reaction is assumed; while in our approach, the different temperatures is assigned to different chemical states, and we assumed local Einstein relation along the continuous transition path. Their assumption is unrealistic if we describe the biochemical cycle kinetics inside a living cell using the master equation approach. The temperature in our master equation model defined for chemical states rather than reaction transitions forces us to consider more details about the reaction path, including the enthalpy and temperature of the transition state.
However, different parameter regions need different approximations, even within the same parameter region, different approximations are still possible. Typically, any of these approaches can be applied in the parameter region in which they are valid, but if we are also interested in the associated thermodynamics, certain existing approximations that are not consistent with thermodynamic laws are not suitable. For examples, reaction-rate formulas with simple expressions are always approximations of the exact values of reaction rates. No matter how good the approximation is, its thermodynamic consistency should be considered, when applying the rate formulas to study nonequilibrium thermodynamics of a real biochemical system.
Lastly, the master equation is a mathematical model developed to describe stochastic biochemical dynamics. If a temperature gradient is present, the disappearance of the entropy production rate of the master equation is not equivalent to the thermodynamic equilibrium: both the chemical and thermal driving forces vanish. The inconsistency between the vanishing of the entropy production rate and thermodynamic equilibrium under the non-isothermal conditions can be attributed to the neglecting of kinetic energy relating to temperature in this model. In a description without kinetic energy, the thermal driving force can hardly be distinguished from other driving forces when applying the mathematical framework of stochastic thermodynamics [43,[47][48][49]. For this reason, we carefully studied the heat dissipation along the reaction coordinates of the chemical reactions in the present work.
We assume A is under heat bath T 1 , and B is under heat bath T 2 , when a single reaction A B  is complete, the heat released to the two heat baths with temperature  Q  T  H  H  T  T  T  H  H  S T  S T  1  1 On the other hand, we may apply equation (32) in the main text for the forward and backward reactions: Notice that the temperature difference is small, ΔT/T=1, so the term k ln

Appendix B. More general form of the new thermodynamic relation
If we consider the temperature as a function of reaction coordinate, we can derive a general thermodynamic relation.
For a single reaction in one dimension (see figure 2), we assume that around x A and x B , ]. Under the overdamped assumption, i.e. the vanishing of kinetic energy, from the time that the trajectory escapes from [x 0 , x i ] to the time escaping from [x 0 , x i+1 ], the change of internal energy is H(x i+1 )−H(x i ), which is exactly the heat absorbed from the heat reservoir with temperature T(x i ). Therefore Then, after taking ensemble average, we have The third equality use the fact that the temperature is T 1 around x A and is T 2 around x B ; the last equality use the l . In the high dimensional case, under the local Gaussian approximation, we first let H(x) denotes the average enthalpy of (x, 0) on the rest (n−1) dimensions that perpendicular to the x coordinate:

Appendix D. Size of the connecting region
We denote the connecting region around the transition state as D, and the region between the two local minima of the enthalpy surface as Ω, then the statement 'narrow connecting region' we said in the main text is equivalent to the inequality in which the temperature T is a function of x, and the local Einstein relation holds. u(x)=E x τ B is the mean first passage time to B, start at x. Then u(x) satisfied the following Dynkin equation .
Assume additionally that f x 0 , then Because the ratio between k B T(x) and the barrier height along the enthalpy surface is rather small, we can get the following asymptotic behavior of the reaction rate: We define We assume that around x x , A B and x C , Recalling the relationship between m x ( ) and Therefore, we can have the following approximation for the transition rate from state A to state B when the fluctuation is rather small, If we assume that in most of the region on the left-hand side of x C , the temperature is always T 1 , while in most of the region on the right-hand side of x C , the temperature is T 2 , and the temperature varies slowly and slightly in a small region around x C , then we have: 16) can be simplified to: Appendix G. Reaction-rate formulas under non-isothermal condition in n dimensions n 1 > ( ) As in the Barezhkovskiiʼs paper [38], we consider the potential of mean force along the reaction coordinate x, the potential of mean force is defined as: where the length L is given by: Recall that in calculating H C the reaction coordinate x is left out, and H A involves one more coordinate than H C [40].   F.16) and (G.6)) for the forward and backward reaction rates e.g. If we apply the revised reaction rate formula (equation (G.12)) for the forward and backward reaction rates, e.g.