The role of non-equilibrium fluxes in the relaxation processes of the linear chemical master equation Luciana Renata de Oliveira, Armando Bazzani, Enrico Giampieri, and Gastone C. Castellani Citation: The Journal of Chemical Physics 141, 065102 (2014); doi: 10.1063/1.4891515 View online: http://dx.doi.org/10.1063/1.4891515 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/141/6?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Role of filament annealing in the kinetics and thermodynamics of nucleated polymerization J. Chem. Phys. 140, 214904 (2014); 10.1063/1.4880121 Mechano-chemical coupling in Belousov-Zhabotinskii reactions J. Chem. Phys. 140, 124110 (2014); 10.1063/1.4869195 Non-invasive estimation of dissipation from non-equilibrium fluctuations in chemical reactions J. Chem. Phys. 139, 124109 (2013); 10.1063/1.4821760 Nucleated polymerization with secondary pathways. II. Determination of self-consistent solutions to growth processes described by non-linear master equations J. Chem. Phys. 135, 065106 (2011); 10.1063/1.3608917 Fluctuation theorem for nonequilibrium reactions J. Chem. Phys. 120, 8898 (2004); 10.1063/1.1688758
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.189.205.30 On: Wed, 10 Dec 2014 20:35:55
THE JOURNAL OF CHEMICAL PHYSICS 141, 065102 (2014)
The role of non-equilibrium fluxes in the relaxation processes of the linear chemical master equation Luciana Renata de Oliveira, Armando Bazzani, Enrico Giampieri, and Gastone C. Castellania) Physics and Astronomy Department, Bologna University and INFN Sezione di Bologna, Italy
(Received 28 March 2014; accepted 16 July 2014; published online 11 August 2014) We propose a non-equilibrium thermodynamical description in terms of the Chemical Master Equation (CME) to characterize the dynamics of a chemical cycle chain reaction among m different species. These systems can be closed or open for energy and molecules exchange with the environment, which determines how they relax to the stationary state. Closed systems reach an equilibrium state (characterized by the detailed balance condition (D.B.)), while open systems will reach a non-equilibrium steady state (NESS). The principal difference between D.B. and NESS is due to the presence of chemical fluxes. In the D.B. condition the fluxes are absent while for the NESS case, the chemical fluxes are necessary for the state maintaining. All the biological systems are characterized by their “far from equilibrium behavior,” hence the NESS is a good candidate for a realistic description of the dynamical and thermodynamical properties of living organisms. In this work we consider a CME written in terms of a discrete Kolmogorov forward equation, which lead us to write explicitly the non-equilibrium chemical fluxes. For systems in NESS, we show that there is a nonconservative “external vector field” whose is linearly proportional to the chemical fluxes. We also demonstrate that the modulation of these external fields does not change their stationary distributions, which ensure us to study the same system and outline the differences in the system’s behavior when it switches from the D.B. regime to NESS. We were interested to see how the non-equilibrium fluxes influence the relaxation process during the reaching of the stationary distribution. By performing analytical and numerical analysis, our central result is that the presence of the non-equilibrium chemical fluxes reduces the characteristic relaxation time with respect to the D.B. condition. Within a biochemical and biological perspective, this result can be related to the “plasticity property” of biological systems and to their capabilities to switch from one state to another as is observed during synaptic plasticity, cell fate determination, and differentiation. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4891515] I. INTRODUCTION
The dynamics of a chemical cycle chain reaction among m different species is relevant for a large number of biochemical systems, because they are fundamental for the kinetic description of biochemical networks; both in the case of small and large network cycles. These reaction cycles are the chemical basis of several cellular signaling processes such as signal transduction,1, 2 biological morphogenesis,3, 4 and cell fate determination.5, 6 These cycles are also relevant in the induction of synaptic plasticity, where a small network of phosphorylation/dephosphorylation reactions is taking place7–9 at the post-synaptic receptors level. How to describe these systems is a very important topic: their dynamics can be well approximated with nonlinear differential equations,9, 10 but in the case of small number of molecules we have to take into account for fluctuations and we need to apply stochastic methods such as the Chemical Master Equation (CME). The reaction cycles are the simplest non-trivial example of dynamical networks for which it is possible to develop an analytical approach. The characterization of the stationary distribution of a) Electronic mail: [email protected]
0021-9606/2014/141(6)/065102/11/$30.00
a CME of a reaction network is still a challenge, because one should deal with the fact that these systems can be open systems (here we understand open, as the possibility of exchanging molecules and energy with the surroundings).11–21 After a sufficient long time, isolated systems reach a thermodynamics equilibrium state (detailed balance (D.B.) equilibrium), where all the chemical currents are null. While an open system (if the exchange of energy with its surroundings is sustained) approaches a non-equilibrium steady state (NESS),13, 19–28 the probability chemical fluxes are not null leading to the breaking down of the D.B. When the D.B. holds, the equilibrium distribution can be written as a Maxwell-Boltzmann distribution, which means there is a unique equilibrium constant (the equilibrium temperature) for every chemical reaction in a system, regardless of how complex the system is. However, in a NESS, the stationary chemical currents are different from zero and this can be associated with the existence of an external non-conservative field interacting with the system11, 19, 25, 26 and with a positive entropy production due to the currents.29 As a matter of fact most of the chemical cycles relax toward NESS so that the problem of understanding the dynamical and thermodynamical properties of NESS could suggest if the energy cost of NESS with re-
141, 065102-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.189.205.30 On: Wed, 10 Dec 2014 20:35:55
065102-2
de Oliveira et al.
spect to the D.B. equilibrium is compensated by other advantages. In this paper we discuss a general approach to study the dynamics of a chemical chain reaction among m different species, based on the results of previous works,11, 13, 19–21, 24, 25 that allows to compute a drift and a diffusion coefficient and to introduce a non-equilibrium thermodynamical description for the stationary states of CME. The determination of stationary fluxes allows us to completely describe systems in equilibrium (D.B.) as well as in NESS.11, 15, 19, 21, 28 We discuss the role of chemical currents in the transient states and, in particular, in the relaxation process towards the stationary distribution. Indeed the relaxation times of biochemical reactions could be related to the plasticity properties of biological systems, where plasticity means the capability to reach the stationary distribution or to jump from one stationary state to another. The paper is organized as follows: Sec. II is devoted to the biochemical relevance of the unimolecular reaction cycle; in Sec. III we described the model and we discuss the meanfield approximation; in Sec. IV the thermodynamical formalism is introduced to characterize the stationary solutions; in Sec. V we explicitly study the effect of chemical fluxes on the relaxation time toward the stationary states in Sec. VI some numerical results are presented to check the analytic and the thermodynamically assumptions.
II. BIOCHEMICAL RELEVANCE OF LINEAR CYCLES
The more general, and probably realistic, form of the Michaelis-Menten reaction is obtained by allowing a certain degree of reversibility in the product formation step. This consideration leads to a model that was originally developed by Haldane in 1930. This mechanism can be used to model enzymatic reactions, and the kinetic scheme is as follows: This chain of reactions can be represented as a cycle between enzymatic states (E, ES, EP). In Fig. 1 the free enzyme E, after binding the substrate S, forms the first reaction intermediate, the enzyme-substrate complex (ES). This complex, in turn, is transformed in the second reaction intermediate, the enzyme-product complex (EP), where the substrate has been transformed in the product P. As final step, the product is released and the enzyme turns back into the free form. The ES complex can dissociate back in the free enzyme E and the unbounded substrate S, and, in the same way, the association reaction between the free enzyme E and the product P can form the EP complex. If the abundance of the substrate and the product are in chemical balance with the environment, the transition rates of the cycle can be approximated by pseudofirst order rate constants (see Fig. 2). From a general point of view, the possible states that the enzymatic complex can assume would be more numer-
FIG. 1. Reversible Michaelis-Menten reaction scheme. The transformation between substrate and product is composed of a chain of reversible steps.
J. Chem. Phys. 141, 065102 (2014)
FIG. 2. Reversible Michaelis-Menten reduced chain reaction. Under the right conditions the transformation can be approximated as a cycle of linear reversible transformations.
ous than two, but are usually still representable as a cyclic enzymatic transformation. A detailed analysis of the thermodynamical properties of these kind of systems has been done in Ref. 14. This kind of cyclic reactions can be founded in several molecular mechanisms underlying important biological processes such as cell differentiation, cell fate determination, and synaptic plasticity. The bidirectional synaptic plasticity (the transition between long-term synaptic potentiation (LTP) and long-term synaptic depression (LTD) ) is generally assumed as the basis of learning and memory. The molecular mechanism is a phosphorylation/dephosphorylation of sites on the postsynaptic glutamate receptor α-amino-3-hydroxy-5-methyl4-isoxazolepropionic acid receptor (AMPA), on the subunit GluR1. This subunit can be phosphorylated by specific protein-kinase enzymes (such as CaMKII, protein kinase C and A). The induction of LTP specifically increases phosphorylation of one site, while the induction of LTD is accompanied by a decrease in the phosphorylation of the other site. Therefore, knowledge of the phosphorylation state of GluR1 may be a strong predictor of synaptic plasticity. Moreover, the phosphorylation state of AMPA receptors in the same subunit is affecting their cytoplasmatic trafficking, specially their insertion/removal from the postsynaptic membrane (the double phosphorylated AMPA receptor is more expressed in the cellular membrane, while the non-phosphorylated is mainly in the cytoplasm). This phosphorylation/dephosphorylation cycle is represented schematically in Fig. 3. The fraction of AMPARs containing GluR1 phosphorylated at one site is denoted by A∗ , whereas A∗ denotes the fraction of AMPARs containing GluR1s phosphorylated on another, non equivalent, site, and A and A∗∗ denote AMPARs phosphorylation at neither and both sites, respectively. EP1 and EP2 are used to denote the protein phosphatases that dephosphorylate these two sites, and EK1 and EK2 are used to denote the protein kinases that phosphorylate these sites. This molecular cycle, if we assume low concentration of substrates (the number of AMPARs, that in a single spine are on the order of hundred molecules) with respect to the saturation limit of the enzymes (the Michaelis-Menten constant Km ), is linear as the cycle in Fig. 2 and the same considerations can be used for both systems. Such a system is, by definition, always far from thermodynamical equilibrium, and the organism has to continuously
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.189.205.30 On: Wed, 10 Dec 2014 20:35:55
065102-3
de Oliveira et al.
J. Chem. Phys. 141, 065102 (2014)
where n ∈ N m and Ek± are the Van Kampen step operators defined according to Ek± f (n) = f (n1 , n2 , .., nk ± 1, .., nm )
n ∈ N m,
for any function f(n). If the reaction chain is a cycle, we impose periodic boundary conditions in the sum m + 1 ∼ 1 and the total number of particles is constant |n| =
m
nk = N.
(2)
k=1
FIG. 3. Representation of phosphorylation/dephosphorylation cycle of AMPA receptors. Under the right conditions the transformation between the receptors’ phosphorylation state can be approximated by a cycle of four linear reversible transformations driven by two couples of kinases and phosphatases.
In a general case the transition probability π k − 1, k may depend on the distribution n of particles into the different chemical species. B. The drift and diffusion coefficients for the discrete case
spend energy to maintain the system in this state and to ensure the capability to switch between different phosphorylation states.
If one introduces the difference operators Dk± defined as + Dk+ ≡ Ek−1 Ek− − 1
− Dk− ≡ Ek−1 Ek+ − 1,
(3)
Eq. (1) can be written as a discrete Kolmogorov forward equation11, 15, 19, 20 (see Appendix A): III. NON-EQUILIBRIUM FLUXES AND STATIONARY STATES FOR THE CME
˙ t) = − p(n,
Let us consider the dynamics of a chemical chain reaction among m different species (see Fig. 4). Under the assumption that each particle is independent the transition rate from the specie k − 1 to the specie k is given by π k − 1, k (n)nk − 1 where nk − 1 denotes the number of particles of the k − 1 specie and π k − 1, k is the probability that a single particle performs a transition in a time unit. Then we write the CME for the evolution of the probability distribution p(n, t) according to22 ˙ t) = p(n,
+ [(Ek−1 Ek−
Dk+ Jk
k=1
A. Definition of the model
m
− 1)πk−1,k (n)nk−1 p(n, t)
m
Dk+ [Ak (n)pn (t) + Dk− Bk (n)pn (t)]. (4)
k=1
In the continuous limit, Eq. (4) becomes a Fokker-Planck equation and the difference operators tend to the partial differential operators on the tangent space of the constraint (2). The vector Jk denotes the chemical currents and we have introduced the drift field Ak (n) = −πk−1,k (n)nk−1 + πk,k−1 (n)nk and the diffusion coefficient
k − Ek+ − 1)πk,k−1 (n)nk p(n, t)], + (Ek−1
=−
(1)
Bk (n) = πk,k−1 (n)nk . The drift field is directly correlated to the average dynamics: let {n ∈ N m / |n| = N } such that the chemical fluxes (or the distribution p(n, t)) almost vanish at its boundary then we have the average equations (see Eq. (A11) in Appendix A) n˙ k = Ak+1 (n) − Ak (n)
k = 1, .., m
(5)
with an error of the order of the chemical flux at the boundary of . We say that the mean-field approximation can be applied to Eq. (5) on , if n˙ k Ak+1 (n ) − Ak (n )
j = 1, .., m.
(6)
This approximation is correctly applied when the fluctuations are small with respect to the average values, and it is exact when the transition probabilities π k − 1, k are constant (i.e., the particles are independent). C. Solution of the linear CME FIG. 4. Chemical chain reaction among m different species. Each transition is regulated by its own transition rate, that, in general, depends on the number of particles in that state.
In the last case we get a linear CME whose timedependent solution can be explicitly computed in the form
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.189.205.30 On: Wed, 10 Dec 2014 20:35:55
065102-4
de Oliveira et al.
J. Chem. Phys. 141, 065102 (2014)
of a multinomial distribution,30–32 assuming that the system starts from a known state (p(n∗ , t = 0) = 1 for some n∗ ), p(n, t) = N !
n m λ k (t)
|n| = N,
k
k=1
nk !
(7)
Appendix A we show that the time derivative of the Gibbs’ entropy for the system can be written in the form33–35 S˙in =
k=1 |n|=N
|λ| =
λk = 1
λk > 0.
k=1
Letting t → ∞ we get the stationary distribution ps (n) with limt→∞ λk (t) = λ∗k where λ∗k satisfies the condition Ak+1 (N λ∗ ) = Ak (N λ∗ ). n∗k
When the average field approximation holds the values Nλ∗k , this corresponds to the local maximum of the stationary distribution ps (n) if the diffusion coefficient can be considered constant. Therefore, in the linear case the critical value of the stationary distribution corresponds to the stable fixed point of the average systems (6). IV. THERMODYNAMICAL PROPERTIES OF CME
A thermodynamical approach has been proposed to characterize the properties of the stationary solutions of CME.24 In particular one distinguishes the equilibrium states at which the chemical fluxes Jks (n) cancel out (D.B. condition for equilibrium states), and the NESS11, 13, 19, 20, 22–24, 26–28 where the weaker condition holds (see Eq. (4)), m
Dk+ Jks (n) = 0.
(8)
k=1
We will explicitly recall some properties of the stationary solution for Eq. (4). When the D.B. holds, the equilibrium distribution can be written as a Maxwell-Boltzmann distribution
Jk2 (n) Bk (n)pn (t)
Ak (n) + Dk− Bk (n) , (12) Jk (n, t) − Bk (n) k=1 |n|=N m
where the quantities λk are the solutions of the linear system (6) with the constraint m
m
when the stationary distribution satisfies Dk+ pn (t) 1 (continuous approximation) and it is possible to linearize the discrete derivative Dk+ ln pn (t) Dk+ pn (t)/pn (t). In the rhs of Eq. (12) one recognizes the total entropy production S˙ and the environment entropy production S˙en according to S˙ =
m k=1 |n|=N
S˙en
Jk2 (n) , Bk (n)pn (t)
(13) Ak (n) + Dk− Bk (n) . =− Jk (n, t) Bk (n) k=1 |n|=N m
In this way we get the mathematical formulation for the thermodynamic variables in terms of the master equation. When the system presents detailed balance condition, for t s s = S˙ s = S˙en = → ∞ it reaches an equilibrium state with S˙in s s 0. In contrast systems in NESS, S˙ and S˙en are equal but not s = 0. null, to ensure that S˙in Therefore, in a NESS it is straightforward to observe that the total entropy production is always positive, so that to maintain the stationary distribution the environment is exchanging energy with the system through the work done by the external forces, which is dissipated. Conversely in a D.B. equilibrium the total entropy production is zero and there is no dissipated work. By using the decomposition (10) it is possible to modulate the external field by changing the drift term according to Aˆ k (n) = Ak (n) − ηAex k (n),
p (n) ∝ exp(−V (n)), s
(9)
where V (n) is an interaction microscopic energy. The expectation value of V (n) can be identified with the internal energy of the system and it is natural to introduce the Gibbs’ entropy to describe the heat exchange with the environment.24 In a NESS it is possible to split the drift field into an internal and an external vector field,11, 14, 19–21, 25 Ak (n) =
Ain k (n)
+
Aext k (n),
(10)
where Aext k (n) is an external non-conservative field which generates the chemical fluxes Aex k (n) =
Jks (n) ps (n)
(14)
where η is a parameter: η = 0 corresponds to the initial case and η = 1 to the D.B. equilibrium when the external field vanishes. Moreover, one proves (see Appendix A) that the stationary distribution does not depend on η, whereas the new fluxes read Jˆks (n) = (1 − η)Jks (n).
(15)
Therefore, from this point of view it seems to be not convenient for any biochemical system to relax towards a NESS due to its energetic and entropic cost. Then the question is to study the effect of chemical fluxes in the transient states and, in particular, in the relaxation process towards the stationary distribution.
(11)
and Ain k (n) is a conservative vector whose potential satisfies V in (n) = − ln ps (n). The external field Aex k (n) is proportional to the thermodynamical force associated to the chemical fluxes. As it is known19–21, 25 in a NESS, the stationary fluxes Jks (n) are related to the total entropy production. In
V. NON-EQUILIBRIUM FLUXES IN THE LINEAR CME
We study the relaxation times of a linear CME in D.B. equilibrium and NESS to understand the influence of the nonequilibrium fluxes. For seek of simplicity we have chosen a chemical reaction with three states, as represented in Fig. 5
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.189.205.30 On: Wed, 10 Dec 2014 20:35:55
065102-5
de Oliveira et al.
J. Chem. Phys. 141, 065102 (2014)
FIG. 5. Unimolecular chemical reaction cycle. The cycle is composed only of reactions with a transition rate that is linear in the number of particles.
FIG. 6. Stationary distribution for CME (16) computed using (17) for N = 50 as a function of nA and nB .
Eq. (16) reads − ˙ x , ny , t) = πAB (E+ p(n x Ey − 1)nx p(nx , ny )
whose CME (see Eq. (1)) is
+ + πBA (E− x Ey − 1)ny p(nx , ny )
˙ t) p(n, =
− (E+ C EA
− 1)πCA nC p(n, t) +
− × (E+ A EB − × (E+ B EC
+ (E− C EA
+ πCA (E− x − 1)(N − nx − ny )p(nx , ny )
− 1)πAC nA p(n, t)
− 1)πAB nA p(n, t) +
+ (E− A EB
− 1)πBA nB p(n, t)
− 1)πBC nB p(n, t) +
+ (E− B EC
− 1)πCB nC p(n, t),
+ πCB (E− y − 1)(N − nx − ny )p(nx , ny ) + πAC (E+ x − 1)nx p(nx , ny ) + πBC (E+ y − 1)ny p(nx , ny ).
(16) where nA , nB , and nC are the number of molecules of the species A, B, and C, respectively, with the constraint nA + nB + nC = N. The particle distribution is a multinomial distribution (see (7)), n
n
n
λ A (t)λBB (t)λCC (t) , p(n, t) = N! A nA !nB !nC !
(17)
where λA (t), λB (t), λC (t) are the solutions of the linear system
(20)
We write CME (20) in the form of a continuity Eq. (4) without any constraint, ˙ x , ny ) = −Dx Jx (nx , ny ) − Dy Jy (nx , ny ) p(n
(21)
by using the discrete difference operators Dx = E+ x −1
Dy = E+ y − 1.
So we identify the fluxes Jx (nx , ny ) = JA (n) − JB (n) and Jy (nx , ny ) = JB (n) − JC (n) as − Jx (nx , ny ) = −(πAB E− y nx − πBA Ex ny + πAC nx
λ˙ A (t) = −(πAB + πAC )λA (t) + πBA λB (t) + πCA λC (t) λ˙ B (t) = πAB λA (t) − (πBA + πBC )λB (t) + πCB λC (t)
− πCA E− x (N − nx − ny ))p(nx , ny ), (18)
λ˙ C (t) = πAC λA (t) + πBC λB (t) − (πCA + πCB )λC (t)
Substituting the explicit form of the steady state solution pns x ,ny in (22), we determine the stationary non-equilibrium currents Jxs (nx , ny ) and Jys (nx , ny ), πBA πAB = − s N − nx − ny + 1 λsA λB λs + Cs πCA nx − πAC nx ps (nx , ny ) λA (23) λsC nx ny πAB πBA s − s Jy (nx , ny ) = N − nx − ny + 1 λsB λA λs + Cs πCB ny − πBC ny ps (nx , ny ) λB
λsA ∝ πBA πCA + πBC πCA + πBA πCB λsC
(22)
− πCB E− y (N − nx − ny ))p(nx , ny ).
with the constraint λA + λB + λC = 1. The steady state solution ps (n) is computed taking the limit t → ∞ in Eq. (17) and the limit values λs can be explicitly computed as follows:
λsB ∝ πAB πCA + πAB πCB + πAC πCB
− Jy (nx , ny ) = −(πBA E− x ny − πAB Ey nx + πBC ny
Jxs (nx , ny ) (19)
∝ πAC πBA + πAB πBC + πAC πBC .
Fig. 6 represents the stationary distribution of CME (16), for N = 50 as a function of nA and nB : the maximal value corresponds to nA = N λsA and nB = N λsB . For an explicit computation of the chemical fluxes we reduce the system’s dimensionality using nA = nx ; nB = ny and nC = N − nx − ny . Then
λsC nx ny
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.189.205.30 On: Wed, 10 Dec 2014 20:35:55
065102-6
de Oliveira et al.
J. Chem. Phys. 141, 065102 (2014)
According to Eq. (4), we write the non-equilibrium fluxes introducing a drift term and a diffusion coefficient Jx (nx , ny ) = Ax (nx , ny )p(nx , ny ) − Dx Bxx p(nx − 1, ny) − Dy Bxy p(nx , ny − 1),
(24) Jy (nx , ny ) = Ay (nx , ny )p(nx , ny ) − Dx Byx p(nx − 1, ny) − Dy Byy p(nx , ny − 1).
B=
+ πCA (N − nx − ny )
(25)
Ay (nx , ny ) = πAB nx − (πBA + πBC )ny + πCB (N − nx − ny )
πBA ny + πCA (N − nx − ny + 1)
−πAB nx
−πBA ny
πAB nx + πCB (N − nx − ny + 1)
=
(26)
s Aex y (nx , ny )p (nx , ny ).
+Jys (nx , ny + 1)Dy ps (nx , ny ) = 0.
(27)
λsC nx ny
+
Aex y (nx , ny ) =
N − nx − ny + 1 +
To compute the relaxation time towards the stationary distribution we take advantage from the fact that a multinomial distribution is completely determined by the average values. Then we study the average dynamics, n˙ x = −(πAB + πAC )nx + πBA ny + πCA (N − nx − ny ), n˙ y = +πAB nx − (πBA + πBC )ny + πCB (N − nx − ny ).
(28)
The relaxation process towards the limit values is an exponential whose exponents are the eigenvalues of Eq. (28). To study the effects of the chemical fluxes we modulate the external field which produces the fluxes by changing the drift term according to Aˆ x,y (n) = Ax,y (n) − ηAex x,y (n) (see Eq. (14)). In this way we change the non-equilibrium fluxes without modifying the stationary distribution (see Appendix A). From Eq. (15) the new non-equilibrium fluxes read
Jˆy (nx , ny , t) = Jy (nx , ny , t) − ηAex y (nx , ny )p(nx , ny , t).
πAB πBA s − λB λsA
λsC πCB s − πBC ny . λB
(30)
(31)
It is convenient to set η = 1 + , so that = 0 corresponds to the D.B. equilibrium and changing , we drive the system to a NESS condition. From (21) we get the modified CME, ˙ x , ny , t) p(n
N λsA , N λsB
Jˆx (nx , ny , t) = Jx (nx , ny , t) − ηAex x (nx , ny )p(nx , ny , t),
πBA πAB s − λA λsB
λsC πCA s − πAC nx , λA
λsC nx ny
VI. RESULTS
N − nx − ny + 1
In Appendix B we show that the fluxes Jxs (nx , ny ) and Jys (nx , ny ) are orthogonal to the discrete gradient of the stationary distribution ps (nx , ny ), Jxs (nx + 1, ny )Dx ps (nx , ny )
.
where we explicitly have Aex x (nx , ny )
s Jxs (nx , ny ) = Aex x (nx , ny )p (nx , ny ),
=
Ax (nx , ny ) = −(πAB + πAC )nx + πBA ny
and the diffusion matrix B is
Then we define an external field Aex (cf. Eq. (11)) related to the non-equilibrium fluxes as
Jys (nx , ny )
The drift vector A(nx , ny ) is
(29)
= −Dx (Jx (nx , ny , t) − (1 + )Aex x (nx , ny , t)p(nx , ny , t)) − Dy (Jy (nx , ny , t) − (1 + )Aex y (nx , ny , t)p(nx , ny , t)). (32) We remark that the new CME is non-linear. From the explicit knowledge of the stationary chemical currents, we are able to control the NESS condition for systems with = 0, using the s conditions Jks = 0 and S˙ s = S˙en
= 0. To confirm the NESS, we calculated the value of the norm of the non-equilibrium currents (||Js ||) as a function of (see Figure 7(a)) and we s estimate the values of S˙ s and S˙en for 0 ≤ ≤ 1 in the stationary state (see Figure 7(b)). In this way we verify for all s are not null and that values of = 0 that ||Js ||, and S˙ s = S˙en the system is in a NESS state.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.189.205.30 On: Wed, 10 Dec 2014 20:35:55
065102-7
de Oliveira et al.
J. Chem. Phys. 141, 065102 (2014) Ss 300 250 200 150 100 (a)
(b)
50 0.2
0.4
0.6
0.8
1.0
Ε
FIG. 7. (a) Norm of the non-equilibrium currents (||Js ||) as a function of , (b) Entropy production rate in the stationary state (Ss ) as a function of . The calculation is performed using N = 50 and π AB = 1; π CA = 1.1; π BC = 1; π BA = 1; π AC = 1; π CB = 1.
To study the relaxation process of CME (32) when the chemical fluxes change we consider the average equations,
∂Aex x (nx , ny )
n n˙ x = − πAB + πAC + πCA + (1 + ) x
∗ ∂nx nx
∂Aex x (nx , ny )
n , + πBA − πCA − (1 + ) y
∗ ∂ny ny (33)
∂Aex y (nx , ny )
n + n˙ y = πAB − πCB − (1 + ) x
∗ ∂nx nx
∂Aex y (nx , ny )
n − (πBA + πBC + πCB + (1 + ) y
∗ ∂ny ny in the mean field approximation. The evolution of the average values n(t) relaxes towards the stable critical point of Eqs. (33) with a characteristic time τ that depends on the minimum of the real part of eigenvalues α k k = 1, 2 computed linearizing the equation at the critical points, τ=
1 . Min|Re(α1 , α2 )|
(34)
The nonlinear character of the external fields Aex implies that the average Eqs. (33) are not exact, but we expect that the mean field approximation holds in the neighborhood of the
FIG. 8. Comparison between the relaxation time for the distance p(nx , ny , t) − ps (nx , ny ) for the solution of CME (32) and the relaxation time estimated by (34) using the average Eqs. (33). The parameters used in the simulations are p(nx, ny, 0) = 1/N2 ,π AB = 1; π CA = 1.1; π BC = 1; π BA = 1; π AC = 1; π CB = 1, which correspond to a NESS case with = .9, and N = 50 molecules. The black line is the norm p(nx , n)y, t) − ps (nx , ny ) and the gray line gives an interpolation with the exponential decaying cexp ( − t/τ ).
maximum of the stationary distribution function ps (nx , ny ) for the modified CME (32) that is well approximated by the stable critical point of the average Eq. (33) when the number of molecules is not too small. To check the validity of this assumption we have numerically computed the evolution of the distance p(nx , ny , t) − ps (nx , ny ) with N = 50 and we have compared it with an exponential decaying ∝exp ( − t/τ ) where τ is defined by Eq. (34). The numerical results are shown in Fig. 8 where we have used the parameter values π AB = 1; π CA = 1.1; π BC = 1; π BA = 1;π AC = 1; π CB = 1 which correspond to a NESS where = .1. We observe, as in the considered case, that the average field approximation gives a very good estimate of the relaxation characteristic time for the CME. Then we have studied the dependence relaxation time τ (see Eq. (34)) as a function of near a D.B. equilibrium using the average Eqs. (33). We have computed the external fields Aex in the NESS case corresponding to the same parameters π AB = 1; π CA = 1.1; π BC = 1; π BA = 1;π AC = 1; π CB = 1 used in the Fig. (8) and we have varied 0 ≤ ≤ 1 according to (32) to modulate the chemical currents. The results are plotted in Fig. 9(a) where we observe that as the relaxation characteristic time decreases with up to a critical value ∗ at which a bifurcation phenomenon occurs as shown by Fig. 9(a) where the -dependence of the real and imaginary parts of the eigenvalues α k is plotted. After this critical value ( ∗ 0.6799) the relaxation time is not affected by the chemical currents.
FIG. 9. (a) τ (), (b) Re(α()), and (c) Im(α()). The calculation is performed using N = 50 and π AB = 1; π CA = 1.1; π BC = 1; π BA = 1; π AC = 1; π CB = 1.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.189.205.30 On: Wed, 10 Dec 2014 20:35:55
065102-8
de Oliveira et al.
(a)
J. Chem. Phys. 141, 065102 (2014)
(b)
(c)
FIG. 10. (a) Before the critical threshold there are two real eigenvalues, the equilibrium is a stable fixed point, the decaying time τ is dominated by the smaller eigenvalue. (b) Approaching the critical threshold, at which we have the bifurcation of the eigenvalues in the complex plane, the decaying time gets closer to an intermediate value between those associated to the initial eigenvalues. (c) After the critical threshold the fixed point becomes a stable focus (the eigenvalues are complex conjugated), the decaying time is an intermediate value between those associated to the initial eigenvalues, hence smaller than the decaying time before the bifurcation.
From Eq. (32) one can easily prove that the parameter is directly proportional to the chemical currents, and our analysis suggests as the currents can reduce the relaxation time. The underlying bifurcation mechanism in the eigenvalue space is illustrated in Fig. 10 where we plot a sketch of the chemical near the critical point.
VII. CONCLUSIONS
In this paper we have studied the dynamical role of chemical fluxes that characterize the NESS of a unimolecular chemical reaction cycle. Using the correspondence between the CME and a discrete Kolmogorov forward equation we were able to show that the chemical fluxes are linearly proportional to a non-conservative “external vector field,” and that the effect of this field on the system is directly related to the entropy production rate in the NESS. As a consequence, by modulating the external field we can change the chemical fluxes without affecting the stationary probability distribution of the chemical species. In such a way it is possible to study the effect of the fluxes on the relaxation characteristic time of the CME in the case of NESS. We have performed explicit calculations on a linear CME for which it is possible to compute explicitly the stationary probability distribution, the chemical fluxes and the external non-linear field. Our main result is that the stationary fluxes reduce the characteristic relaxation time in comparison to the detailed balance condition. This effect is explained by a bifurcation phenomenon for the eigenvalues of the linearized dynamics around a local maximum of the probability distribution. Even if we have explicitly considered a simple system, we expect that this is a generic result generalizable to a non-linear CME that can be linearized around local maxima. From a biochemical and biological point of view, our study can be relevant for a series of reasons: First, all the biological systems are open systems, hence they exchange materials with their environment, and they consume chemical energy, usually in the form of adenosine triphosphate (ATP) hydrolysis. Second, the phosphorylation (and dephosphorylation) reactions are widely spread into biological systems; the role of these reactions ranges from “energization” of substrates such as in the first stage of sugar catabolism, to post-translational protein modification (the phosphoric group can bee seen as a “tag” that will
drive proteins to different pathways such as cellular compartmentalization or degradation), to enzymatic switching (one enzyme can be turned on or off by reversible phosphorylation associated to conformational changes). Third, cyclic protein phosphorylation and dephosphorylation are relevant for crucial biological processes such as synaptic plasticity, cell fate determination, rate of transcription, and many other. The protein phosphorylation state is central for the cellular signaling processes and its role is under extensive physico-chemical studies both experimental and theoretical. Despite the unquestionable relevance of reversible phosphorylation, we want to remark that the chemical cycle studied here can be applied also to other proteins post-translational modifications such as acetylation, methylation, glycosylation, and many other. All these reactions are based on the reversible addition of a functional groups to proteins, and these modifications lead to an extension of the proteins biological role. Hence, a deeper characterization of biochemical signaling cycles, that are the molecular basis for complex cellular behaviors implemented as a switch between states, is crucial for a better comprehension of biochemical and biological systems’ thermodynamical properties. We can refer to this property, as a “plastic” or “adaptable” behaviour of biological systems, where the term “plastic” means the capability to reach a given molecular state as a function of the environmental conditions. Our finding, that the time necessary to reach the stable distribution depends from the “external field” hence, from the interaction with the environment, without altering the shape of such distribution, can help to better understand how biological systems respond to external perturbations and how they adapt and use their “thermodynamic openness.”
ACKNOWLEDGMENTS
E.G., A.B., and G.C. acknowledge support by the Italian Ministry of Education and Research through the Flagship (PB05) InterOmics, EU projects FibeBiotics (289517), Mission-T2D (600803) and MIMOmics (305280). L.R.O. would like to acknowledge the MONESIA: Erasmus Mundus External Cooperation Window Lot 17.
APPENDIX A: DISCRETE KOLMOGOROV FORWARD EQUATION AND RELATED PROPERTIES
CME (1) can be interpreted as a discrete continuity equation in the hyperplane {|n| = N} and the difference operators Dk± consider all the possible elementary change n of the dynamics. The operator Dk± has the following properties: 1. We can write the relation, − + − Ek+ (1 − Ek−1 Ek− ) = −Dk+ Ek−1 Ek+ . (A1) Dk− = Ek−1
2. We can write a discretized Laplacian, + − −Dk+ Dk− = −(Ek−1 Ek− − 1)(Ek−1 Ek+ − 1) − + = Ek−1 Ek+ − 2 + Ek−1 Ek− .
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.189.205.30 On: Wed, 10 Dec 2014 20:35:55
065102-9
de Oliveira et al.
J. Chem. Phys. 141, 065102 (2014)
3. The operator Dk± has similar properties to a derivative and in particular the product rules read Dk+ f (n)g(n)
=
+ Ek−1 Ek− f (n)Dk+ g(n)
+
g(n)Dk+ f (n),
vanish at Dk , we have n˙ k =
(A2) − − + − − Dk f (n)g(n) = Ek−1 Ek f (n)Dk g(n) + g(n)Dk f (n),
=− −
4. The commutativity property holds
=− 5. We can extend the definition of the operator subset ⊆{|n| = N} according to Dk+
=
+ {Ek−1 Ek− \}
∪
+ {\Ek−1 Ek− },
− − Dk− = {Ek−1 Ek+ \} ∪ {\Ek−1 Ek+ },
and we have the relation Dk± f (n) = f (n)
∀ k.
(A4)
Using the definition (3), CME (1) can be written in the form of discrete continuity equation in the hyperplane {|n| = N}, Dk+ Jk (n, t),
(A7)
k=1
where Jk are the chemical fluxes and are defined according to Jk (n, t) = −πk−1,k (n)nk−1 pn (t) − + Ek−1 Ek+ πk,k−1 (n)nk pn (t)
m
Dk+ (Ak (n)pn (t) + Dk− Bk (n)pn (t))
(A8)
(A9)
k=1
and the fluxes Jk (n, t) can be written in the form Jk (n, t) = Ak (n)pn (t) +
Dk− Bk (n)pn (t),
+ Dk+1 (nk − 1)Jk+1 (n) +
Ak (n)p(n, t) +
Jk+1 (n)
Ak+1 (n)p(n, t). (A11)
Dk− Bk (n)ps (n) = −Ak (n)ps (n).
Dk−
(A12)
− Ek+ ps (n) Ek−1 Dk− ps (n) ln p (n) = ln = ln 1 + ps (n) ps (n) s
from the D.B. equilibrium (A12) and the property (A2), we get the equation Ak (n) + Dk− Bk (n) − s Dk ln p (n) = ln 1 − (A13) − Ek−1 Ek+ Bk (n) that allows to compute the distribution ps (n) in a recursive way using any path connecting a fixed point n0 with a generic point n in the surface |n| = N. Then we define an internal interaction energy V (n), Ak (n) + Dk− Bk (n) − (A14) Dk V (n) = − ln 1 − − Ek−1 Ek+ Bk (n)
E(t) =
Eq. (A7) has the form of a discrete Kolmogorov forward equation, p˙ n (t) = −
Jk (n)
and an internal energy as
= (−πk−1,k (n)nk−1 + πk,k−1 (n)nk )pn (t) + Dk− πk,k−1 (n)nk pn (t).
By using the relation
p˙ n (t) = −
This remark allows to apply the average field approximation in a local way when we have multimodal distribution. To characterize the properties of the stationary solution of the CME we first analyze the D.B. case. Using the definition (A10), the D.B. equilibrium can be written in the form
(A5)
As a consequence if f(n) vanishes at the boundary points we have Dk± f (n) = 0. (A6)
m
nk Dj+ Jj
Dk+ (nk + 1)Jk (n) −
to any
Dk±
j =1
m j =1
(A3) Dk±
m
for any real function f(n) and g(n) with n ∈ N m . Dk± Dh± = Dh± Dk± .
˙ t) = − nk p(n,
V (n)p(n, t)
so that the equilibrium distribution has the MaxwellBoltzmann form (9). In a NESS, the stationary chemical currents Jks (n) are different from zero and this can be associated with the existence of an external non-conservative field interacting with the system, that satisfies the relation Aex k (n) =
(A10)
where Ak (n) = −πk−1,k (n)nk−1 + πk,k−1 (n)nk , Bk (n) = πk,k−1 (n)nk . Eq. (A6) is useful to compute the average dynamics for CME (4): if we choose such that the chemical fluxes Jk (n)
(A15)
|n|=N
Jks (n) , ps (n)
(A16)
i.e., the chemical fluxes are proportional to an external non conservative field. Indeed one can also define the internal field ex Ain k (n) = Ak (n) − Ak (n) that satisfies the conditions s Dk− Bk (n)p s (n) = Ain k (n)p (n),
which follows from the definition (A10) and can be associated to the internal interactions. In such a case the stationary
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.189.205.30 On: Wed, 10 Dec 2014 20:35:55
de Oliveira et al.
065102-10
J. Chem. Phys. 141, 065102 (2014)
solution satisfies
Dk− ln ps (n) = −Dk− V in (n) = ln 1 −
− Ain k (n) + Dk Bk (n) − Ek−1 Ek+ Bk (n)
(A17)
using (A2) we write − Ek+ pn (t))Dk− Bk (n) + Bk (n)Dk− pn (t), Dk− Bk (n)pn (t) = (Ek−1 (A25) putting (A25) in (A24),
and corresponds to the Maxwell-Boltzmann distribution associated to the internal energy V in (n). The stationary condition for the modulated drift (defined in Eq. (14)) reads m
Dk+ Aˆ k (n)p s (n) + Dk− B)k(n)ps (n)
=
+ Bk (n)Dk− pn (t).
Dk+ Ak (n)p s (n) + Dk− B)k(n)pˆ s (n) m
ˆ s (n) = 0 Dk+ Aex k (n)p
(A18)
−
k=1
and it is satisfied if we set pˆ s (n) = ps (n) since the first term vanishes as a consequence of Eq. (8), whereas the second term becomes m m s Dk+ Aex (n)p (n) = λ Dk+ Jks (n) = 0 λ k k=1
due to the definition of vector reads
S˙in = −
The new stationary current
|n|=N
S˙in = −
p˙ n (t) −
|n|=N
=−
m
=
m
+ Ek−1 Ek− Jk (n, t)
(A27)
.
Dk+ pn pn
+ Ek−1 Ek− Jk (n, t)
× =
+ − Ek−1 Ek (Jk (n) − Ak (n)pn (t)) + pn (t)Ek−1 Ek− Bk (n)
m k=1 |n|=N
taking its time derivative, we obtain S˙in = − p˙ n (t)(1 + ln pn (t)),
+ Ek−1 Ek− Bk (n)
k=1 |n|=N
The fluxes Jk can be related to the entropy production, Sin = − pn (t) ln pn (t) (A19)
+ Ek− Dk− Bk (n) Ek−1
k=1 |n|=N
s s Jˆks (n) = Jks (n) − λAex k (n)p (n) = (1 − λ)Jk (n).
|n|=N
we use the
Therefore, we can write the entropy production (A21) as
k=1
Aex k (n).
(A26)
E + E − (J (n) − A (n)pn (t)) Dk+ pn (t) = − k−1 k k + − k pn (t) pn (t)Ek−1 Ek Bk (n)
k=1
−λ
Dk+ pn (t) pn (t)
To find an explicit expression for relation (A1)
k=1 m
− Jk (n) − Ak (n)pn (t) = (Ek−1 Ek+ pn (t))Dk− Bk (n)
+ Ek−1 Ek− Bk (n)
+ (Ek−1 Ek− Jk (n))2
+ pn (t)Ek−1 Ek− Bk (n)
m
−
−
+ Ek−1 Ek− Dk− Bk (n)
+ Ek−1 Ek− Jk (n, t)
k=1 |n|=N
p˙ n (t) ln pn (t)
(A20)
|n|=N
×
+ + Ek−1 Ek− Ak (n)pn (t) + pn (t)Ek−1 Ek− Dk− Bk (n)
+ pn (t)Ek−1 Ek− Bk (n)
(A28)
p˙ n (t) ln pn (t),
|n|=N
and shifting the sum indexes, we get the expression
replacing p˙ n (t) by CME (A7) we have S˙in =
m
Dk+ Jk (n, t) ln pn (t),
(A21)
k=1 |n|=N
m
+ Ek−1 Ek− Jk (n, t)Dk+ ln pn (t).
Dk+
1, we can approximate ln pn (t) as + − E E p D+p D+p Dk+ ln pn = ln k−1 k n = ln 1 + k n k n . pn pn pn (A23) To determine Dk+ ln pn (t) we consider the definition of currents (A10), If
Jk (n) − Ak (n)pn (t) = Dk− Bk (n)pn (t),
×
Jk2 (n) − Jk (n, t) − Bk (n)Ek−1 Ek+ pn (t)
− Ek+ pn (t) Ak (n)pn (t) + Dk− Bk (n)Ek−1 − Bk (n)Ek−1 Ek+ pn (t)
(A22)
k=1 |n|=N
Dk+ pn
m k=1 |n|=N
where we can use the property (A2) and get the explicit form S˙in = −
S˙in =
(A24)
. (A29)
− If N is not too small we can approximate Ek−1 Ek+ pn (t) pn (t) and we get the relation
S˙in =
m k=1 |n|=N
Jk2 (n) Bk (n)pn (t)
Ak (n) + Dk− Bk (n) . Jk (n, t) − Bk (n) k=1 |n|=N m
(A30)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.189.205.30 On: Wed, 10 Dec 2014 20:35:55
de Oliveira et al.
065102-11
J. Chem. Phys. 141, 065102 (2014) 1 H.
APPENDIX B: PROPERTIES OF THE CHEMICAL CURRENTS FOR THE LINEAR CME
We have Dx Jxs (nx , ny )ps (nx , ny ) + Dy Jys (nx , ny )ps (nx , ny ) = 0 (B1) and the divergence of non-equilibrium fluxes are orthogonal to the gradient of probability if s s Dx ps (nx , ny ) + Jy+1 Dy ps (nx , ny ) = 0. Jx+1
(B2)
So we should compute how the operators Dx and Dy act on the stationary solution ps (nx , ny ), (N − nx − ny )λA λ−1 C s − 1 ps (nx , ny ) Dx p (nx , ny ) = nx + 1 (B3) and (N − nx − ny )λB λ−1 C s − 1 ps (nx , ny ). Dy p (nx , ny ) = ny + 1 (B4) We finally can determine that (B2) is verified, s s Dx ps (nx , ny ) + Jy+1 Dy ps (nx , ny ) Jx+1 (N − nx − ny )λA λ−1 C = −1 nx + 1 πAB nx ny λ−1 πBA nx ny λ−1 A λC B λC − × (N − nx − ny ) (N − nx − ny ) −1 + πCA nx λA λC − πAC nx ps (nx , ny )
+ ×
(N − nx − ny )λB λ−1 C ny + 1 πAB nx ny λ−1 B λC (N − nx − ny )
−
−1
πBA nx ny λ−1 A λC
(N − nx − ny ) s + πCB ny λ−1 λ − π n BC y p (nx , ny ) = 0. B C
(B5)
V. Westerhoff and B. O. Palsson, Nat. Biotechnol. 22, 1249 (2004). 2 H. Qian and T. C. Reluga, Phys. Rev. Lett. 94, 028101 (2005). 3 I. R. Epstein and K. Showalter, J. Chem. Phys. 100, 13132 (1996). 4 J. D. Murray, Mathematical Biology (Springer, 2002), Vol. 2. 5 J. Elf and M. Ehrenberg, Syst. Biol. 1, 230 (2004). 6 E. Giampieri, D. Remondini, L. de Oliveira, G. Castellani, and P. Lió, Mol. BioSyst. 7, 2796 (2011). 7 E. L. Bienenstock, L. N. Cooper, and P. W. Munro, J. Neurosci. 10, 79 (1995). 8 N. Intrator and L. N. Cooper, Neural Networks 5, 3 (1992). 9 G. Castellani, A. Bazzani, and L. Cooper, Proc. Natl. Acad. Sci. U.S.A. 106, 14091 (2009). 10 A. Bazzani, G. C. Castellani, E. Giampieri, D. Remondini, and L. N. Cooper, J. Chem. Phys. 136, 235102 (2012). 11 L. Xu, H. Shi, H. Feng, and J. Wang, J. Chem. Phys. 136, 165102 (2012). 12 X. Wang and X. Wang, Nucleic Acids Res. 34, 1646 (2006). 13 H. Ge and H. Qian, Phys. Rev. E 81, 051133 (2010). 14 H. Qian, Annu. Rev. Phys. Chem. 58, 113 (2007). 15 H. Qian, J. Phys. Chem. B 110, 15063 (2006). 16 H. Qian and L. M. Bishop, Int. J. Mol. Sci. 11, 3472 (2010). 17 H. Qian, Phys. Rev. E 65, 016102 (2002). 18 H. Qian and E. L. Elson, Biophys. Chem. 101, 565 (2002). 19 U. Seifert, Rep. Prog. Phys. 75, 126001 (2012). 20 R. Zia and B. J. Schmittmann, J. Stat. Mech. 2007, P07012. 21 R. Zia and B. J. Schmittmann, J. Phys. A: Math. Gen. 39, L407 (2006). 22 N. van Kampen, Stochastic Processes in Physics and Chemistry (NorthHolland Personal Library, 2007). 23 C. Gardiner, Handbook of Stochastic Methods: For Physics, Chemistry & The Natural Sciences (Springer, Berlin, 2004). 24 J. Schnakenberg, Rev. Mod. Phys. 48, 571 (1976). 25 U. Seifert, Phys. Rev. Lett. 95, 040602 (2005). 26 T. Schmiedl and U. Seifert, J. Chem. Phys. 126, 044101 (2007). 27 A. De Martino, D. De Martino, R. Mulet, and G. Uguzzoni, PLoS One 7, e39849 (2012). 28 M. Polettini, Phys. Rev. E 84, 051117 (2011). 29 B. R. Irvin and J. Ross, J. Chem. Phys. 89, 1064 (1988). 30 T. L. Hill, J. Theor. Biol. 10, 442 (1965). 31 T. L. Hill, J. Chem. Phys. 54, 34 (1971). 32 T. L. Hill, Proc. Natl. Acad. Sci. U.S.A. 85, 5345 (1988). 33 S. R. de Groot and P. Mazur, Non-Equilibrium Thermodynamics (Dover Books on Physics, 1984). 34 L. Jiu-Li, C. Van den Broeck, and G. Nicolis, Z. Phys. B: Condens. Matter 56, 165 (1984). 35 P. Gaspard, C. R. Phys. 8, 598 (2007).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.189.205.30 On: Wed, 10 Dec 2014 20:35:55
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.189.205.30 On: Wed, 10 Dec 2014 20:35:55
THE JOURNAL OF CHEMICAL PHYSICS 141, 065102 (2014)
The role of non-equilibrium fluxes in the relaxation processes of the linear chemical master equation Luciana Renata de Oliveira, Armando Bazzani, Enrico Giampieri, and Gastone C. Castellania) Physics and Astronomy Department, Bologna University and INFN Sezione di Bologna, Italy
(Received 28 March 2014; accepted 16 July 2014; published online 11 August 2014) We propose a non-equilibrium thermodynamical description in terms of the Chemical Master Equation (CME) to characterize the dynamics of a chemical cycle chain reaction among m different species. These systems can be closed or open for energy and molecules exchange with the environment, which determines how they relax to the stationary state. Closed systems reach an equilibrium state (characterized by the detailed balance condition (D.B.)), while open systems will reach a non-equilibrium steady state (NESS). The principal difference between D.B. and NESS is due to the presence of chemical fluxes. In the D.B. condition the fluxes are absent while for the NESS case, the chemical fluxes are necessary for the state maintaining. All the biological systems are characterized by their “far from equilibrium behavior,” hence the NESS is a good candidate for a realistic description of the dynamical and thermodynamical properties of living organisms. In this work we consider a CME written in terms of a discrete Kolmogorov forward equation, which lead us to write explicitly the non-equilibrium chemical fluxes. For systems in NESS, we show that there is a nonconservative “external vector field” whose is linearly proportional to the chemical fluxes. We also demonstrate that the modulation of these external fields does not change their stationary distributions, which ensure us to study the same system and outline the differences in the system’s behavior when it switches from the D.B. regime to NESS. We were interested to see how the non-equilibrium fluxes influence the relaxation process during the reaching of the stationary distribution. By performing analytical and numerical analysis, our central result is that the presence of the non-equilibrium chemical fluxes reduces the characteristic relaxation time with respect to the D.B. condition. Within a biochemical and biological perspective, this result can be related to the “plasticity property” of biological systems and to their capabilities to switch from one state to another as is observed during synaptic plasticity, cell fate determination, and differentiation. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4891515] I. INTRODUCTION
The dynamics of a chemical cycle chain reaction among m different species is relevant for a large number of biochemical systems, because they are fundamental for the kinetic description of biochemical networks; both in the case of small and large network cycles. These reaction cycles are the chemical basis of several cellular signaling processes such as signal transduction,1, 2 biological morphogenesis,3, 4 and cell fate determination.5, 6 These cycles are also relevant in the induction of synaptic plasticity, where a small network of phosphorylation/dephosphorylation reactions is taking place7–9 at the post-synaptic receptors level. How to describe these systems is a very important topic: their dynamics can be well approximated with nonlinear differential equations,9, 10 but in the case of small number of molecules we have to take into account for fluctuations and we need to apply stochastic methods such as the Chemical Master Equation (CME). The reaction cycles are the simplest non-trivial example of dynamical networks for which it is possible to develop an analytical approach. The characterization of the stationary distribution of a) Electronic mail: [email protected]
0021-9606/2014/141(6)/065102/11/$30.00
a CME of a reaction network is still a challenge, because one should deal with the fact that these systems can be open systems (here we understand open, as the possibility of exchanging molecules and energy with the surroundings).11–21 After a sufficient long time, isolated systems reach a thermodynamics equilibrium state (detailed balance (D.B.) equilibrium), where all the chemical currents are null. While an open system (if the exchange of energy with its surroundings is sustained) approaches a non-equilibrium steady state (NESS),13, 19–28 the probability chemical fluxes are not null leading to the breaking down of the D.B. When the D.B. holds, the equilibrium distribution can be written as a Maxwell-Boltzmann distribution, which means there is a unique equilibrium constant (the equilibrium temperature) for every chemical reaction in a system, regardless of how complex the system is. However, in a NESS, the stationary chemical currents are different from zero and this can be associated with the existence of an external non-conservative field interacting with the system11, 19, 25, 26 and with a positive entropy production due to the currents.29 As a matter of fact most of the chemical cycles relax toward NESS so that the problem of understanding the dynamical and thermodynamical properties of NESS could suggest if the energy cost of NESS with re-
141, 065102-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.189.205.30 On: Wed, 10 Dec 2014 20:35:55
065102-2
de Oliveira et al.
spect to the D.B. equilibrium is compensated by other advantages. In this paper we discuss a general approach to study the dynamics of a chemical chain reaction among m different species, based on the results of previous works,11, 13, 19–21, 24, 25 that allows to compute a drift and a diffusion coefficient and to introduce a non-equilibrium thermodynamical description for the stationary states of CME. The determination of stationary fluxes allows us to completely describe systems in equilibrium (D.B.) as well as in NESS.11, 15, 19, 21, 28 We discuss the role of chemical currents in the transient states and, in particular, in the relaxation process towards the stationary distribution. Indeed the relaxation times of biochemical reactions could be related to the plasticity properties of biological systems, where plasticity means the capability to reach the stationary distribution or to jump from one stationary state to another. The paper is organized as follows: Sec. II is devoted to the biochemical relevance of the unimolecular reaction cycle; in Sec. III we described the model and we discuss the meanfield approximation; in Sec. IV the thermodynamical formalism is introduced to characterize the stationary solutions; in Sec. V we explicitly study the effect of chemical fluxes on the relaxation time toward the stationary states in Sec. VI some numerical results are presented to check the analytic and the thermodynamically assumptions.
II. BIOCHEMICAL RELEVANCE OF LINEAR CYCLES
The more general, and probably realistic, form of the Michaelis-Menten reaction is obtained by allowing a certain degree of reversibility in the product formation step. This consideration leads to a model that was originally developed by Haldane in 1930. This mechanism can be used to model enzymatic reactions, and the kinetic scheme is as follows: This chain of reactions can be represented as a cycle between enzymatic states (E, ES, EP). In Fig. 1 the free enzyme E, after binding the substrate S, forms the first reaction intermediate, the enzyme-substrate complex (ES). This complex, in turn, is transformed in the second reaction intermediate, the enzyme-product complex (EP), where the substrate has been transformed in the product P. As final step, the product is released and the enzyme turns back into the free form. The ES complex can dissociate back in the free enzyme E and the unbounded substrate S, and, in the same way, the association reaction between the free enzyme E and the product P can form the EP complex. If the abundance of the substrate and the product are in chemical balance with the environment, the transition rates of the cycle can be approximated by pseudofirst order rate constants (see Fig. 2). From a general point of view, the possible states that the enzymatic complex can assume would be more numer-
FIG. 1. Reversible Michaelis-Menten reaction scheme. The transformation between substrate and product is composed of a chain of reversible steps.
J. Chem. Phys. 141, 065102 (2014)
FIG. 2. Reversible Michaelis-Menten reduced chain reaction. Under the right conditions the transformation can be approximated as a cycle of linear reversible transformations.
ous than two, but are usually still representable as a cyclic enzymatic transformation. A detailed analysis of the thermodynamical properties of these kind of systems has been done in Ref. 14. This kind of cyclic reactions can be founded in several molecular mechanisms underlying important biological processes such as cell differentiation, cell fate determination, and synaptic plasticity. The bidirectional synaptic plasticity (the transition between long-term synaptic potentiation (LTP) and long-term synaptic depression (LTD) ) is generally assumed as the basis of learning and memory. The molecular mechanism is a phosphorylation/dephosphorylation of sites on the postsynaptic glutamate receptor α-amino-3-hydroxy-5-methyl4-isoxazolepropionic acid receptor (AMPA), on the subunit GluR1. This subunit can be phosphorylated by specific protein-kinase enzymes (such as CaMKII, protein kinase C and A). The induction of LTP specifically increases phosphorylation of one site, while the induction of LTD is accompanied by a decrease in the phosphorylation of the other site. Therefore, knowledge of the phosphorylation state of GluR1 may be a strong predictor of synaptic plasticity. Moreover, the phosphorylation state of AMPA receptors in the same subunit is affecting their cytoplasmatic trafficking, specially their insertion/removal from the postsynaptic membrane (the double phosphorylated AMPA receptor is more expressed in the cellular membrane, while the non-phosphorylated is mainly in the cytoplasm). This phosphorylation/dephosphorylation cycle is represented schematically in Fig. 3. The fraction of AMPARs containing GluR1 phosphorylated at one site is denoted by A∗ , whereas A∗ denotes the fraction of AMPARs containing GluR1s phosphorylated on another, non equivalent, site, and A and A∗∗ denote AMPARs phosphorylation at neither and both sites, respectively. EP1 and EP2 are used to denote the protein phosphatases that dephosphorylate these two sites, and EK1 and EK2 are used to denote the protein kinases that phosphorylate these sites. This molecular cycle, if we assume low concentration of substrates (the number of AMPARs, that in a single spine are on the order of hundred molecules) with respect to the saturation limit of the enzymes (the Michaelis-Menten constant Km ), is linear as the cycle in Fig. 2 and the same considerations can be used for both systems. Such a system is, by definition, always far from thermodynamical equilibrium, and the organism has to continuously
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.189.205.30 On: Wed, 10 Dec 2014 20:35:55
065102-3
de Oliveira et al.
J. Chem. Phys. 141, 065102 (2014)
where n ∈ N m and Ek± are the Van Kampen step operators defined according to Ek± f (n) = f (n1 , n2 , .., nk ± 1, .., nm )
n ∈ N m,
for any function f(n). If the reaction chain is a cycle, we impose periodic boundary conditions in the sum m + 1 ∼ 1 and the total number of particles is constant |n| =
m
nk = N.
(2)
k=1
FIG. 3. Representation of phosphorylation/dephosphorylation cycle of AMPA receptors. Under the right conditions the transformation between the receptors’ phosphorylation state can be approximated by a cycle of four linear reversible transformations driven by two couples of kinases and phosphatases.
In a general case the transition probability π k − 1, k may depend on the distribution n of particles into the different chemical species. B. The drift and diffusion coefficients for the discrete case
spend energy to maintain the system in this state and to ensure the capability to switch between different phosphorylation states.
If one introduces the difference operators Dk± defined as + Dk+ ≡ Ek−1 Ek− − 1
− Dk− ≡ Ek−1 Ek+ − 1,
(3)
Eq. (1) can be written as a discrete Kolmogorov forward equation11, 15, 19, 20 (see Appendix A): III. NON-EQUILIBRIUM FLUXES AND STATIONARY STATES FOR THE CME
˙ t) = − p(n,
Let us consider the dynamics of a chemical chain reaction among m different species (see Fig. 4). Under the assumption that each particle is independent the transition rate from the specie k − 1 to the specie k is given by π k − 1, k (n)nk − 1 where nk − 1 denotes the number of particles of the k − 1 specie and π k − 1, k is the probability that a single particle performs a transition in a time unit. Then we write the CME for the evolution of the probability distribution p(n, t) according to22 ˙ t) = p(n,
+ [(Ek−1 Ek−
Dk+ Jk
k=1
A. Definition of the model
m
− 1)πk−1,k (n)nk−1 p(n, t)
m
Dk+ [Ak (n)pn (t) + Dk− Bk (n)pn (t)]. (4)
k=1
In the continuous limit, Eq. (4) becomes a Fokker-Planck equation and the difference operators tend to the partial differential operators on the tangent space of the constraint (2). The vector Jk denotes the chemical currents and we have introduced the drift field Ak (n) = −πk−1,k (n)nk−1 + πk,k−1 (n)nk and the diffusion coefficient
k − Ek+ − 1)πk,k−1 (n)nk p(n, t)], + (Ek−1
=−
(1)
Bk (n) = πk,k−1 (n)nk . The drift field is directly correlated to the average dynamics: let {n ∈ N m / |n| = N } such that the chemical fluxes (or the distribution p(n, t)) almost vanish at its boundary then we have the average equations (see Eq. (A11) in Appendix A) n˙ k = Ak+1 (n) − Ak (n)
k = 1, .., m
(5)
with an error of the order of the chemical flux at the boundary of . We say that the mean-field approximation can be applied to Eq. (5) on , if n˙ k Ak+1 (n ) − Ak (n )
j = 1, .., m.
(6)
This approximation is correctly applied when the fluctuations are small with respect to the average values, and it is exact when the transition probabilities π k − 1, k are constant (i.e., the particles are independent). C. Solution of the linear CME FIG. 4. Chemical chain reaction among m different species. Each transition is regulated by its own transition rate, that, in general, depends on the number of particles in that state.
In the last case we get a linear CME whose timedependent solution can be explicitly computed in the form
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.189.205.30 On: Wed, 10 Dec 2014 20:35:55
065102-4
de Oliveira et al.
J. Chem. Phys. 141, 065102 (2014)
of a multinomial distribution,30–32 assuming that the system starts from a known state (p(n∗ , t = 0) = 1 for some n∗ ), p(n, t) = N !
n m λ k (t)
|n| = N,
k
k=1
nk !
(7)
Appendix A we show that the time derivative of the Gibbs’ entropy for the system can be written in the form33–35 S˙in =
k=1 |n|=N
|λ| =
λk = 1
λk > 0.
k=1
Letting t → ∞ we get the stationary distribution ps (n) with limt→∞ λk (t) = λ∗k where λ∗k satisfies the condition Ak+1 (N λ∗ ) = Ak (N λ∗ ). n∗k
When the average field approximation holds the values Nλ∗k , this corresponds to the local maximum of the stationary distribution ps (n) if the diffusion coefficient can be considered constant. Therefore, in the linear case the critical value of the stationary distribution corresponds to the stable fixed point of the average systems (6). IV. THERMODYNAMICAL PROPERTIES OF CME
A thermodynamical approach has been proposed to characterize the properties of the stationary solutions of CME.24 In particular one distinguishes the equilibrium states at which the chemical fluxes Jks (n) cancel out (D.B. condition for equilibrium states), and the NESS11, 13, 19, 20, 22–24, 26–28 where the weaker condition holds (see Eq. (4)), m
Dk+ Jks (n) = 0.
(8)
k=1
We will explicitly recall some properties of the stationary solution for Eq. (4). When the D.B. holds, the equilibrium distribution can be written as a Maxwell-Boltzmann distribution
Jk2 (n) Bk (n)pn (t)
Ak (n) + Dk− Bk (n) , (12) Jk (n, t) − Bk (n) k=1 |n|=N m
where the quantities λk are the solutions of the linear system (6) with the constraint m
m
when the stationary distribution satisfies Dk+ pn (t) 1 (continuous approximation) and it is possible to linearize the discrete derivative Dk+ ln pn (t) Dk+ pn (t)/pn (t). In the rhs of Eq. (12) one recognizes the total entropy production S˙ and the environment entropy production S˙en according to S˙ =
m k=1 |n|=N
S˙en
Jk2 (n) , Bk (n)pn (t)
(13) Ak (n) + Dk− Bk (n) . =− Jk (n, t) Bk (n) k=1 |n|=N m
In this way we get the mathematical formulation for the thermodynamic variables in terms of the master equation. When the system presents detailed balance condition, for t s s = S˙ s = S˙en = → ∞ it reaches an equilibrium state with S˙in s s 0. In contrast systems in NESS, S˙ and S˙en are equal but not s = 0. null, to ensure that S˙in Therefore, in a NESS it is straightforward to observe that the total entropy production is always positive, so that to maintain the stationary distribution the environment is exchanging energy with the system through the work done by the external forces, which is dissipated. Conversely in a D.B. equilibrium the total entropy production is zero and there is no dissipated work. By using the decomposition (10) it is possible to modulate the external field by changing the drift term according to Aˆ k (n) = Ak (n) − ηAex k (n),
p (n) ∝ exp(−V (n)), s
(9)
where V (n) is an interaction microscopic energy. The expectation value of V (n) can be identified with the internal energy of the system and it is natural to introduce the Gibbs’ entropy to describe the heat exchange with the environment.24 In a NESS it is possible to split the drift field into an internal and an external vector field,11, 14, 19–21, 25 Ak (n) =
Ain k (n)
+
Aext k (n),
(10)
where Aext k (n) is an external non-conservative field which generates the chemical fluxes Aex k (n) =
Jks (n) ps (n)
(14)
where η is a parameter: η = 0 corresponds to the initial case and η = 1 to the D.B. equilibrium when the external field vanishes. Moreover, one proves (see Appendix A) that the stationary distribution does not depend on η, whereas the new fluxes read Jˆks (n) = (1 − η)Jks (n).
(15)
Therefore, from this point of view it seems to be not convenient for any biochemical system to relax towards a NESS due to its energetic and entropic cost. Then the question is to study the effect of chemical fluxes in the transient states and, in particular, in the relaxation process towards the stationary distribution.
(11)
and Ain k (n) is a conservative vector whose potential satisfies V in (n) = − ln ps (n). The external field Aex k (n) is proportional to the thermodynamical force associated to the chemical fluxes. As it is known19–21, 25 in a NESS, the stationary fluxes Jks (n) are related to the total entropy production. In
V. NON-EQUILIBRIUM FLUXES IN THE LINEAR CME
We study the relaxation times of a linear CME in D.B. equilibrium and NESS to understand the influence of the nonequilibrium fluxes. For seek of simplicity we have chosen a chemical reaction with three states, as represented in Fig. 5
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.189.205.30 On: Wed, 10 Dec 2014 20:35:55
065102-5
de Oliveira et al.
J. Chem. Phys. 141, 065102 (2014)
FIG. 5. Unimolecular chemical reaction cycle. The cycle is composed only of reactions with a transition rate that is linear in the number of particles.
FIG. 6. Stationary distribution for CME (16) computed using (17) for N = 50 as a function of nA and nB .
Eq. (16) reads − ˙ x , ny , t) = πAB (E+ p(n x Ey − 1)nx p(nx , ny )
whose CME (see Eq. (1)) is
+ + πBA (E− x Ey − 1)ny p(nx , ny )
˙ t) p(n, =
− (E+ C EA
− 1)πCA nC p(n, t) +
− × (E+ A EB − × (E+ B EC
+ (E− C EA
+ πCA (E− x − 1)(N − nx − ny )p(nx , ny )
− 1)πAC nA p(n, t)
− 1)πAB nA p(n, t) +
+ (E− A EB
− 1)πBA nB p(n, t)
− 1)πBC nB p(n, t) +
+ (E− B EC
− 1)πCB nC p(n, t),
+ πCB (E− y − 1)(N − nx − ny )p(nx , ny ) + πAC (E+ x − 1)nx p(nx , ny ) + πBC (E+ y − 1)ny p(nx , ny ).
(16) where nA , nB , and nC are the number of molecules of the species A, B, and C, respectively, with the constraint nA + nB + nC = N. The particle distribution is a multinomial distribution (see (7)), n
n
n
λ A (t)λBB (t)λCC (t) , p(n, t) = N! A nA !nB !nC !
(17)
where λA (t), λB (t), λC (t) are the solutions of the linear system
(20)
We write CME (20) in the form of a continuity Eq. (4) without any constraint, ˙ x , ny ) = −Dx Jx (nx , ny ) − Dy Jy (nx , ny ) p(n
(21)
by using the discrete difference operators Dx = E+ x −1
Dy = E+ y − 1.
So we identify the fluxes Jx (nx , ny ) = JA (n) − JB (n) and Jy (nx , ny ) = JB (n) − JC (n) as − Jx (nx , ny ) = −(πAB E− y nx − πBA Ex ny + πAC nx
λ˙ A (t) = −(πAB + πAC )λA (t) + πBA λB (t) + πCA λC (t) λ˙ B (t) = πAB λA (t) − (πBA + πBC )λB (t) + πCB λC (t)
− πCA E− x (N − nx − ny ))p(nx , ny ), (18)
λ˙ C (t) = πAC λA (t) + πBC λB (t) − (πCA + πCB )λC (t)
Substituting the explicit form of the steady state solution pns x ,ny in (22), we determine the stationary non-equilibrium currents Jxs (nx , ny ) and Jys (nx , ny ), πBA πAB = − s N − nx − ny + 1 λsA λB λs + Cs πCA nx − πAC nx ps (nx , ny ) λA (23) λsC nx ny πAB πBA s − s Jy (nx , ny ) = N − nx − ny + 1 λsB λA λs + Cs πCB ny − πBC ny ps (nx , ny ) λB
λsA ∝ πBA πCA + πBC πCA + πBA πCB λsC
(22)
− πCB E− y (N − nx − ny ))p(nx , ny ).
with the constraint λA + λB + λC = 1. The steady state solution ps (n) is computed taking the limit t → ∞ in Eq. (17) and the limit values λs can be explicitly computed as follows:
λsB ∝ πAB πCA + πAB πCB + πAC πCB
− Jy (nx , ny ) = −(πBA E− x ny − πAB Ey nx + πBC ny
Jxs (nx , ny ) (19)
∝ πAC πBA + πAB πBC + πAC πBC .
Fig. 6 represents the stationary distribution of CME (16), for N = 50 as a function of nA and nB : the maximal value corresponds to nA = N λsA and nB = N λsB . For an explicit computation of the chemical fluxes we reduce the system’s dimensionality using nA = nx ; nB = ny and nC = N − nx − ny . Then
λsC nx ny
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.189.205.30 On: Wed, 10 Dec 2014 20:35:55
065102-6
de Oliveira et al.
J. Chem. Phys. 141, 065102 (2014)
According to Eq. (4), we write the non-equilibrium fluxes introducing a drift term and a diffusion coefficient Jx (nx , ny ) = Ax (nx , ny )p(nx , ny ) − Dx Bxx p(nx − 1, ny) − Dy Bxy p(nx , ny − 1),
(24) Jy (nx , ny ) = Ay (nx , ny )p(nx , ny ) − Dx Byx p(nx − 1, ny) − Dy Byy p(nx , ny − 1).
B=
+ πCA (N − nx − ny )
(25)
Ay (nx , ny ) = πAB nx − (πBA + πBC )ny + πCB (N − nx − ny )
πBA ny + πCA (N − nx − ny + 1)
−πAB nx
−πBA ny
πAB nx + πCB (N − nx − ny + 1)
=
(26)
s Aex y (nx , ny )p (nx , ny ).
+Jys (nx , ny + 1)Dy ps (nx , ny ) = 0.
(27)
λsC nx ny
+
Aex y (nx , ny ) =
N − nx − ny + 1 +
To compute the relaxation time towards the stationary distribution we take advantage from the fact that a multinomial distribution is completely determined by the average values. Then we study the average dynamics, n˙ x = −(πAB + πAC )nx + πBA ny + πCA (N − nx − ny ), n˙ y = +πAB nx − (πBA + πBC )ny + πCB (N − nx − ny ).
(28)
The relaxation process towards the limit values is an exponential whose exponents are the eigenvalues of Eq. (28). To study the effects of the chemical fluxes we modulate the external field which produces the fluxes by changing the drift term according to Aˆ x,y (n) = Ax,y (n) − ηAex x,y (n) (see Eq. (14)). In this way we change the non-equilibrium fluxes without modifying the stationary distribution (see Appendix A). From Eq. (15) the new non-equilibrium fluxes read
Jˆy (nx , ny , t) = Jy (nx , ny , t) − ηAex y (nx , ny )p(nx , ny , t).
πAB πBA s − λB λsA
λsC πCB s − πBC ny . λB
(30)
(31)
It is convenient to set η = 1 + , so that = 0 corresponds to the D.B. equilibrium and changing , we drive the system to a NESS condition. From (21) we get the modified CME, ˙ x , ny , t) p(n
N λsA , N λsB
Jˆx (nx , ny , t) = Jx (nx , ny , t) − ηAex x (nx , ny )p(nx , ny , t),
πBA πAB s − λA λsB
λsC πCA s − πAC nx , λA
λsC nx ny
VI. RESULTS
N − nx − ny + 1
In Appendix B we show that the fluxes Jxs (nx , ny ) and Jys (nx , ny ) are orthogonal to the discrete gradient of the stationary distribution ps (nx , ny ), Jxs (nx + 1, ny )Dx ps (nx , ny )
.
where we explicitly have Aex x (nx , ny )
s Jxs (nx , ny ) = Aex x (nx , ny )p (nx , ny ),
=
Ax (nx , ny ) = −(πAB + πAC )nx + πBA ny
and the diffusion matrix B is
Then we define an external field Aex (cf. Eq. (11)) related to the non-equilibrium fluxes as
Jys (nx , ny )
The drift vector A(nx , ny ) is
(29)
= −Dx (Jx (nx , ny , t) − (1 + )Aex x (nx , ny , t)p(nx , ny , t)) − Dy (Jy (nx , ny , t) − (1 + )Aex y (nx , ny , t)p(nx , ny , t)). (32) We remark that the new CME is non-linear. From the explicit knowledge of the stationary chemical currents, we are able to control the NESS condition for systems with = 0, using the s conditions Jks = 0 and S˙ s = S˙en
= 0. To confirm the NESS, we calculated the value of the norm of the non-equilibrium currents (||Js ||) as a function of (see Figure 7(a)) and we s estimate the values of S˙ s and S˙en for 0 ≤ ≤ 1 in the stationary state (see Figure 7(b)). In this way we verify for all s are not null and that values of = 0 that ||Js ||, and S˙ s = S˙en the system is in a NESS state.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.189.205.30 On: Wed, 10 Dec 2014 20:35:55
065102-7
de Oliveira et al.
J. Chem. Phys. 141, 065102 (2014) Ss 300 250 200 150 100 (a)
(b)
50 0.2
0.4
0.6
0.8
1.0
Ε
FIG. 7. (a) Norm of the non-equilibrium currents (||Js ||) as a function of , (b) Entropy production rate in the stationary state (Ss ) as a function of . The calculation is performed using N = 50 and π AB = 1; π CA = 1.1; π BC = 1; π BA = 1; π AC = 1; π CB = 1.
To study the relaxation process of CME (32) when the chemical fluxes change we consider the average equations,
∂Aex x (nx , ny )
n n˙ x = − πAB + πAC + πCA + (1 + ) x
∗ ∂nx nx
∂Aex x (nx , ny )
n , + πBA − πCA − (1 + ) y
∗ ∂ny ny (33)
∂Aex y (nx , ny )
n + n˙ y = πAB − πCB − (1 + ) x
∗ ∂nx nx
∂Aex y (nx , ny )
n − (πBA + πBC + πCB + (1 + ) y
∗ ∂ny ny in the mean field approximation. The evolution of the average values n(t) relaxes towards the stable critical point of Eqs. (33) with a characteristic time τ that depends on the minimum of the real part of eigenvalues α k k = 1, 2 computed linearizing the equation at the critical points, τ=
1 . Min|Re(α1 , α2 )|
(34)
The nonlinear character of the external fields Aex implies that the average Eqs. (33) are not exact, but we expect that the mean field approximation holds in the neighborhood of the
FIG. 8. Comparison between the relaxation time for the distance p(nx , ny , t) − ps (nx , ny ) for the solution of CME (32) and the relaxation time estimated by (34) using the average Eqs. (33). The parameters used in the simulations are p(nx, ny, 0) = 1/N2 ,π AB = 1; π CA = 1.1; π BC = 1; π BA = 1; π AC = 1; π CB = 1, which correspond to a NESS case with = .9, and N = 50 molecules. The black line is the norm p(nx , n)y, t) − ps (nx , ny ) and the gray line gives an interpolation with the exponential decaying cexp ( − t/τ ).
maximum of the stationary distribution function ps (nx , ny ) for the modified CME (32) that is well approximated by the stable critical point of the average Eq. (33) when the number of molecules is not too small. To check the validity of this assumption we have numerically computed the evolution of the distance p(nx , ny , t) − ps (nx , ny ) with N = 50 and we have compared it with an exponential decaying ∝exp ( − t/τ ) where τ is defined by Eq. (34). The numerical results are shown in Fig. 8 where we have used the parameter values π AB = 1; π CA = 1.1; π BC = 1; π BA = 1;π AC = 1; π CB = 1 which correspond to a NESS where = .1. We observe, as in the considered case, that the average field approximation gives a very good estimate of the relaxation characteristic time for the CME. Then we have studied the dependence relaxation time τ (see Eq. (34)) as a function of near a D.B. equilibrium using the average Eqs. (33). We have computed the external fields Aex in the NESS case corresponding to the same parameters π AB = 1; π CA = 1.1; π BC = 1; π BA = 1;π AC = 1; π CB = 1 used in the Fig. (8) and we have varied 0 ≤ ≤ 1 according to (32) to modulate the chemical currents. The results are plotted in Fig. 9(a) where we observe that as the relaxation characteristic time decreases with up to a critical value ∗ at which a bifurcation phenomenon occurs as shown by Fig. 9(a) where the -dependence of the real and imaginary parts of the eigenvalues α k is plotted. After this critical value ( ∗ 0.6799) the relaxation time is not affected by the chemical currents.
FIG. 9. (a) τ (), (b) Re(α()), and (c) Im(α()). The calculation is performed using N = 50 and π AB = 1; π CA = 1.1; π BC = 1; π BA = 1; π AC = 1; π CB = 1.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.189.205.30 On: Wed, 10 Dec 2014 20:35:55
065102-8
de Oliveira et al.
(a)
J. Chem. Phys. 141, 065102 (2014)
(b)
(c)
FIG. 10. (a) Before the critical threshold there are two real eigenvalues, the equilibrium is a stable fixed point, the decaying time τ is dominated by the smaller eigenvalue. (b) Approaching the critical threshold, at which we have the bifurcation of the eigenvalues in the complex plane, the decaying time gets closer to an intermediate value between those associated to the initial eigenvalues. (c) After the critical threshold the fixed point becomes a stable focus (the eigenvalues are complex conjugated), the decaying time is an intermediate value between those associated to the initial eigenvalues, hence smaller than the decaying time before the bifurcation.
From Eq. (32) one can easily prove that the parameter is directly proportional to the chemical currents, and our analysis suggests as the currents can reduce the relaxation time. The underlying bifurcation mechanism in the eigenvalue space is illustrated in Fig. 10 where we plot a sketch of the chemical near the critical point.
VII. CONCLUSIONS
In this paper we have studied the dynamical role of chemical fluxes that characterize the NESS of a unimolecular chemical reaction cycle. Using the correspondence between the CME and a discrete Kolmogorov forward equation we were able to show that the chemical fluxes are linearly proportional to a non-conservative “external vector field,” and that the effect of this field on the system is directly related to the entropy production rate in the NESS. As a consequence, by modulating the external field we can change the chemical fluxes without affecting the stationary probability distribution of the chemical species. In such a way it is possible to study the effect of the fluxes on the relaxation characteristic time of the CME in the case of NESS. We have performed explicit calculations on a linear CME for which it is possible to compute explicitly the stationary probability distribution, the chemical fluxes and the external non-linear field. Our main result is that the stationary fluxes reduce the characteristic relaxation time in comparison to the detailed balance condition. This effect is explained by a bifurcation phenomenon for the eigenvalues of the linearized dynamics around a local maximum of the probability distribution. Even if we have explicitly considered a simple system, we expect that this is a generic result generalizable to a non-linear CME that can be linearized around local maxima. From a biochemical and biological point of view, our study can be relevant for a series of reasons: First, all the biological systems are open systems, hence they exchange materials with their environment, and they consume chemical energy, usually in the form of adenosine triphosphate (ATP) hydrolysis. Second, the phosphorylation (and dephosphorylation) reactions are widely spread into biological systems; the role of these reactions ranges from “energization” of substrates such as in the first stage of sugar catabolism, to post-translational protein modification (the phosphoric group can bee seen as a “tag” that will
drive proteins to different pathways such as cellular compartmentalization or degradation), to enzymatic switching (one enzyme can be turned on or off by reversible phosphorylation associated to conformational changes). Third, cyclic protein phosphorylation and dephosphorylation are relevant for crucial biological processes such as synaptic plasticity, cell fate determination, rate of transcription, and many other. The protein phosphorylation state is central for the cellular signaling processes and its role is under extensive physico-chemical studies both experimental and theoretical. Despite the unquestionable relevance of reversible phosphorylation, we want to remark that the chemical cycle studied here can be applied also to other proteins post-translational modifications such as acetylation, methylation, glycosylation, and many other. All these reactions are based on the reversible addition of a functional groups to proteins, and these modifications lead to an extension of the proteins biological role. Hence, a deeper characterization of biochemical signaling cycles, that are the molecular basis for complex cellular behaviors implemented as a switch between states, is crucial for a better comprehension of biochemical and biological systems’ thermodynamical properties. We can refer to this property, as a “plastic” or “adaptable” behaviour of biological systems, where the term “plastic” means the capability to reach a given molecular state as a function of the environmental conditions. Our finding, that the time necessary to reach the stable distribution depends from the “external field” hence, from the interaction with the environment, without altering the shape of such distribution, can help to better understand how biological systems respond to external perturbations and how they adapt and use their “thermodynamic openness.”
ACKNOWLEDGMENTS
E.G., A.B., and G.C. acknowledge support by the Italian Ministry of Education and Research through the Flagship (PB05) InterOmics, EU projects FibeBiotics (289517), Mission-T2D (600803) and MIMOmics (305280). L.R.O. would like to acknowledge the MONESIA: Erasmus Mundus External Cooperation Window Lot 17.
APPENDIX A: DISCRETE KOLMOGOROV FORWARD EQUATION AND RELATED PROPERTIES
CME (1) can be interpreted as a discrete continuity equation in the hyperplane {|n| = N} and the difference operators Dk± consider all the possible elementary change n of the dynamics. The operator Dk± has the following properties: 1. We can write the relation, − + − Ek+ (1 − Ek−1 Ek− ) = −Dk+ Ek−1 Ek+ . (A1) Dk− = Ek−1
2. We can write a discretized Laplacian, + − −Dk+ Dk− = −(Ek−1 Ek− − 1)(Ek−1 Ek+ − 1) − + = Ek−1 Ek+ − 2 + Ek−1 Ek− .
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.189.205.30 On: Wed, 10 Dec 2014 20:35:55
065102-9
de Oliveira et al.
J. Chem. Phys. 141, 065102 (2014)
3. The operator Dk± has similar properties to a derivative and in particular the product rules read Dk+ f (n)g(n)
=
+ Ek−1 Ek− f (n)Dk+ g(n)
+
g(n)Dk+ f (n),
vanish at Dk , we have n˙ k =
(A2) − − + − − Dk f (n)g(n) = Ek−1 Ek f (n)Dk g(n) + g(n)Dk f (n),
=− −
4. The commutativity property holds
=− 5. We can extend the definition of the operator subset ⊆{|n| = N} according to Dk+
=
+ {Ek−1 Ek− \}
∪
+ {\Ek−1 Ek− },
− − Dk− = {Ek−1 Ek+ \} ∪ {\Ek−1 Ek+ },
and we have the relation Dk± f (n) = f (n)
∀ k.
(A4)
Using the definition (3), CME (1) can be written in the form of discrete continuity equation in the hyperplane {|n| = N}, Dk+ Jk (n, t),
(A7)
k=1
where Jk are the chemical fluxes and are defined according to Jk (n, t) = −πk−1,k (n)nk−1 pn (t) − + Ek−1 Ek+ πk,k−1 (n)nk pn (t)
m
Dk+ (Ak (n)pn (t) + Dk− Bk (n)pn (t))
(A8)
(A9)
k=1
and the fluxes Jk (n, t) can be written in the form Jk (n, t) = Ak (n)pn (t) +
Dk− Bk (n)pn (t),
+ Dk+1 (nk − 1)Jk+1 (n) +
Ak (n)p(n, t) +
Jk+1 (n)
Ak+1 (n)p(n, t). (A11)
Dk− Bk (n)ps (n) = −Ak (n)ps (n).
Dk−
(A12)
− Ek+ ps (n) Ek−1 Dk− ps (n) ln p (n) = ln = ln 1 + ps (n) ps (n) s
from the D.B. equilibrium (A12) and the property (A2), we get the equation Ak (n) + Dk− Bk (n) − s Dk ln p (n) = ln 1 − (A13) − Ek−1 Ek+ Bk (n) that allows to compute the distribution ps (n) in a recursive way using any path connecting a fixed point n0 with a generic point n in the surface |n| = N. Then we define an internal interaction energy V (n), Ak (n) + Dk− Bk (n) − (A14) Dk V (n) = − ln 1 − − Ek−1 Ek+ Bk (n)
E(t) =
Eq. (A7) has the form of a discrete Kolmogorov forward equation, p˙ n (t) = −
Jk (n)
and an internal energy as
= (−πk−1,k (n)nk−1 + πk,k−1 (n)nk )pn (t) + Dk− πk,k−1 (n)nk pn (t).
By using the relation
p˙ n (t) = −
This remark allows to apply the average field approximation in a local way when we have multimodal distribution. To characterize the properties of the stationary solution of the CME we first analyze the D.B. case. Using the definition (A10), the D.B. equilibrium can be written in the form
(A5)
As a consequence if f(n) vanishes at the boundary points we have Dk± f (n) = 0. (A6)
m
nk Dj+ Jj
Dk+ (nk + 1)Jk (n) −
to any
Dk±
j =1
m j =1
(A3) Dk±
m
for any real function f(n) and g(n) with n ∈ N m . Dk± Dh± = Dh± Dk± .
˙ t) = − nk p(n,
V (n)p(n, t)
so that the equilibrium distribution has the MaxwellBoltzmann form (9). In a NESS, the stationary chemical currents Jks (n) are different from zero and this can be associated with the existence of an external non-conservative field interacting with the system, that satisfies the relation Aex k (n) =
(A10)
where Ak (n) = −πk−1,k (n)nk−1 + πk,k−1 (n)nk , Bk (n) = πk,k−1 (n)nk . Eq. (A6) is useful to compute the average dynamics for CME (4): if we choose such that the chemical fluxes Jk (n)
(A15)
|n|=N
Jks (n) , ps (n)
(A16)
i.e., the chemical fluxes are proportional to an external non conservative field. Indeed one can also define the internal field ex Ain k (n) = Ak (n) − Ak (n) that satisfies the conditions s Dk− Bk (n)p s (n) = Ain k (n)p (n),
which follows from the definition (A10) and can be associated to the internal interactions. In such a case the stationary
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.189.205.30 On: Wed, 10 Dec 2014 20:35:55
de Oliveira et al.
065102-10
J. Chem. Phys. 141, 065102 (2014)
solution satisfies
Dk− ln ps (n) = −Dk− V in (n) = ln 1 −
− Ain k (n) + Dk Bk (n) − Ek−1 Ek+ Bk (n)
(A17)
using (A2) we write − Ek+ pn (t))Dk− Bk (n) + Bk (n)Dk− pn (t), Dk− Bk (n)pn (t) = (Ek−1 (A25) putting (A25) in (A24),
and corresponds to the Maxwell-Boltzmann distribution associated to the internal energy V in (n). The stationary condition for the modulated drift (defined in Eq. (14)) reads m
Dk+ Aˆ k (n)p s (n) + Dk− B)k(n)ps (n)
=
+ Bk (n)Dk− pn (t).
Dk+ Ak (n)p s (n) + Dk− B)k(n)pˆ s (n) m
ˆ s (n) = 0 Dk+ Aex k (n)p
(A18)
−
k=1
and it is satisfied if we set pˆ s (n) = ps (n) since the first term vanishes as a consequence of Eq. (8), whereas the second term becomes m m s Dk+ Aex (n)p (n) = λ Dk+ Jks (n) = 0 λ k k=1
due to the definition of vector reads
S˙in = −
The new stationary current
|n|=N
S˙in = −
p˙ n (t) −
|n|=N
=−
m
=
m
+ Ek−1 Ek− Jk (n, t)
(A27)
.
Dk+ pn pn
+ Ek−1 Ek− Jk (n, t)
× =
+ − Ek−1 Ek (Jk (n) − Ak (n)pn (t)) + pn (t)Ek−1 Ek− Bk (n)
m k=1 |n|=N
taking its time derivative, we obtain S˙in = − p˙ n (t)(1 + ln pn (t)),
+ Ek−1 Ek− Bk (n)
k=1 |n|=N
The fluxes Jk can be related to the entropy production, Sin = − pn (t) ln pn (t) (A19)
+ Ek− Dk− Bk (n) Ek−1
k=1 |n|=N
s s Jˆks (n) = Jks (n) − λAex k (n)p (n) = (1 − λ)Jk (n).
|n|=N
we use the
Therefore, we can write the entropy production (A21) as
k=1
Aex k (n).
(A26)
E + E − (J (n) − A (n)pn (t)) Dk+ pn (t) = − k−1 k k + − k pn (t) pn (t)Ek−1 Ek Bk (n)
k=1
−λ
Dk+ pn (t) pn (t)
To find an explicit expression for relation (A1)
k=1 m
− Jk (n) − Ak (n)pn (t) = (Ek−1 Ek+ pn (t))Dk− Bk (n)
+ Ek−1 Ek− Bk (n)
+ (Ek−1 Ek− Jk (n))2
+ pn (t)Ek−1 Ek− Bk (n)
m
−
−
+ Ek−1 Ek− Dk− Bk (n)
+ Ek−1 Ek− Jk (n, t)
k=1 |n|=N
p˙ n (t) ln pn (t)
(A20)
|n|=N
×
+ + Ek−1 Ek− Ak (n)pn (t) + pn (t)Ek−1 Ek− Dk− Bk (n)
+ pn (t)Ek−1 Ek− Bk (n)
(A28)
p˙ n (t) ln pn (t),
|n|=N
and shifting the sum indexes, we get the expression
replacing p˙ n (t) by CME (A7) we have S˙in =
m
Dk+ Jk (n, t) ln pn (t),
(A21)
k=1 |n|=N
m
+ Ek−1 Ek− Jk (n, t)Dk+ ln pn (t).
Dk+
1, we can approximate ln pn (t) as + − E E p D+p D+p Dk+ ln pn = ln k−1 k n = ln 1 + k n k n . pn pn pn (A23) To determine Dk+ ln pn (t) we consider the definition of currents (A10), If
Jk (n) − Ak (n)pn (t) = Dk− Bk (n)pn (t),
×
Jk2 (n) − Jk (n, t) − Bk (n)Ek−1 Ek+ pn (t)
− Ek+ pn (t) Ak (n)pn (t) + Dk− Bk (n)Ek−1 − Bk (n)Ek−1 Ek+ pn (t)
(A22)
k=1 |n|=N
Dk+ pn
m k=1 |n|=N
where we can use the property (A2) and get the explicit form S˙in = −
S˙in =
(A24)
. (A29)
− If N is not too small we can approximate Ek−1 Ek+ pn (t) pn (t) and we get the relation
S˙in =
m k=1 |n|=N
Jk2 (n) Bk (n)pn (t)
Ak (n) + Dk− Bk (n) . Jk (n, t) − Bk (n) k=1 |n|=N m
(A30)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.189.205.30 On: Wed, 10 Dec 2014 20:35:55
de Oliveira et al.
065102-11
J. Chem. Phys. 141, 065102 (2014) 1 H.
APPENDIX B: PROPERTIES OF THE CHEMICAL CURRENTS FOR THE LINEAR CME
We have Dx Jxs (nx , ny )ps (nx , ny ) + Dy Jys (nx , ny )ps (nx , ny ) = 0 (B1) and the divergence of non-equilibrium fluxes are orthogonal to the gradient of probability if s s Dx ps (nx , ny ) + Jy+1 Dy ps (nx , ny ) = 0. Jx+1
(B2)
So we should compute how the operators Dx and Dy act on the stationary solution ps (nx , ny ), (N − nx − ny )λA λ−1 C s − 1 ps (nx , ny ) Dx p (nx , ny ) = nx + 1 (B3) and (N − nx − ny )λB λ−1 C s − 1 ps (nx , ny ). Dy p (nx , ny ) = ny + 1 (B4) We finally can determine that (B2) is verified, s s Dx ps (nx , ny ) + Jy+1 Dy ps (nx , ny ) Jx+1 (N − nx − ny )λA λ−1 C = −1 nx + 1 πAB nx ny λ−1 πBA nx ny λ−1 A λC B λC − × (N − nx − ny ) (N − nx − ny ) −1 + πCA nx λA λC − πAC nx ps (nx , ny )
+ ×
(N − nx − ny )λB λ−1 C ny + 1 πAB nx ny λ−1 B λC (N − nx − ny )
−
−1
πBA nx ny λ−1 A λC
(N − nx − ny ) s + πCB ny λ−1 λ − π n BC y p (nx , ny ) = 0. B C
(B5)
V. Westerhoff and B. O. Palsson, Nat. Biotechnol. 22, 1249 (2004). 2 H. Qian and T. C. Reluga, Phys. Rev. Lett. 94, 028101 (2005). 3 I. R. Epstein and K. Showalter, J. Chem. Phys. 100, 13132 (1996). 4 J. D. Murray, Mathematical Biology (Springer, 2002), Vol. 2. 5 J. Elf and M. Ehrenberg, Syst. Biol. 1, 230 (2004). 6 E. Giampieri, D. Remondini, L. de Oliveira, G. Castellani, and P. Lió, Mol. BioSyst. 7, 2796 (2011). 7 E. L. Bienenstock, L. N. Cooper, and P. W. Munro, J. Neurosci. 10, 79 (1995). 8 N. Intrator and L. N. Cooper, Neural Networks 5, 3 (1992). 9 G. Castellani, A. Bazzani, and L. Cooper, Proc. Natl. Acad. Sci. U.S.A. 106, 14091 (2009). 10 A. Bazzani, G. C. Castellani, E. Giampieri, D. Remondini, and L. N. Cooper, J. Chem. Phys. 136, 235102 (2012). 11 L. Xu, H. Shi, H. Feng, and J. Wang, J. Chem. Phys. 136, 165102 (2012). 12 X. Wang and X. Wang, Nucleic Acids Res. 34, 1646 (2006). 13 H. Ge and H. Qian, Phys. Rev. E 81, 051133 (2010). 14 H. Qian, Annu. Rev. Phys. Chem. 58, 113 (2007). 15 H. Qian, J. Phys. Chem. B 110, 15063 (2006). 16 H. Qian and L. M. Bishop, Int. J. Mol. Sci. 11, 3472 (2010). 17 H. Qian, Phys. Rev. E 65, 016102 (2002). 18 H. Qian and E. L. Elson, Biophys. Chem. 101, 565 (2002). 19 U. Seifert, Rep. Prog. Phys. 75, 126001 (2012). 20 R. Zia and B. J. Schmittmann, J. Stat. Mech. 2007, P07012. 21 R. Zia and B. J. Schmittmann, J. Phys. A: Math. Gen. 39, L407 (2006). 22 N. van Kampen, Stochastic Processes in Physics and Chemistry (NorthHolland Personal Library, 2007). 23 C. Gardiner, Handbook of Stochastic Methods: For Physics, Chemistry & The Natural Sciences (Springer, Berlin, 2004). 24 J. Schnakenberg, Rev. Mod. Phys. 48, 571 (1976). 25 U. Seifert, Phys. Rev. Lett. 95, 040602 (2005). 26 T. Schmiedl and U. Seifert, J. Chem. Phys. 126, 044101 (2007). 27 A. De Martino, D. De Martino, R. Mulet, and G. Uguzzoni, PLoS One 7, e39849 (2012). 28 M. Polettini, Phys. Rev. E 84, 051117 (2011). 29 B. R. Irvin and J. Ross, J. Chem. Phys. 89, 1064 (1988). 30 T. L. Hill, J. Theor. Biol. 10, 442 (1965). 31 T. L. Hill, J. Chem. Phys. 54, 34 (1971). 32 T. L. Hill, Proc. Natl. Acad. Sci. U.S.A. 85, 5345 (1988). 33 S. R. de Groot and P. Mazur, Non-Equilibrium Thermodynamics (Dover Books on Physics, 1984). 34 L. Jiu-Li, C. Van den Broeck, and G. Nicolis, Z. Phys. B: Condens. Matter 56, 165 (1984). 35 P. Gaspard, C. R. Phys. 8, 598 (2007).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.189.205.30 On: Wed, 10 Dec 2014 20:35:55
