Kinetics and thermodynamics of reversible polymerization in closed systems

Motivated by a recent work on the metabolism of carbohydrates in bacteria, we study the kinetics and thermodynamics of two classic models for reversible polymerization, one preserving the total polymer concentration and the other one not. The chemical kinetics is described by rate equations following the mass-action law. We consider a closed system and nonequilibrium initial conditions and show that the system dynamically evolves towards equilibrium where detailed balance is satisfied. The entropy production during this process can be expressed as the time derivative of a Lyapunov function. When the solvent is not included in the description and the dynamics conserves the total concentration of polymer, the Lyapunov function can be expressed as a Kullback-Leibler divergence between the nonequilibrium and the equilibrium polymer length distribution. The same result holds true when the solvent is explicitly included in the description and the solution is assumed dilute, whether or not the total polymer concentration is conserved. Furthermore, in this case a consistent nonequilibrium thermodynamic formulation can be established and the out-of-equilibrium thermodynamic enthalpy, entropy and free energy can be identified. Such a framework is useful in complementing standard kinetics studies with the dynamical evolution of thermodynamic quantities during polymerization.


Introduction
The processes of aggregation and polymerization are ubiquitous in nature, for instance they are present in the polymerization of proteins, the coagulation of blood, or even in the formation of stars. They are often modeled using the classic coagulation equation derived by Smoluchowsky [1,2]. In the forties, Flory [3,4] developed his own approach for reactive polymers, with an emphasis on their thermodynamic properties and on their most likely (equilibrium) size distribution. At this time, only irreversible polymerization was considered. The first study of the kinetics of reversible polymerization combining aggregation and fragmentation processes was carried out by Tobolsky et al [5].
In the seventies, reversible polymerization became a central topic in studies on the association of amino acids into peptides and on the self-assembly of actin. The thermodynamics of assembly of these polymers forms the topic of the now classical treaty by Oosawa and Asakura [6]. At about the same time, Hill made many groundbreaking contributions to non-equilibrium statistical physics and thermodynamics, which allowed him to describe not only the self-assembly of biopolymers like actin and microtubules, but also their free energy transduction by biopolymers and complex chemical networks [7]. In the eighties, Cohen and Benedek [8] revisited the work of Flory and Stockmayer, by showing that the Flory polymer length distribution is obtained under the assumption of equal free energies of bond formation for all bonds of the same type. They also showed that the kinetically evolving polymer distribution does not have the Flory form in general, and they analyzed the irreversible kinetics of the sol-gel transition.

Reversible polymerization with one conservation law
We consider a reversible polymerization process made of the following elementary reactions Therefore in this case, there is only one conservation law, that of M, and one can assume M = 1. Note also that equation (2) generalizes the Becker-Döring equations which describe the dynamic evolution of clusters that gain or lose only one unit at a time [34,35] Assuming that the reaction (1) can be treated as elementary, the entropy production rate is given by [36,37] ( ) where R is the gas constant. The entropy production rate vanishes when the system reaches equilibrium, i.e. when and only when the detailed balance is satisfied: nm n m nm n m eq eq eq = + − + Using the inequality x x ln 1 ⩽ − , which holds for all x 0 > , one easily proves that the following quantity is non-negative, convex and vanishes only at equilibrium Indeed, by taking the time derivative of L(t) and using the definition of c, one obtains dL dt R cl l ∑ = c c log l l eq . Now using (2), we obtain two terms. The first term is By adding these two contributions and using the detailed balance condition of (4) as well as (3), one finds that This proves that L is a Lyapunov function, i.e. a non-negative, monotonically decreasing function, which vanishes at and only at equilibrium. The existence of such a function implies that the dynamics will always relax to a unique equilibrium state. In [38], the authors show that for a chemical system containing a finite number of homogeneous phases, a Gibbs free energy function exists that is minimum at equilibrium. In the context of reacting polymers, a similar free energy function has been derived in [9]. In view of the non-increasing property of this function, the authors have coined the term 'F-Theorem'. We discuss below its relation to the more usual H theorem.

One-fluid version
Until now, our model was exclusively expressed in terms of concentrations and could be used to study nonchemical dynamics such as population dynamics. We now introduce the one-fluid model, where the solvent is not described either by choice or because it is absent. The molar fractions of the polymers of length l, with l 1 ⩾ , are It is important to note that this Lyapunov function is in general distinct from the relative entropy or Kullback-Leibler divergence between the distribution x l and x l eq , which represents only the first term in equation (5). The reason for this difference is that c l , contrary to x l , cannot always be interpreted as a probability since its norm is not always conserved. The two quantities however become equivalent (i.e. they only differ by a constant) when the total concentration of polymers c(t) is constant in time. If furthermore − x x ln l l l eq ∑ is constant in time (the meaning of this assumption will become clear in the two-fluid model), then the negative of the Shannon entropy (which is related to the classic notion of free energy of mixing introduced in [3,4] as explained in [39])

