Glucans monomer-exchange dynamics as an open chemical network
GGlucans monomer-exchange dynamics as an open chemical network
Riccardo Rao, David Lacoste, and Massimiliano Esposito Complex Systems and Statistical Mechanics, Physics and Materials Science Research Unit, University of Luxembourg,L-1511 Luxembourg, Luxembourg Laboratoire de Physico-Chimie Théorique, UMR CNRS Gulliver 7083, ESPCI - 10 rue Vauquelin, F-75231 Paris,France (Dated: 24 September 2018. Published in
J. Chem. Phys.
DOI: 10.1063/1.4938009)
We describe the oligosaccharides-exchange dynamics performed by so-called D-enzymes on polysaccharides.To mimic physiological conditions, we treat this process as an open chemical network by assuming some ofthe polymer concentrations fixed (chemostatting). We show that three different long-time behaviors mayensue: equilibrium states, nonequilibrium steady states, and continuous growth states. We dynamically andthermodynamically characterize these states and emphasize the crucial role of conservation laws in identifyingthe chemostatting conditions inducing them.PACS numbers: 05.70.Ln, 05.70.-a, 82.20.-w
I. INTRODUCTION
Biological systems use large and branched chains ofbasic sugars, called polysaccharides, to store energy . Glucans such as glycogen and starch are polysaccharideswhose building blocks are
D-glucose monosaccharides.Despite the apparent simplicity of their constituents,their metabolism involves several chemical steps, eachperformed by a specific set of enzymes . Interestingly,some of these catalysts lack specificity regarding the re-action they catalyze or the substrates they act on . Anexample is provided by (1 → (EC2.4.1.25), also called D-enzymes , which act on pairs ofglucans regardless of their size . Specifically, D-enzymescatalyze the transfer of groups of glycosyl residues froma donor glucan to an acceptor glucan . Experimentalevidences highlight the presence of bonds between gly-cosyl residues which are not cleaved by D-enzymes —atleast not over physiological time scales . These bondsare called forbidden linkages . In this way, D-enzymestransfer segments of glucan chains containing one or moreforbidden linkages, and the transfer of segments contain-ing one forbidden linkage are the most probable . Also,each glucan chain is characterized by a reducing-end glu-cose which is not transferred by D-enzymes . Hence,glucans made of just the reducing end can act only asacceptor in the transfer.Qualitatively, D-enzymes process medium-size glucansby disproportionating them into unit-size and big-sizeglucans . Since their transfers reactions are neutralenergetically , entropy is the main driving force in thissystem. In closed conditions, this system evolves towardsan equilibrium state maximizing the entropy .In this paper we consider a simplified kinetic descriptionof the D-enzyme’s action on glucans, which we treat as achemical network. Since metabolic processes should bethought of as part of an open system continuously fed fromthe environment, we mimic these physiological conditionsby introducing chemostats (i.e. species whose concentra-tions are kept constant by the environment). Our goalis to characterize the dynamical and thermodynamical implications of treating the action of the D-enzymes onglucans as an open chemical network. In the frameworkof deterministic chemical networks endowed with massaction kinetics , we prove that chemostatting can inducethree different types of long-time behaviors: equilibrium,non-equilibrium steady state, and continuous growth. Theequilibrium state corresponds to the stationary concen-tration distribution in which the concentration currentsalong each reaction pathway vanishes (detailed balanceproperty ). Non-equilibrium steady states refer to sta-tionary distributions violating detailed balance. Hence,contrary to equilibrium states, a continuous and steadyflow of mass circulates across the network. Finally, thecontinuous growth regime we observed corresponds toa non-stationary state characterized by continuous andsteady flow of mass entering the network, and resulting inits continuous growth. We emphasize the dynamical andthermodynamical role of conservation laws and emergentcycles in identifying the chemostatting conditions leadingto these states. We are thus able to confirm the generalrelation between the number of chemostatted species andthe number of independent thermodynamical forces—or affinities —found in Ref. . Despite the simplicity of ourdescription, the closed system results found in Ref. arereproduced and the qualitative disproportionating behav-ior of D-enzymes is captured by our (chemostatted) opensystem description.The plan of the paper is as follows: in sec. II the kineticmodel is established and the related rate equation descrip-tion for the concentration of polysaccharides is introduced.In sec. III the chemostatting conditions leading to non-equilibrium steady states rather than equilibrium onesare found. For this purpose, both the conservation laws ofthe dynamics and the emergent cycles of the network areanalyzed. The dissipation of the non-equilibrium steadystate is also studied. The network’s conservation lawsidentified in sec. III A are used in sec. IV to derive thesteady-state concentration distributions for different num-bers of chemostats. The explosive asymptotic behavioris described in sec. V. Conclusions are drawn in sec. VI.Some technical derivations and proofs are provided in the a r X i v : . [ phy s i c s . c h e m - ph ] J a n κκ (a)(b) FIG. 1: (a) The typical monomer-exchange reaction de-scribing the action of D-enzymes on glucan chains. (b)The attachment of free monomers to other species is notallowed. appendices.
II. THE KINETIC MODEL
The action of D-enzymes is modeled as follows (see alsoFig. 1). Glucans are treated as polymers whose monomersrepresent single transferable segments. Hence, each glucanis identified by its number of monomers, or equivalentlyby its monomeric mass. The enzymatic steps performedby the D-enzymes in order to achieve the transfer are notexplicitly described—they are coarse-grained—, and wedescribe the interaction between two polymers of mass n and m as a mass-exchange process :( n ) + ( m ) κ nm −→ ( n + 1) + ( m − , for n ≥ , m ≥ , (1)where κ nm denotes the related coarse-grained rate con-stant. Transfers of oligosaccharides longer than onemonomeric unit are less probable , and are not consideredin our description. We take into account the presence ofnon-transferable units by imposing the size of the donorglucan ( m ) to be greater than one .Let us note that each reaction is reversible because thebackward path is already included in (1) (it is realizedby replacing n → m − m → n + 1 in the aboveexpression). Furthermore, the constraint on the minimalsize of the donor molecules imposes that m ≥
2. Sincewe describe the glucans as linear polymers, and sinceD-enzymes do not discriminate the size of the polymers,we assume a constant kernel for the reactions: κ mn = κ , ∀ n ≥ , ∀ m ≥
2. This assumption is based on the evidencethat the free-enthalpy release resulting from any reactionis almost vanishing . Indeed, for any bond cleaved, anew one of the same kind will be formed.Assuming a large and well stirred pool of interactingpolymers, the evolution of the system is well described byreaction rate equations . According to this mean-fielddescription, the molar concentration of polymers of mass k at time t , Z k = Z k ( t ), satisfies the following first orderdifferential equations˙ Z k = 12 X n ≥ m ≥ ∇ knm (cid:0) J + nm − J − nm (cid:1)| {z } ≡ J nm , for k ≥ . (2) The factor in front of the summation takes into accountthat summing over all n ≥ m ≥ . ∇ knm represents the element ofthe stoichiometric matrix related to the species of mass k and to the reaction involving an acceptor and a donorpolymer of mass n and m , respectively. The reactionscheme in (1) implies that ∇ knm = δ kn +1 + δ km − − δ kn − δ km , (3)where δ ji represents the Kronecker delta. Assuming a mass action kinetics , the forward (denoted by +) and thebackward fluxes ( − ) can be written as: J + nm = κZ n Z m , J − nm = κZ n +1 Z m − , (4)where Z n , denotes the concentrations of the polymersof size n . To simplify the following discussion, we willuse the Einstein summation notation: upper indexedquantities represent vectors, lower indexed ones covectors,and repeated indexes implies the summation over all theallowed values for those indexes (1 ≤ n ≤ n max and2 ≤ m ≤ m max , or 1 ≤ k ≤ k max , where n max , m max and k max are finite in closed systems but infinite in openones). To avoid confusion, exponents will always act onparenthesis ( e.g. ( a ) n denotes the quantity a to the power n ).The rate equations (2) assume the following form whenthe expressions for both the stoichiometric matrix (3) andthe fluxes (4) are considered˙ Z = κZ (cid:0) Z − Z (cid:1) + κZ Z , ˙ Z k = κZ (cid:0) Z k +1 − Z k + Z k − (cid:1) ++ κZ (cid:0) Z k − Z k − (cid:1) , for k ≥ , (5)where Z ≡ P k max k =1 Z k denotes the total concentration.The second term in the right hand side of (5) arisesfrom the constraint that the donor species cannot bemonomers (see Fig. 1b).To model the open system we now assume that theenvironment keeps the concentrations of some speciesconstant by refilling the consumed ones and eliminatingthe produced ones, see Fig. 2. We call these species chemostats and we denote them with the indices k y ∈ Ω Y , where Ω Y ⊂ N represents a subset of all species. Theremaining (variable) species are explicitly denoted by k x .By definition, the chemostats’ concentrations mustremain constant, ˙ Z k y = 0. The rate of chemostattedmolecules consumed by the reactions in the networkmust therefore be balanced by the rate of chemostattedmolecules injected/rejected from the system. The rate ofinjection/rejection of the k y -th chemostat is quantified bythe external currents , whose expression is I k y = 12 ∇ k y nm (cid:0) J + nm − J − nm (cid:1) == κZ (cid:0) Z − Z (cid:1) + κZ Z if k y = 1= κZ (cid:0) Z k y +1 − Z k y + Z k y − (cid:1) + κZ (cid:0) Z k y − Z k y − (cid:1) if k y ≥ . (6) system environment κκ FIG. 2: Pictorial representation of a reaction involving achemostat. When a reaction produces a chemostat (herea dimer), the environment extracts one molecule of thisspecies from the system (dotted light green reaction). Onthe other hand, when a chemostat reacts, a new moleculeis injected into the system (dashed dark green reaction).
III. STEADY STATES: CONSERVATION LAWS,CYCLES AND DISSIPATION
Three different types of long-time behaviors havebeen identified for our kinetic model: equilibrium, non-equilibrium steady state and continuous growth. We startby focusing on the chemostatting conditions leading toequilibrium or non-equilibrium steady states. The exis-tence and uniqueness of the steady state is currently apriori assumed.Closed systems always reach an equilibrium steadystate defined by ˙ Z k x eq = 0 , ∀ k x and J nm eq = 0 , ∀ n, m .Their dynamics is constrained by conservation laws ,which fully characterize the equilibrium concentration dis-tribution. Chemostatting generic chemical species maybreak these conservation laws and may create chemicalforces—also called affinities . The appearance of affini-ties is directly related to that of so-called emergent cycles ,through which the external chemical forces can act. Infinite chemical networks, if no emergent affinity arisesfrom the chemostatting procedure, the system will al-ways relax to a unique equilibrium state compatible withthe chemostats and the non-broken conservation laws .When emergent cycles—or equivalently affinities—are gen-erated, the system may evolve towards a non-equilibriumsteady state defined by ˙¯ Z k x = 0 , ∀ k x and ¯ J nm = 0 (non-equilibrium steady state quantities are denoted by anoverbar in the text). In the following subsections weanalyze how the closed system’s conservation laws andemergent cycles are modified by the gradual increase ofthe number of chemostatted chemical species. In the lastsubsection we relate these to the dissipation in the system. A. Conservation Laws
Conservation laws denote the presence of physical quan-tities which are conserved during the evolution of thesystem, the so-called components . In general, they canbe identified from the cokernel space of the stoichiomet-ric matrix . Indeed, if l k ∈ coker ∇ , namely if l k ∇ knm = 0, the scalar l k Z k is conserveddd t (cid:0) l k Z k (cid:1) = l k ˙ Z k = 12 l k ∇ knm (cid:0) J + nm − J − nm (cid:1) = 0 . (7)For the closed system, the equation leading to theconservation laws is l kn +1 − l kn = l km − l km − , for 1 ≤ n ≤ n max = k max − ≤ m ≤ m max = k max . It exhibitsthe following solutions: l (1) k = α and l (2) k = α · k (where α is an arbitrary constant, which is taken as one whenexpressing the components), which correspond to theconservations of the total concentration Z ≡ P k max k =1 Z k and the total mass M ≡ P k max k =1 kZ k , respectively. Hence, k max = M − Z + 1.However, when the system is opened by settingchemostats, the relevant stoichiometric matrix becomesthe stoichiometric submatrix of the variable species: ∇ k x nm .Also, k max = ∞ . No matter what the sizes of thechemostatted glucans are, neither the total concentrationconservation law l k x = α nor the total mass conservationlaw l k x = αk x survives (i.e. they are not anymore ele-ments of the cokernel space of ∇ k x nm ). We therefore saythat the total mass and the total concentration are brokenconservation laws . Nevertheless, when just one chemostatis present, Ω Y ≡ { k y } , a new conservation law emerges: l (3) k x = α ( k x − k y ) . (8)Hence, the system exhibits just one (net) broken conser-vation law. It corresponds to the component q = M − k y Z, (9)which can assume any value in R and takes into accountthat the total mass can change in the system only bymultiples of the chemostat mass, k y . In presence of morethan one chemostat, no conservation law survives.The components derived in this section— M and Z for the closed system and q for the network with onechemostat—will be used to characterize the steady statedistribution in sec. IV. B. Emergent cycles
A cycle represents a finite set of reactions which leavethe state of the network unchanged. Algebraically they arerepresented as vectors c nm and they belong by definitionto the kernel space of the stoichiometric matrix ( c nm ∈ ker ∇ ): ∇ knm c nm = 0.The steady-state currents satisfy ∇ knm ¯ J nm = 0 and canalways be written as linear combinations of cycles. Thecycle space of our polymers system is however infinite di-mensional and its complete characterization is of little use.However, in order to characterize non-equilibrium steadystates only the emergent cycles —those cycles that mayappear when chemostatted species are introduced—areneeded . Physically, they represent cyclic transforma-tions leaving the variable species k x unchanged, but whichwould change the concentrations of the chemostats k y ifthey were not kept constant and contribute to the externalcurrents.An emergent cycle ( γ nm ) is thus defined by ( ∇ k x nm γ nm = 0 , ∇ k y nm γ nm = ν k y γ = 0 for at least one k y , (10)where n ν k y γ o k y ∈ Ω Y denotes the amount of chemostatsof mass k y injected (minus sign) or rejected (plus sign)from the chemical network during the transformation γ nm .These quantities cannot take arbitrary values, due to theconstraints imposed by the conservation laws of ∇ knm .Indeed, for any conservation law, l ( i ) k , a constraint of thefollowing form holds l ( i ) k y ν k y γ = l ( i ) k y ∇ k y nm γ nm = 0 . (11)Taking into account the total concentration l (1) k = α and total mass l (2) k = αk conservation laws, derived insec. III A (the emergent conservation law l (3) k is a linearcombination of the first two on the whole set of speciesindexes), we obtain the following constraints P k y ν k y γ = 0 P k y k y ν k y γ = 0 . (12)Non-trivial solutions of this set of equations signal thepresence of emergent cycles, and thus of independentaffinities, which read A γ = 12 X nm γ nm ln Y k y ( Z k y ) −∇ nmk y . (13)The set of linearly independent solutions of (12) gives thenumber of independent emergent cycles in the chemostat-ted chemical network. If we normalize this set so to havethe smallest non-vanishing integer values for ν k y γ , thesevalues indicate the number of chemostatted species whichare introduced in or rejected from the system in preciselyone (emergent) cyclic transformation.For less than three chemostats, only trivial solutions of(12) exist and therefore no emergent cycle appears. Forthree chemostats, we obtain one emergent cycle charac-terized by the following normalized values for ν k y : ν k y1 = k y3 − k y2 ,ν k y2 = k y1 − k y3 ,ν k y3 = k y2 − k y1 , (14)where k y1 , k y2 and k y3 represent the masses of thechemostats. For any additional chemostat we obtainan additional emergent cycle, each characterized by itsvalue for the coefficients ν k y . C. External Currents and Dissipation
We now show that at steady state, the emergent cy-cles determine the external currents ¯ I k y and the entropyproduction rate Σ.We first observe that the steady-state external currents¯ I k y are in general linear combination of the coefficients ν k y γ i and must satisfy the same constraints (eq. 12). Indeed,the steady-state equations in presence of chemostats ( ∇ k x nm ¯ J nm = 0 , ∇ k y nm ¯ J nm = ¯ I k y (15)are equivalent to eq. (10): the emergent cycles γ nm aresubstituted by the steady state currents ¯ J nm and thecoefficients ν k y by the steady-state external currents ¯ I k y .Thereby, if no cycle emerges due to the chemostats, thesteady-state external currents ¯ I k y are vanishing, providedthat the steady state exists. The system is then at equi-librium.The dissipation at steady state is intimately relatedto the external currents . Indeed, the (non-negative)entropy production rate for our chemical reaction networkcan be written asΣ ≡ X nm J nm R ln J + nm J − nm == − X k x ˙ Z k x R ln Z k x Z k x eq | {z } ≡ Σ X − X k y I k y R ln Z k y Z k y eq | {z } ≡ Σ Y , (16)where R is the gas constant. At the steady state, theinternal species’ contribution Σ X always vanishes. Hence,the dissipation is characterized by the contribution dueto the chemostats Σ Y , which is non-vanishing if the set ofsteady-state external currents ¯ I k y is also non-vanishing.We also mention that, the steady state entropy productioncan be expressed as the sum along a set of independentemergent cycles of products of affinities (13) and emergentcycle currents J γ : ¯Σ = P γ A γ J γ .Summarizing, the conservation laws provide us withboth the components— which are useful for expressing thesteady state distributions—and the constraints (eq. 12) onthe emergent cycles of the network (eq. 10). Due to theseconstraints, the first emergent cycle appears in the systemwith three chemostats. For any additional chemostat anadditional independent cycle emerges. Through thesecycles the environment exerts chemical forces, which aregenerated by the chemostats concentrations. The above-defined external currents result from these forces andcharacterize the dissipation.We emphasize that the relation between the numberof chemostats s Y , of net broken conservation laws b , andof emergent cycles a , is in perfect agreement with thegeneral result obtained for finite-dimensional phase spacein Ref. stating that s Y = b + a (17) number of broken independent asymptoticchemostats, s Y c. laws, b affinities, a behavior0 0 0 ES1 1 0 ES2 2 0 ES/growth3 2 1 NESS/growth4 2 2 NESS/growthTABLE I: Summary of the behaviors of our model for dif-ferent numbers of chemostats (ES stands for “equilib-rium state” whereas NESS for “non-equilibrium steadystate”). The number of broken conservation laws and in-dependent affinities are also reported. The growth stateoccurs whenever the concentration of the largest chemo-stat is larger than the concentration of the smallest one:( Z k y larger ≥ Z k y1 ). These results are summarized in Tab. I.
IV. THE STATIONARY DISTRIBUTIONS
We now use the components introduced in the previoussection to derive the steady-state concentration distribu-tion for different number of chemostats. The conditions onthe chemostats’ concentrations not leading to the steadystate solution are also identified.From the steady-state equations corresponding to (5)and from the equations for the external currents (6), wecan write a general expression for the steady-state concen-trations as a function of the concentration of monomers,¯ Z , the fraction of polymers larger than monomers,¯ r ≡ − ¯ Z / ¯ Z , and the chemostats fluxes, ¯ I k y :¯ Z k = ¯ Z (¯ r ) k − ++ X k y ∈ Ω Y ¯ I k y κ − (¯ r ) k − k y − ¯ r Θ ( k − k y − , (18)where Θ( · ) represents the discrete step function (we referthe reader to appendix A for details). Here, the number ofchemostats is arbitrary, and since the external currents atsteady state satisfy the same constraints as in (12), only s Y − Z , ¯ r and¯ I k y will be expressed in terms of the components and ofthe chemostats’ concentrations. A. Closed system
As previously discussed, the closed system exhibitsthe following components: Z = P k max k =1 Z k and M = P k max k =1 kZ k . In order to express the equilibrium distribu-tion algebraically as function of Z and M we considerthe following limit M (cid:29) Z . In this way k max ∼ ∞ and k s t ead y (cid:45) s t a t e c onen t r a t i on , Z k FIG. 3: Equilibrium concentration distribution for theclosed system of monomers-exchanging polymers at dif-ferent values of the total concentration Z and total mass M . The dark blue bar plot refers to the choice Z = 10 and M = 15, while the light blue one to Z = 10 and M = 55. imposing Z = P ∞ k =1 Z k and M = P ∞ k =1 kZ k on the ex-pression (18) we can write ¯ Z and ¯ r as functions of Z and M . Hence¯ Z k = ( Z ) M (cid:18) − ZM (cid:19) k − . (19)Fig. 3 shows the behavior of this distribution for differentvalues of Z and M . As expected, the higher the ratiobetween the mass and the concentration M (cid:29) Z , thebroader the distribution. Remark.
The equilibrium distribution we obtainedfrom our dynamical description is equivalent to the resultobtained using maximum entropy approaches and is con-sistent with experimental observations . The equivalenceis inferred by comparing eq. (19) with eq. (1), (3) and (4)in Ref. . B. Open system: 1 chemostat
Introducing a chemostat breaks the concentration andmass conservation laws, but a new one arises (8). Asa result, no affinity appears ( s Y = 1, b = 1 and a =0) and the system evolves towards an equilibrium statecompatible with the chemostat concentration Z k y and thevalue of the component q (9) (the steady-state externalcurrent vanishes, ¯ I k y = 0). Also, since the system is nowopen, k max is infinite.Imposing the constraints on the expression of the steadystate (18), namely q = ¯ Z − k y (1 − ¯ r )(1 − ¯ r ) Z k y = ¯ Z (¯ r ) k y − , (20)we can express the variables ¯ Z and ¯ r numerically asfunctions of q and Z k y , and obtain the equilibrium—exponential—distribution as a function of q and Z k y . C. Open system: 2 chemostats
From two chemostats on, the infinite dimension of thesystem starts to play a role. As discussed in the previoussection, two chemostats are not enough to drive the net-work towards a non-equilibrium steady state ( s Y = 2, b = 2 and a = 0): I k y1 = 0 and I k y2 = 0, where k y1 and k y2 represent the masses of the two chemostats( k y1 < k y2 ). Thus, imposing the known values of thechemostat concentrations on the expression (18) leads to ( Z k y1 = ¯ Z (¯ r ) k y1 − Z k y2 = ¯ Z (¯ r ) k y2 − , (21)which only admits physical solutions if Z k y1 > Z k y2 . Inthis case, from (21) we obtain the equilibrium distribution¯ Z k = Z k y1 (cid:18) Z k y2 Z k y1 (cid:19) k − k y1 k y2 − k y1 , (22)which is broader the smaller Z k y1 − Z k y2 is or the larger k y2 − k y1 is. When Z k y1 ≤ Z k y2 the equilibrium concen-tration distribution becomes an increasing exponentialwhich cannot be reached. As a result the system willenter a regime of continuous growth aimed at reachingthat state (which we analyze in sec. V). D. Open system: 3 chemostats
Three is the minimum number of chemostats ableto drive the system in a non-equilibrium steady state(sec. III B). Indeed, a class of emergent cycles appears( s Y = 3, b = 2 and a = 1) and the system exhibits aset of non-vanishing external currents. If we impose thevalues for the chemostats’ concentrations on the generalexpression for the steady state (18), we obtain: ¯ Z k y1 = ¯ Z (¯ r ) k y1 − ¯ Z k y2 = ¯ Z (¯ r ) k y2 − + ¯ I k y1 κ − (¯ r ) k y2 − k y1 − ¯ r ¯ Z k y3 = ¯ Z (¯ r ) k y3 − + ¯ I k y1 κ − (¯ r ) k y3 − k y1 − ¯ r ++ ¯ I k y2 κ − (¯ r ) k y3 − k y2 − ¯ r . (23)As discussed in sec. III C, the external currents ¯ I k y aresubject to the same constraints as the emergent cycles,and can be written as linear combinations of them. Sincewe have one class of emergent cycles, characterized by the ν k y values in (14), we have that¯ I k y i = ¯ Iν k y i , i = 1 , , , (24) where ¯ I ∈ R determines the exact value of the fluxes. Asfor two chemostats, the set of equations in (23), in thevariables ¯ Z , ¯ r and ¯ I , does not exhibit physical solutionsif the concentration of the largest chemostat is higher thanthe one of the smallest one, i.e. Z k y1 ≤ Z k y3 . On theother hand, whenever the above condition is not fulfilled,the stationary solution is unique and stable (appendix B).Solving the system (23) numerically, we obtain the valuesof ¯ Z , ¯ r and ¯ I given Z k y1 , Z k y2 and Z k y3 . In Fig. 4the distribution is shown for different values of theseconcentrations.The chemostat concentrations also determine the signof the related fluxes: if the concentration of the secondchemostat lies above the equilibrium distribution obtainedby the first and third one, we have a continuous flow ofmass from the intermediate chemostat towards the exter-nal ones ( ¯ I >
0, Fig. 4a).
Vice versa , if the concentrationof the second chemostat lies below the equilibrium distri-bution obtained by the first and the third one, we havea continuous flow of mass from the smallest and largestchemostats towards the intermediate one ( ¯
I <
0, Fig. 4b).Importantly, whatever physical value Z k y1 , Z k y2 and Z k y3 assume, the system cannot exhibit a condition in which anet flux of matter from the largest species to the smallestone occurs. This is clear by looking at the ν k y -values in(14) used to express ¯ I k y i , eq. (24): the sign of ν k y1 and ν k y3 are always the same, and opposite to the one of ν k y2 . E. Open system: more chemostats
Going on adding chemostats, new independent classesof emergent cycles appear. The procedure for determin-ing the steady-state distribution is equivalent to thatdiscussed is subsec. IV C and IV D. In these two caseswe proved that when the largest chemostat has a con-centration greater or equal to that of the smallest one,the system does not reach a steady state. The sameexact behavior has been observed numerically for morechemostats, hence we speculate that this property holdsfor an arbitrary number of chemostats.As a final remark, we point out that the steady-statedistributions do not depend on the value of the rateconstant κ . Indeed, solving the equations (20), (21), and(23) for ¯ Z , ¯ r and ¯ I k y /κ , we obtain them as functionsof the components and the chemostats’ concentrations.Since the latter do not depend on κ , the same holds for¯ Z , ¯ r and ¯ I k y /κ . As a corollary ¯ I k y is proportional to κ and the same holds true for the entropy production (16). V. ASYMPTOTIC GROWTH REGIME
We mentioned in the previous section that the systemdoes not exhibit a steady state when the concentrationof the largest chemostat exceeds that of the smallest k s t ead ys t a t e c onen t r a t i on , Z k (a) k s t ead ys t a t e c onen t r a t i on , Z k (b) FIG. 4: Non-equilibrium steady-state distributions forthe system of monomer-exchanging polymers withthree chemostatted species. In both of the plots, thechemostats—highlighted in green and by the arrows—are k y1 = 2, k y2 = 5 and k y3 = 10. The orientation of the ar-rows denotes the sign of the external fluxes of chemostats:arrows pointing up means chemostats leaving the system, i.e. I k y >
0. The chosen chemostat’s concentrations are:plot (a) Z k y1 = 5, Z k y2 = 7 and Z k y3 = 2; plot (b) Z k y1 = 5, Z k y2 = 1 and Z k y3 = 2. one, Z k y1 ≤ Z k y last —we refer in the text to this con-figuration of chemostats leading to continuous growthas “unbalanced”. The dynamical fixed point moves out-side the region of physical solutions—namely to ¯ r ≥ r →
1. This indicates that the concentration of the singlemonomer species becomes negligible compared to the restof the species. Hence the system grows towards an un-reachable steady state with an exponentially increasingconcentration distribution.Fig. 5a shows the concentration distributions of anunbalanced system at different times before the numericalcut-off (more details are given in the related caption) isreached. These different distributions show that whilethe concentrations of the species between two chemostats t (cid:61) (cid:61) (cid:61) (cid:61) .5t (cid:61) .10 10 20 30 40 50 60 7002468 polymer mass, k c on c en t r a t i on , Z k (a) t (cid:61) (cid:61) (cid:61) (cid:61) .5t (cid:61) .10 10 20 30 40 50 60 70012345 polymer mass, k c onen t r a t i on , Z k (b) FIG. 5: Concentration distributions at different times areshown for system in unbalanced conditions. Different col-ors from red to violet correspond to exponentially increas-ing times. The set of plots is obtained by numerical solu-tion of the differential equation (5). Absorbing boundaryconditions have been chosen, meaning that the concentra-tion at the cut-off—here set to k cut − off = 1000—is zero.We point out that this prescription is safe before the cut-off is reached. In plot (a) we report a system with threechemostats. The chemostat’s masses and the related con-centrations chosen are: Z = 1, Z = 7 and Z = 2. Theconcentrations of the species between the chemostats ba-sically overlap at times t (cid:38) Z = 3 and Z = 4. stabilize to steady values, the concentrations of the specieslarger than the biggest chemostat do not. Hence, thesystem continuously grows trying to populate the infinitesize polymer. This behavior has been observed taking intoaccount different number of chemostats and chemostat’concentrations.In order to characterize this growth algebraically, weconsider a system with monomer and dimer chemostats( k y1 = 1 and k y2 = 2) such that Z k y1 ≤ Z k y2 . (The typi-cal growth obtained numerically in this scenario is shownin Fig. 5b). Since the growth dynamics cannot be solvedexactly, we assume that the asymptotic concentration dis-tribution can be parametrized by the (equilibrium) steadystate expression (18) with time dependent parameters, i.e. Z k ( t ) ’ A ( t ) (cid:0) a ( t ) (cid:1) k − , for k ≥ . (25)where A ( t ) and a ( t ) are unknown real functions of time.To simplify the notation, let us denote the concentrationsof the chemostats by Y ≡ Z k y1 and Y ≡ Z k y2 . Thefunctions A ( t ) and a ( t ) can be determined by means ofthe differential equations for the total concentration Z and the total mass M :˙ Z = − I − I = − κZ ( Z − Y ) − κY Y ˙ M = − I − I = − κZ (2 Z − Y + Y ) − κ Y Y + κY Y , (26)where, Z , M and the concentration of trimers Z assumethe following form when the ansatz (25) is taken intoaccount Z ( t ) ’ A ( t )1 − a ( t ) + Y + Y ,M ( t ) ’ − a ( t ) (cid:0) − a ( t ) (cid:1) A ( t ) + Y + 2 Y ,Z ( t ) ’ A ( t ) . (27)When the equations are expressed in terms of A ( t )and a ( t ), the stream plots for different values of thechemostats’ concentrations show that the ansatz capturesthe non-equilibrium phase transition occurring when thechemostats become unbalanced, Fig. 6. Indeed, for bal-anced chemostats, the system evolves towards a fixedpoint with a lying in ]0 , a = 1 signaling an asymptotic growth regime, seeFig. 6b.The numerical solution for A ( t ) and a ( t ) obtained using(26) and (27) accurately characterizes the asymptoticgrowth. Indeed, as seen in Fig. 7, when comparing theevolution of Z and M obtained from A ( t ) and a ( t ) withthat obtained by solving numerically the rate equations,the former solution overlaps with the latter before thecut-off used in the numerics is reached. We find that thetotal concentration grows linearly with time whereas themass quadratically.Taking into account the ansatz (25), the entropy pro-duction rate (16) becomesΣ ’ RI ln A ( t ) Y ( a ( t )) + RI ln A ( t ) Y a ( t ) , (28)where, I and I can be written in terms of Y , Y , A ( t ) and a ( t ) using eq. (26). The latter is plotted in (cid:72) t (cid:76) A (cid:72) t (cid:76) (a) balanced chemostats (cid:72) t (cid:76) A (cid:72) t (cid:76) (b) unbalanced chemostats FIG. 6: Stream plot of the differential equations (26) ex-pressed in terms of the ansatz functions a ( t ) (abscissa) and A ( t ) (ordinate). When balanced chemostat concentrationsare used, the fixed point lies for values of a ( t ) in ]0 , Y = 4 and Y = 2. Vice versa , when the chemostats are unbalanced( Y = 2 and Y = 4) the fixed point moves outside fromthe physical region ( a ( t ) > Fig. 8, where it is compared with the numerical solutionsfor different cut-offs. The agreement with the numericalsolution is not perfect but captures the linear asymptoticgrowth of the entropy production rate reasonably well.Also, we point out that the unbalanced dynamics shownin Fig. 8 exhibits an initial transient relaxation stageshown in inset.We conclude mentioning that the same ansatz could beused for systems characterized by more chemostats withunbalanced concentrations. Indeed, the growth alwaysinvolves the species larger than the biggest chemostat,whereas the species between chemostats converge fasterto proper steady values. Hence, fixing the concentrationof these latter species, we could assume a growth like(25) for the species larger then the biggest chemostat and t o t a l c on c en t r a t i on , Z (a) time, t t o t a l m a ss , M (b) FIG. 7: Total concentration (a) and total mass (b) as func-tions of time in the asymptotic growth regime. The nu-merical solution obtained using the ansatz (25) is plottedin green (dashed). These plots are compared with numer-ical solutions of the system of differential equations (5)with different cut-offs (blue curves). The chosen chemostatconcentrations are: Y = 3 and Y = 4 while the initialcondition imposed is Z k ( t = 0) = ( ) k . Finally, the cho-sen cut-off concentrations are: k c = 200 (dark blue curve), k c = 500 (blue curve) and k c = 1000 (light blue curve). perform the same analysis. VI. CONCLUSIONS
This paper provides a kinetic description of systemsmade of glucans and processed by the class of enzymesknown as D-enzymes. The action of the enzyme inducesa monomer-exchange process between pairs of glucanswhich are distinguished by their mass or degree of poly-merization. Free monomers are not allowed to attach toother polymers implying that the total concentrationand the total mass are conserved when the system isclosed. The system’s dynamics is ruled by rate equationsfor the polymer concentrations endowed with mass actionkinetics. We mimic physiological conditions by introduc- en t r op y p r odu c t i on r a t e , (cid:83)
0. 1.180220
FIG. 8: Entropy production rate as a function of time inthe asymptotic growth regime. The numerical solution ob-tained using the ansatz (25) is plotted in green (dashed).This plot is compared with numerical solutions of the sys-tem of differential equations (5) with different cut-offs(blue curves). In all the plot, the entropy production rateis given in units of R . The chosen chemostat concentra-tions are: Y = 3 and Y = 4 while the initial conditionimposed is Z k ( t = 0) = ( ) k . The chosen cut-offs k c are:200 (dark blue curve), 500 (blue curve) and 1000 (lightblue curve). Also, the inset shows in greater details theinitial transient relaxation stage. ing chemostats which effectively describe the action ofthe environment by fixing the concentrations of certainglucans. In this scenario, chemostats represent speciesprocessed by the environment. For example, they mayrepresent species which need to be processed and injectedby the environment in the system, analogously, they mayrepresent final products of the metabolic processes whichare taken out of the system. Importantly, chemostattingthe system amounts to open it and introduce drivingforces on the non-chemostatted species.Our main results are summarized in Table I. We identi-fied three types of different long-time behaviors dependingon the chemostatting conditions: equilibrium state, non-equilibrium steady state, and continuous growth of thesystem. The closed system as well as the open systemwith a single chemostat always relax to an equilibriumstate. In presence of two chemostats the system will eitherrelax to equilibrium or turn into a state of continuousgrowth depending on whether or not the concentrationof the largest chemostat is lower than the concentrationof the smallest one. We proved that this latter conditionfor growth holds true for up to three chemostats andconjectured that it is generally true based on numericalevidence. For more than two chemostats, if the concen-tration of the largest chemostat is lower than that ofthe smallest one, the system will reach a nonequilibriumsteady state where the chemostats continuously exchangesmatter across the system. Our results confirm that, evenin the infinite-dimensional chemical network consideredhere, the number of chemostats equals to the number ofbroken conservation laws plus the number of emergent0cycles (see Tab. I). A proof of this equality for finite dimen-sional chemical networks is provided in Ref. . We alsoemphasized the role of the emergent cycles in driving thechemostatted chemical networks towards nonequilibriumsteady states rather than equilibrium states .The metabolism of polysaccharides is a complex processinvolving many steps and several enzymes and its com-plete dynamical characterization is beyond the scope ofthe present paper. We focused on the dynamical character-ization of the disproportionating action of D-enzymes inthe breakdown and synthesis processes of glucans . Underphysiological conditions, it has been pointed out that oneof the possible role of D-enzymes in these processes is toproduce glucans of large sizes (which are then processedby other enzymes) starting from medium sized ones .Importantly, a production of glucose (monomers in ourdescriptions) is expected, too . This disproportionatingbehavior can be reproduced in a (nonequilibrium) steadystate by the three chemostats system depicted in Figure4a. The intermediate high concentration chemostattedglucans represent the species to be processed, while thelow concentration chemostatted glucans represent thespecies to be produced—in this case the small and largeglucans. In this scenario, a continuous flow of interme-diate glucans enters the system and consequently boththe smaller and the larger glucans are steadily producedand expelled from the system (§ IV D). We stress thatthe production of the small glucans follows from the totalconcentration conservation law (§ III A), i.e. the fact thatfree monomers cannot attach to other glucans. As seen insec. IV C, two chemostats are not sufficient to reproducea nonequilibrium steady state.Also, under closed in vitro conditions, the equilibriumdistribution (which has also been analyzed in Ref. andcan be equivalently obtained by means of Maximum En-tropy methods ) agrees with experiments . This meansthat if chemostatting conditions could be implemented in vitro , our predictions could be verified experimentally.Such a procedure would also enable to engineer differentpolymer concentration distributions.The approach we developed could be easily extendedto describe the behavior of more sophisticated forms ofD-enzymes embedding further conservation laws. It isalso relevant to study any type of exchange process oraggregation–fragmentation dynamics in an open systemframework , emphasizing the importance of conserva-tion laws and providing more insights into the mechanismsdriving these processes out of equilibrium. ACKNOWLEDGMENTS
R.R. is grateful to A. Wachtel for valuable discussionsand suggestions. The present project was supported bythe National Research Found, Luxembourg, in the frameof project FNR/A11/02 and of the AFR PhD Grant 2014-2, No. 9114110.
Appendix A: Steady-state distributions
The generic expression for the steady-state distribu-tion (18) can be obtained as follows. The steady-stateequations can be expressed as¯ Z (cid:8) ¯ Z − ¯ Z (cid:9) + ¯ Z ¯ Z = 0 , ¯ Z (cid:8) ¯ Z k +1 − Z k + ¯ Z k − (cid:9) + ¯ Z (cid:8) ¯ Z k − ¯ Z k − (cid:9) = ¯ I k κ δ k k y ∈ Ω Y , for k ≥ . (A1)Defining the variable ∆ ¯ Z k ≡ ¯ Z k − ¯ Z k − they become¯ Z ∆ ¯ Z + ¯ Z ¯ Z = 0 , ¯ Z (cid:8) ∆ ¯ Z k +1 − ∆ ¯ Z k (cid:9) + ¯ Z ∆ ¯ Z k = ¯ I k κ δ k k y ∈ Ω Y , for k ≥ . (A2)Hence, by hierarchically substituting these expression oneinto the other, and using the variable ¯ r ≡ − ¯ Z / ¯ Z , weobtain∆ ¯ Z k = − (1 − ¯ r ) ¯ Z ¯ r k − ++ X k y ∈ Ω Y ¯ I k y κ ¯ r k − k y − Θ ( k − k y − , (A3)where Θ( · ) represents the discrete step function:Θ( k ) = ( k < , k ≥ . (A4)Finally,¯ Z k = k X i =1 ∆ ¯ Z i = ¯ Z (¯ r ) k − ++ X k y ∈ Ω Y ¯ I k y κ − (¯ r ) k − k y − ¯ r Θ ( k − k y − , (A5)which corresponds to the equation (18) in the main text. Appendix B: Three chemostats steady state
We discuss the uniqueness and stability conditions forthe steady state when three chemostats are present.From the constraints on the steady state (23) and fromthe condition for the external currents (24), we can writea single steady state condition involving just ¯ r as variable: (cid:0) ν k y3 ¯ Z k y1 + ν k y1 ¯ Z k y2 (cid:1) (¯ r ) ν k y1 + ν k y3 − (cid:0) ν k y1 + ν k y3 (cid:1) ¯ Z k y2 (¯ r ) ν k y1 − (cid:0) ν k y3 ¯ Z k y1 + ν k y1 ¯ Z k y3 (cid:1) (¯ r ) ν k y3 + (cid:0) ν k y3 ¯ Z k y2 + ν k y1 ¯ Z k y3 (cid:1) = 0 . (B1)1 y (a) x < y < y (b) x < y > y (c) x > y < y (d) x > y > FIG. 9: Plots of the hyperbola (dark purple curve) andpower law (light purple curve) in (B2) for different con-figurations of parameters. The center of the hyperbola ishighlighted by a dark purple dot, while the physical regionby the dashed orange lines.
Let us define the variables x ≡ (¯ r ) ν k y3 and y ≡ (¯ r ) ν k y1 ,so that the above-expressed steady-state condition can bewritten as the intersection of two curves: a rectangularhyperbola and a power law function y h = y − z x − x y p = ( x ) ν k y1 /ν k y3 , (B2)where the coefficients are given by x = (cid:0) ν k y1 + ν k y3 (cid:1) ¯ Z k y2 ν k y3 ¯ Z k y1 + ν k y1 ¯ Z k y2 ,y = ν k y3 ¯ Z k y1 + ν k y1 ¯ Z k y3 ν k y3 ¯ Z k y1 + ν k y1 ¯ Z k y2 ,z = ν k y1 ν k y3 (cid:0) ¯ Z k y2 − ¯ Z k y1 (cid:1) (cid:0) ¯ Z k y2 − ¯ Z k y3 (cid:1)(cid:0) ν k y3 ¯ Z k y1 + ν k y1 ¯ Z k y2 (cid:1) . (B3)[The subscripts h and p simply help us to distinguish thetwo functions.] From a geometrical point of view, physicalsolutions are represented by those intersection points lyingin ( x, y ) ∈ (0 , × (0 , Z k y1 > ¯ Z k y3 we observe that all of thepossible configurations of chemostat concentrations aredescribed by the following four cases for the parameters x and y . • x < y < z < Z k y1 > ¯ Z k y2 > ¯ Z k y3 .In this case we have always one and only one solu-tion. Indeed, the center of the hyperbola ( x , y )takes (0 , × (0 , x = 1 (which is non-physical). The left lower one,instead, always intersects the power law for valuesin (0 , × (0 ,
1) since y h ( x = 0) > • x < y > z > Z k y1 > ¯ Z k y2 and ¯ Z k y3 > ¯ Z k y2 .In this case we have one solution if and only if¯ Z k y1 > ¯ Z k y3 . The center of the hyperbola lies in(0 , × (1 , ∞ ) and the upper left branch of the hy-perbola never intersects the power law. The rightlower one, instead, always intersects the power lawin x = 1 , y = 1 (fig. 9b). We have a furtherintersection in the physical region if and only if d y p d x (cid:12)(cid:12)(cid:12) x =1 > d y p d x (cid:12)(cid:12)(cid:12) x =1 , which holds iff ¯ Z k y1 > ¯ Z k y3 —indeed, x ∗ : y h ( x ∗ ) = 0 is such that x ∗ >
0, for anychoice of the chemostats. • x > y < z > Z k y1 < ¯ Z k y2 and ¯ Z k y3 < ¯ Z k y2 .Once again, we have one solution if and only if¯ Z k y1 > ¯ Z k y3 . The center of the hyperbola liesin (1 , ∞ ) × (0 ,
1) and the right lower branch ofthe hyperbola never intersects the power law. Theupper left one, instead, always intersects the powerlaw in x = 1 , y = 1 (fig. 9c). We have a furtherintersection in the physical region if and only if d y p d x (cid:12)(cid:12)(cid:12) x =1 > d y p d x (cid:12)(cid:12)(cid:12) x =1 , which holds iff ¯ Z k y1 > ¯ Z k y3 —indeed, y h (0) > • x > y > z < Z k y1 < ¯ Z k y2 < ¯ Z k y3 .In this case we have no solutions. Indeed, the centerof the hyperbola lies in ( x, y ) ∈ (1 , ∞ ) × (1 , ∞ ) andneither the upper right nor the lower left branch ofthe hyperbola intersects the power law in the physi-cal region. The left lower one, indeed, always inter-sects the power law in (1 ,
1) which is non-physical(fig. 9c).Summarizing, we have a unique steady state wheneverthe concentration of the largest chemostat is higher thanthe concentration of the smallest one: ¯ Z k y1 > ¯ Z k y3 . Stability.
In order to prove the stability of the fixedpoint we resort to the following Lyapunov function: L = X k Z k ln Z k Z k s − ( Z − Z s ) . (B4)It is easy to prove that this function is always positiveand vanishes only for Z k = Z k s , where Z k s represents the2steady-state solution. If the steady-state solution exists,namely if exists Z k s : ˙ Z k s = 0, the time derivative of theLyapunov function (B4) can be written asd L d t = X k x ˙ Z k x ln Z k x Z k x s . (B5)Close to the steady state the above derivative is negative.Spanning the phase space with small perturbations onevery concentration, we always obtain d L d t ≤
0, where theequal sign is reached only at the steady state. Disregard-ing the infinite dimension of the phase space, we considerthe independent set of perturbations labeled with theindex k x and quantified by the small real value (cid:15)Z k x = Z k x s + (cid:15)δ k x k x , | (cid:15) | (cid:28) min k x Z k x s . (B6)Embedding these perturbation in (B5) and using the rateequations (5) we obtaind L d t ’ − κZ (cid:0) Z s − Z − Z (cid:1) (cid:15) , for k x = 1 , d L d t ’ − κZ k x s (cid:16) Z s + 2 Z k x s + − Z k x +1s − Z k x − − Z (cid:17) (cid:15) , for k x = 1 , (B7)which are always negative, no matter the sign of theperturbation. D. Nelson and M. Cox,
Lehninger Principles of Biochemistry (W.H. Freeman, 2008), pp. 244–246. S. G. Ball and M. K. Morell, Annu. Rev. Plant Biol. , 207(2003). Enzymes’ lack of specificity is nowadays getting increased at-tention. The most notorious examples are enzymes involved ingenetic processes and which are unspecific because they need toact on different substrates, e.g. polymerases and synthetases act-ing on the four nucleo-bases . More recently, enzymes involvedin metabolic processes have also been shown to lack specificity .These observations have led to various theoretical studies on thethermodynamical cost of error correction . G. Jones and W. Whelan, Carbohydr. Res. , 483 (1969). T. Takaha and S. M. Smith, Biotechnol. Genet. Eng. Rev. ,257 (1999). C. Colleoni, D. Dauvillée, G. Mouille, M. Morell, M. Samuel, M.-C. Slomiany, L. Liénard, F. Wattebled, C. d’Hulst, and S. Ball,Plant Physiol. , 1005 (1999). Ö. Kartal, S. Mahlow, A. Skupin, and O. Ebenhöh, Mol. Syst.Biol. , 542 (2011). R. N. Goldberg, D. Bell, Y. B. Tewari, and M. A. McLaughlin,Biophys. Chem. , 69 (1991). S. Lahiri, Y. Wang, M. Esposito, and D. Lacoste, New J. Phys. , 085008 (2015). See for instance: D. Kondepudi and I. Prigogine,
Modern Thermo-dynamics: From Heat Engines to Dissipative Structures (Wiley,2014), § 9.4. M. Polettini and M. Esposito, J. Chem. Phys. , 024117 (2014). P. Krapivsky, E. Ben-Naim, and S. Redner,
A Kinetic View ofStatistical Physics (Cambridge University Press, 2010), ch. 5. Indeed, the symmetry of the reaction scheme (1) is inherited bythe stoichiometric matrix, which satisfies ∇ knm = −∇ km − n +1 . We refer to eq. 5.99 in Ref. for the monomers-exchange processeswithout constraints. To be precise, chemostats represent reservoirs of particles towhom the system is connected. However, in order to simplify thenomenclature we refer to the chemical species exchanged by thereservoirs as chemostat, as well. S. Schuster and R. Schuster, J. Math. Chem. , 25 (1989). R. A. Alberty,
Thermodynamics of biochemical reactions (JohnWiley & Sons, 2003). B. O. Palsson,
Systems Biology: properties of reconstructed net-works. (Cambridge University Press, 2006). Z. Budrikis, G. Costantini, C. La Porta, and S. Zapperi, Nat.Commun. , 3620 (2014). T. P. J. Knowles, W. Shu, G. L. Devlin, S. Meehan, S. Auer,C. M. Dobson, and M. E. Welland, Proc. Natl. Acad. Sci. U.S.A. , 10016 (2007). T. P. J. Knowles, C. A. Waudby, G. L. Devlin, S. I. A. Cohen,A. Aguzzi, M. Vendruscolo, E. M. Terentjev, M. E. Welland, andC. M. Dobson, Science , 1533 (2009). S. I. A. Cohen, S. Linse, L. M. Luheshi, E. Hellstrand, D. A.White, L. Rajah, D. E. Otzen, M. Vendruscolo, C. M. Dobson,and T. P. J. Knowles, Proc. Natl. Acad. Sci. U.S.A. , 9758(2013). W. Bialek,
Biophysics: Searching for Principles (Princeton Uni-versity Press, 2012), § 4.5. C. L. Linster, E. van Schaftingen, and A. D. Hanson, Nat. Chem.Biol. , 72 (2013). J. J. Hopfield, Proc. Natl. Acad. Sci. U.S.A. , 4135 (1974). J. Ninio, Biochimie , 587 (1975). C. H. Bennett, Biosystems , 85 (1979). D. Andrieux and P. Gaspard, Proc. Natl. Acad. Sci. U.S.A. ,9516 (2008). G. Bel, B. Munsky, and I. Nemenman, Phys. Biol. , 016003(2010). A. Murugan, D. H. Huse, and S. Leibler, Proc. Natl. Acad. Sci.U.S.A. , 12034 (2012). P. Sartori and S. Pigolotti, Phys. Rev. Lett. , 188101 (2013). R. Rao and L. Peliti, J. Stat. Mech. Theor. Exp.2015