Two-fluid version
In order to make contact with thermodynamics, we now introduce the two-fluid model which includes the solvent explicitly in the list of chemical species. Then, the molar fraction of the polymer of length l becomes y t c t C t ( ) ( ) ( ) l l = , which importantly is now defined with respect to the total concentration of all species including the (time-independent) solvent concentration c 0 (water for instance): C t c t ( ) ( ) If there is no solvent, c 0 0 = and one recovers the molar fraction of the previous section which was denoted by x l . In dilute solutions, since y 0 is very close to one and the other y l are much smaller, C(t) becomes almost constant: C c 0 ≈ . The chemical potential of a polymer of length l in a dilute solution is defined by RT y ln l l l is the standard reference chemical potential and h l 0 and s 0 l are the standard enthalpy and entropy, respectively. We restrict ourselves here to ideal solutions, and by this we assume that this form of chemical potential applies not only to the polymers (l 0 ≠ ) but also to the solvent. An interesting study of the effect of non-ideality on the time evolution thermodynamic quantities during reversible polymerization can be found in [40].
Let us define the intensive enthalpy function as and CS  = . In equation (10), the first term proportional to s l 0 therefore represents the entropic contribution due to the disorder in the internal degrees of freedom of each polymer, while the second term represents the nonequilibrium entropy in the distribution of the variables y l [26]. Let us introduce the intensive free enthalpy where we use the definition of the chemical potential in the last equality, and its extensive counterpart CG  = . Since the change of chemical potential associated with each reaction must vanish at equilibrium, i.e. Combining (12) with (4), we find that the kinetic constants must satisfy the local detailed balance where C eq denotes the equilibrium value taken by C. We note that since the chemical reactions do not involve a solvent, we formally define the rate constants with any zero subscript (n or m) to be zero. Naturally, equation (13) is not applicable in this case. The validity of equation (13) relies on two main assumptions: the first one is that of dilute solutions, while the second one is the ideality of the heat bath. The latter assumption means that there are no hidden degrees of freedom which can dissipate energy in the chemical reactions under consideration, which implies in particular that these reactions must be elementary. Using the definitions of the enthalpy and entropy functions, we find that The entropy production (3), using (13) and the chemical potential definition, may be written as Using equations (14) and (15), it can be rewritten as The first term is the heat flow, the second is the entropy change and the third and fourth terms represent a contribution due to the change in the total concentration. It is important to note that within the two-fluid model with ideal solutions, since C(t) is essentially constant, these latter two contributions are negligible. Neglecting these terms, the entropy production can be expressed as the change in free energy which is also equal to a change in Kulback-Leibler divergence between the nonequilibrium and the equilibrium polymer distribution We finally note that when all the polymerization reactions are neutral from a standard chemical potential standpoint i.e. when n m n m Since C can be assumed constant, the Lyapunov function can be expressed only in terms of the Shannon entropy constructed from y l instead of the full KL divergence. The dynamics can then be compared to a Boltzmann equation where the relaxation to equilibrium is purely driven by the maximization of the Shannon entropy, as in the H-Theorem.

Application to the String model
As a simple realization of the reversible polymerization given by equation (1), we now consider the String model which assumes constant rates of aggregation and fragmentation, independent of the length of the reacting polymers. Following [13], we choose k  The evolution equation (22) can be solved using an exponential ansatz of the form [13] c t a t a t ( ) (1 ( )) ( ), , which translates into the condition a (0) 0 = . Unfortunately, the exponential ansatz cannot be used to describe more general initial conditions, which cannot be accounted for by such a simple l-independent condition on a(0). For the monomer-only initial condition, the following explicit solution is obtained: while the entropy production rate given by equation (6) is For completeness, we also discuss an approach using generating functions to study the String model without resorting to the exponential ansatz which is restricted to the monomer-only initial condition. The full dynamics remains nevertheless complex to solve within this approach. Introducing the generating function z t c t z ( , ) ( ) l l l 1 ∑ χ = ⩾ and using equation (22), we obtain the dynamical equation The stationary solution of this differential equation, i.e. the solution of t 0 χ ∂ ∂ = has the following form: Besides the stationary solution, it is also possible to obtain analytically the evolution of the total concentration c(t). To show this, we take the limit z 1 → in equation (32). Using the l'Hospital rule, we obtain The explicit solution of equation (35) for c (0) 1 = , as imposed by our choice of initial conditions of the form c t ( 0) One can verify that equation (27) is recovered for the monomer only initial condition, M = 1, as expected. Furthermore, one recovers that c(t) tends towards c ( ) 2 eq Δ λ = − as t → ∞, which is the equilibrium concentration entering equations (33) and (34).
Incidentally, one may wonder whether this behavior of the total concentration agrees with the predictions of [6] regarding the notion of critical concentration in reversible polymerization. This is indeed the case: if one evaluates the concentration of monomers c 1 as a function of M with equation (34) and eliminating c eq using the above expression, one finds a function of M which first increases rapidly and then reaches a plateau for M λ ⩾ . Naturally, this is not a sharp transition but rather a cross-over between two regimes. One can also look at the average length of the polymer M c eq which increases significantly when M becomes larger than λ. Both features indicate that λ represents the critical concentration of this model [6].
As shown in the right part of figure 2, which has been obtained by explicit numerical integration, the polymer concentration c(t) either decreases as a function of time for the monomer-only initial condition (M = 1) or increases as a function of time when the initial condition corresponds to polymers of length 3 or above (M 3 ⩾ ), while it remains constant in the case of dimers (M = 2). Intuitively, when the initial condition is monomer-only (M = 1), there is mainly aggregation of monomers, so that the net concentration must decrease with time. On the other hand, if the initial solution consists of long polymers (M 2 > ), the probability of fragmentation is higher than that of aggregation. As a result, the total concentration must increase with time. For the dimer-only initial condition (M = 2), the probability of fragmentation is the same as that of aggregation, so that the net concentration stays constant. As a result of this time dependence of c(t), we also see in the left part of figure 2, that the Shannon entropy S Sh does not always increase monotonically as a function of time. It does so for M 2 ⩽ but not for M 2 > , where it presents an overshoot before reaching its equilibrium value. Such an overshoot reveals that the Shannon entropy is not a Lyapunov function L as discussed in the previous section.

Two-fluid version of the String model
One of the main differences between the two-fluid approach as compared with the previous one with a single fluid, is the existence of the local detailed balance condition, namely equation (13), which connects the rate constants to the difference of standard chemical potentials. Furthermore, the specific form of the standard chemical potentials enters in the equilibrium length distribution of the polymers and into the kinetics of the selfassembly process.
For instance, if the polymers self-assemble linearly, the standard chemical potential of a polymer of length l, . From equation (13), it follows that > , the entropy increases monotonically whereas we see that it decreases for M = 1. As in the case of the one-fluid model, this decrease is not inconsistent since the entropy is not a Lyapunov function. We note that the non-monotonicity that was present in the one-fluid model in figure 2 for M 1 > is absent in the two-fluid case.
If the monomers were to assemble in the polymer in a different way, for instance in the form of disks instead of linear chains, the standard chemical potentials would be different. In such a case, under similar assumptions to those above, these chemical potentials would be of the form l l RT ( ) , where α is again some constant characteristic of the monomer-monomer and monomer-solvent interaction [41]. The term in l represents the contribution of the surface energy of the cluster of size l. This term necessarily implies that the rate constants k nm + , k nm − must depend on n and m in order to satisfy the local detailed balance condition equation (13). For such a case, the above derivation of a simple exponential for the equilibrium distribution would no longer hold and both the equilibrium and the dynamics will be more complex.
After symmetrizing this sum, we recover that This result is equivalent to equation (6) in the presence of the additional conservation law c˙0 = . This system will therefore relax to a unique equilibrium state, where Σ vanishes.
We now turn to the two-fluid version of the model. As in section 2.
we can express it, as in (20), as To summarize, we recover exactly the same results as in section 2, provided we treat the total polymer concentration c as a constant.

Application to the kinetics of glucanotransferases DPE1 and DPE2
In this section, we consider the polymerization of glycans by two enzymes studied by Kartal et al [18], namely the glucanotransferases DPE1 and DPE2. We show how to construct dynamical models that are compatible with the equilibrium polymer length distributions that they found and we study the dynamics of the Shannon entropy for various initial conditions.

Kinetics of glucanotransferases DPE1
Let us assume that the initial condition is not purely made of monomers, since the solution of [18] becomes singular in that limit (see equation (4) on P3 of [18] when the parameter DP 1 in = for instance), and let us construct an appropriate dynamics, choosing for simplicity constant rates k nm κ = independent of n and m. Using equation (41), we have for l 2 ⩾ , ( )

Conclusion
In this paper, we have considered two classic models for reversible polymerization in closed systems following the mass-action law, one preserving the total polymer concentration and the other one not. In both cases, the entropy production can be written as the time derivative of a Lyapunov function which guarantees the relaxation of any initial condition to a unique equilibrium satisfying a detailed balance. As such, these models could also describe non-chemical systems undergoing aggregation-fragmentation dynamics. When considering the polymerization dynamics in dilute solutions, we have shown that a consistent nonequilibrium thermodynamics can be established for both models. We find that entropy production is the negative of the time derivative of the nonequilibrium free energy of the system, which is a Lyapunov function and takes the form of a Kullback-Leibler divergence between the nonequilibrium and the equilibrium distribution of polymer length. A related result was found for the cyclical work performed by chemical machines feeding on polymers in [42]. Similar relations expressing the entropy production as a Kullback-Leibler divergence between the nonequilibrium and equilibrium distributions have also been found or used in many studies on stochastic thermodynamics [43][44][45].
As an application of reversible polymerization models which do not preserve the total polymer concentration, we have studied the String model. In this model, the rates of aggregation and fragmentation are constants, which leads to an exponential equilibrium distribution of polymer length. At the one-fluid level, we have observed that the Shannon entropy is non-monotonic, which is allowed since it differs from the Lyapunov function. At the two-fluid level where there is proper nonequilibrium thermodynamics, no such nonmonotonicity arises.
As an application of reversible polymerization models preserving the total polymer concentration in addition to the total number of monomers, we have studied two specific examples named DPE1 or DPE2 after [18]. We have shown how to construct dynamics which converge over long times to the expected form and we