In modeling the interior of cells by simulating a reaction–diffusion master equation over a grid of compartments, one employs the assumption that the copy numbers of various chemical species are small, discrete quantities. We show that, in this case, textbook expressions for the change in Gibbs free energy accompanying a chemical reaction or diffusion between adjacent compartments are inaccurate. We derive exact expressions for these free energy changes for the case of discrete copy numbers and show how these expressions reduce to traditional expressions under a series of successive approximations leveraging the relative sizes of the stoichiometric coefficients and the copy numbers of the solutes and solvent. Numerical results are presented to corroborate the claim that if the copy numbers are treated as discrete quantities, then only these more accurate expressions lead to correct behavior. Thus, the newly derived expressions are critical for correctly computing entropy production in mesoscopic simulations based on the reaction–diffusion master equation formalism.

In recent years, the coarse-grained computational modeling of intracellular environments has enjoyed significant advances. An important paradigm shared by many such models is to treat the evolution of reacting chemical species’ copy numbers and spatial distributions by simulating a reaction–diffusion master equation (RDME).1 In this approach, the system volume is divided into compartments, each with local values of the copy numbers and chemical potentials of the different chemical species (Fig. 1). The RDME is a differential equation describing the evolution of the probability PN,t of observing the vector of copy numbers N={Ni,A}iL,AΩ of chemical species i in compartment A at time t, where L is the set of solute species and Ω is the set of compartments in the system. The RDME can be written schematically as

dPN,tdt=M^+D^PN,t,
(1)

where M^ and D^ represent operators describing chemical reactions and inter-compartment diffusion, respectively.2 An in-depth description of the RDME approach can be found in Refs. 3 and 4, where the forms of the operators are discussed. Rather than directly solving Eq. (1), one often simulates trajectories of the vector N obeying the stochastic dynamics encoded in the RDME using a variant of the Gillespie algorithm.5 The sizes of the compartments are commonly determined by the Kuramoto lengths, the mean free diffusional path for a species before it participates in a chemical reaction. Within each compartment, the spatial distributions of the reacting species are assumed to be homogeneous, allowing the use of mass-action kinetics with the compartment’s local values of the species’ concentrations to describe the stochastic chemical reaction propensities. Molecules can additionally jump between adjacent compartments in “diffusion events” (whose propensities also depend on the compartments’ local concentrations of species) to give rise to concentration gradients on the scale of the compartment length. This modeling approach is appropriate when the Kuramoto length is small compared to the system size (i.e., the assumption of homogeneity over the system volume fails), yet large compared to the intermolecular distance scale.1 Examples of simulation platforms based on such an approach include Virtual Cell,6,7 lattice microbes,8,9 MesoRD,10 MEDYAN (Mechanochemical Dynamics of Active Networks),11,12 and others.13–17 

FIG. 1.

An example of a cubic compartment grid used in the simulation of a RDME. Each compartment, labeled with letters A, B, …, has local values of the quantities NR, NG, referring to the copy numbers of the red and green molecules in the compartment, and of μR, μG, referring to the chemical potentials of those molecules. Molecules can react with each other within compartments (long dashed arrow), as well as hop between adjacent compartments, representing diffusion (short dashed arrow).

FIG. 1.

An example of a cubic compartment grid used in the simulation of a RDME. Each compartment, labeled with letters A, B, …, has local values of the quantities NR, NG, referring to the copy numbers of the red and green molecules in the compartment, and of μR, μG, referring to the chemical potentials of those molecules. Molecules can react with each other within compartments (long dashed arrow), as well as hop between adjacent compartments, representing diffusion (short dashed arrow).

Close modal

One important aspect of simulating non-equilibrium biological systems is the computation of thermodynamic forces that drive the observed flux on the network of chemical reactions.18–20 Determining these forces can allow for the quantification of entropy production in chemically reactive systems.21,22 In several research groups, measuring entropy production in biological active matter has been a recent goal.23,24 For instance, in recent work, we quantified the entropy production rates of self-organizing non-equilibrium actomyosin networks in MEDYAN using the expressions derived here as a first step.12 The ability to measure dissipation in active matter systems will allow us to test the applicability of different physical organizing principles relating the production of entropy to the likelihood of observing certain trajectories.25,26 For isothermal, isobaric, chemically reactive solutions, which include most mesoscopic biological systems, measuring the total entropy production amounts to determining the change in Gibbs free energy that accompanies chemical reactions and diffusion down concentration gradients.27–29 A ubiquitous textbook expression for the change in Gibbs free energy G accompanying a chemical reaction is

ΔG=kBTlogKeq+kBTlogQ,
(2)

where Keq is the equilibrium constant, Q is the reaction quotient (we give definitions of these quantities below), kB is Boltzmann’s constant, and T is the temperature.30–35 At equilibrium, Q=Keq1, and as a result, ΔG = 0. However, in this paper, we argue that Eq. (2) is a biased approximation to the exact value of ΔG accompanying a chemical reaction that holds when the copy numbers of the reacting molecules are large, such as on the order of Avogadro’s number. When the system is small, such as when copy numbers are on the order of 100 as is often the case in RDME simulations of intracellular environments, thermodynamic expressions such as Eq. (2) require corrections.36,37

As a simple motivating example, consider a mixture of an even total number of two chemical species, red and green, which inter-convert at equal rates. At equilibrium, the copy numbers of these molecules will be equal, and Q=Keq1=1. Now, if a reaction were to occur at equilibrium to produce one additional red molecule and to eliminate one green molecule. Then we would expect that the Gibbs free energy of the system had increased, since we have left equilibrium where the free energy attains its minimum. However, the prediction of Eq. (2) is that ΔG = 0 for this reaction. The assumption whose violation leads to Eq. (2) being incorrect is that the copy number of chemical species is a continuous quantity. When these variables are considered as discrete, a different expression for ΔG must be used to give correct behavior.

Similarly, for diffusion between adjacent compartments, a common expression for the change in Gibbs free energy accompanying the jump of a molecule i from compartment A with copy number Ni,A to compartment B, where its copy number is Ni,B, is

ΔG=kBTlogNi,BNi,A.
(3)

Imagine, however, we have a situation where Ni,A = Ni,B and a molecule jumps from compartment A to B. The Gibbs free energy should have increased since we have departed from the highest entropy distribution of the molecules over the two compartments; however, Eq. (3) will predict that ΔG = 0.

In this paper, we derive exact expressions for the change in Gibbs free energy accompanying chemical reactions within compartments and diffusion events between compartments, and we further show how these expressions relate to the familiar textbook formulas [Eqs. (2) and (3)] through a series of approximations. We also discuss the assumptions involved in defining the Gibbs free energy of a grid of homogeneously mixed compartments that can exchange energy and particles, such as that used in a simulation of a RDME. Finally, we present numerical simulations using MEDYAN to demonstrate the need to use these more exact expressions for ΔG in order to obtain sensible results when copy numbers are treated as discrete variables. Only these more exact expressions will give correct, unbiased behavior when measuring entropy production in mesoscopic in silico studies of biological non-equilibrium systems that rely on the RDME formalism.

Here, we make successive approximations to the formula for ΔG accompanying a chemical reaction, and our notation reflects the level of approximation in which certain quantities are being used: when appropriate, we subscript quantities with a parenthesized number, i.e., ΔG(0), where increasing numbers represent more approximate versions. The symbol     ̃ will indicate that the quantities of chemical species are being represented by mole fractions χi, rather than concentrations Ci. In this section, we treat the case that our system comprises a single closed compartment of a homogeneous dilute solution in which a chemical reaction has occurred, and we derive an expression for ΔG. In this system, the number of solvent molecules is fixed and the solute molecules participate in chemical reactions, causing their copy numbers to change. In Secs. II B and II C, we consider a system with multiple weakly interacting compartments (which may, in general, comprise an irregular grid rather than the regular Cartesian grid visualized in Fig. 1), each of which contains a homogeneous solution with local copy numbers of solvent and solutes and between which both solvent and solutes can diffuse. The nearly exact result [Eq. (21)] obtained in this section will also apply to those systems, as we argue below.

Before restricting to the case of a single closed compartment, we establish notation for properties of the chemical species in a compartment grid. The chemical potential for species i in compartment A can be expressed as depending either upon the mole fraction, χi,A, or upon the concentration, Ci,A, of that species in the compartment,

μi,A=μi0̃(T,p)+kBTlogχi,A=μi0(T,p)+kBTlogCi,A,
(4)

where kB is Boltzmann’s constant, μi0̃(T,p) is the standard state chemical potential at temperature T and pressure p when working with the dimensionless χi,A, and μi0(T,p) is the same when working with Ci,A. Ci,A and χi,A both play the role of the “composition variable” leading to these different, yet equivalent expressions for the chemical potential.38 We make the distinction between dependence upon copy number and dependence upon concentration in order to establish parameters that can be used in simulation, which commonly works with copy numbers, based on those given in the literature, which typically use units of concentration. Here, we make the assumption of an ideal-dilute solution and neglect the coefficient of activity of the solute species.39 The mole fraction can be written as

χi,A=Ni,ANA=Ni,AjMNj,A,
(5)

where NA is the total copy number of molecules in compartment A and M is the set of all species including the solvent in the system. Similarly, the concentration can be written as

Ci,A=Ni,AΘA,
(6)

where

ΘA=NAvVA
(7)

is a conversion factor, NAv is Avogadro’s number, and VA is the compartment volume (we assume constant pressure and that the fluctuations in volume are negligible for these liquid systems, allowing the use of Gibbs free energy). Using Eqs. (4)–(6), we can write40 

μi0̃(T,p)=μi0(T,p)+kBTlogNAΘA.
(8)

Using standard arguments concerning thermodynamic extensivity, it is possible to establish that the Gibbs free energy of the solution in compartment A can be written as a weighted sum over the chemical potentials of the species,

GA(NA)=iMNi,Aμi,A(χi,A),
(9)

where NA={Ni,A}iM is the vector of species copy numbers Ni,A, and where we have explicitly indicated the dependency of μi,A upon mole fraction χi,A via Eq. (5).28,41 One may be concerned that Eq. (9) fails to apply when the copy numbers of solutes are small. As explained below, we make the assumption that boundary effects are still negligible, which allows us to treat G as a first order homogeneous function of the number of copies of the system.37 This is the necessary property to establish Eq. (9), so the small copy numbers of solutes do not render this approach invalid. We rely on Eqs. (4) and (9) to derive changes in Gibbs free energy accompanying chemical reactions and inter-compartment diffusion.

Consider a reaction of the general form

ν1X1+ν2X2+kk+υ1Y1+υ2Y2+,
(10)

where Xi represent reactants, Yj represent products, νi and υj are stoichiometric coefficients, and the rate of reaction is k+ to the right and k to the left. We have dropped the subscript A indicating the compartment in which the reaction takes place and now restrict to the case that our system is a single compartment. When this reaction has occurred once to the right, the copy numbers of reactants have changed NiNiνi and those of the products have changed NjNj + υj. We calculate the change in Gibbs free energy accompanying this reaction by considering it as resulting from these finite, discrete changes in copy numbers,37 not from infinitesimal changes. Using Eqs. (4), (5), and (9), the Gibbs free energy before the reaction has occurred can be written as

Ginitial=iRNiμi0̃+kBTlogNiN+jPNjμj0̃+kBTlogNjN+Nsμs0*̃+kBTlogNsN,
(11)

where R is the set of reactants, P is the set of products, the subscript s refers to the solvent, and where we have dropped the dependence of the standard state chemical potential on T and p. μi0̃ and μj0̃ describe the chemical potential at a reference concentration of the solute in the solvent (also referred to as the solute standard state), whereas μs0*̃ describes the chemical potential at a reference state of pure solvent (also referred to as the solvent standard state).39 We assume here, for simplicity and without loss of generality, that there are no solute species that have not participated in the reaction (i.e., spectator solute species). These species would also contribute terms to Eq. (11), but when we subtract the initial from the final Gibbs free energy, their inclusion would not yield a different result. The final Gibbs free energy is

Gfinal=iR(Niνi)μi0̃+kBTlogNiνiN+σ+jP(Nj+υj)μj0̃+kBTlogNj+υjN+σ+Nsμs0*̃+kBTlogNsN+σ,
(12)

where

σ=jPυjiRνi,
(13)

which is the amount by which the total species copy number N has changed. As described in Refs. 38 and 42, when σ ≠ 0, it is important to account for the solvent species in the calculation of free energy differences. This is because the free energy of the solvent, which is the most abundant species in the reaction volume, will not be the same after the reaction has taken place since its mole fraction will change as N changes to N + σ. Neglecting the solvent species when σ ≠ 0 leads to expressions for ΔG that are off by an amount σkBT.38 Whereas the authors of Refs. 38 and 42 describe the appearance of this erroneous term while formulating the Gibbs free energy as a function of a continuous degree of advancement of reaction = −dNi/νi = dNj/υj, here we treat the extent of reaction as a discrete quantity. In the limit that νi/Ni → 0, υj/Nj → 0 for each reactant and product, the discrete case passes into the continuum case; however, under the assumption of small copy numbers, we do not satisfy this limit. In  Appendix A, we discuss further differences between the continuum treatment and the discrete treatment, as well as the relation to the Gibbs–Duhem equation.

After some algebra (using the fact that ∑iRNi + ∑jPNj + Ns = N), we can write the change in Gibbs free energy as

ΔG(0)=GfinalGinitial=ΔG0̃+kBTlogiR(Niνi)NiνiNiNi×jP(Nj+υj)Nj+υjNjNjNN(N+σ)N+σ,
(14)

where

ΔG0̃=jPυjμj0̃jRνiμi0̃.
(15)

Equation (14) is exact, but we would like to avoid specifying N in simulation since the solvent is typically not modeled explicitly, as we elaborate on in Sec. II B. We would also like to determine ΔG0̃ from literature values of ΔG0. To these ends, we first rewrite Eq. (14) as

ΔG(0)=ΔG0̃+kBTlogQ̃(0)+kBTlogNN(N+σ)N+σ,
(16)

where

Q̃(0)=iR(Niνi)NiνiNiNijP(Nj+υj)Nj+υjNjNj.
(17)

From Eqs. (8) and (15), we can write

ΔG0̃=ΔG0+σkBTlogNΘ,
(18)

where

ΔG0=jPυjμj0jRνiμi0.
(19)

Inserting this into Eq. (16), we have

ΔG(0)=ΔG0+kBTlogNΘσNN(N+σ)N+σ+kBTlogQ̃(0)=ΔG0σkBTlogΘ+kBTlogNN+σN+σ+kBTlogQ̃(0).
(20)

We now make the approximation that Nσ, which is certainly reasonable for most realistic parameterizations of the compartment grid (in the example of a 0.125 µm3 compartment filled with water, N ∼ 109 while σ ∼ 1). With this, we can write the third term in Eq. (20) as −σkBT,43 giving

ΔG(1)=ΔG0σkBTlogΘσkBT+kBTlogQ̃(1),
(21)

where Q̃(1)=Q̃(0) (we updated the subscript to indicate the use of this quantity in a more approximate version of the formula for ΔG). We recommend using Eq. (21) in simulation because it allows us to incorporate literature values for ΔG0 and avoids specification of N. To understand the term −σkBT in Eq. (21) and to understand the relationship between Eqs. (21) and (16) and the textbook expression for ΔG, we proceed by making further approximations leveraging the large sizes of the solute copy numbers compared to their stoichiometric coefficients. First, we rewrite Q̃(1) as

Q̃(1)=iR1νiNiNi(Niνi)νijP1+υjNjNj(Nj+υj)υj.
(22)

Assuming Niνi and Njυj, and using the limit

limx1±yxx=e±y,
(23)

we obtain

kBTlogQ̃(1)σkBT+kBTlogiR(Niνi)νijP(Nj+υj)υj.
(24)

Inserting this into Eq. (21), canceling the term σkBT, gives

ΔG(2)=ΔG0σkBTlogΘ+kBTlogQ̃(2),
(25)

where

Q̃(2)=iR(Niνi)νijP(Nj+υj)υj.
(26)

Finally, introducing

Q̃(3)=iRNiνijPNjυj
(27)

and discarding terms in ΔG(2) that scale like νiNi or vjNj, we arrive at

ΔG(3)=ΔG0σkBTlogΘ+kBTlogQ̃(3)=ΔG0+kBTlogQ,
(28)

where

Q=iRCiνijPCjυj.
(29)

We see that Eq. (28) is obtained from Eq. (25) upon making the approximations NiνiNi and Nj + υjNj. This final result in Eq. (28) is a standard textbook expression for the change in Gibbs free energy.30–35 

We, thus, have four expressions for ΔG of chemical reactions:

  • Equation (14) for ΔG(0) is exact; however, it requires specifying N.

  • Equation (21) for ΔG(1) uses the approximation Nσ. We recommend the use of this expression because it is the most exact expression for which we need not specify N, and since it is written in terms of ΔG0, for which literature values can be found.

  • Equation (25) for ΔG(2) uses the approximations Niνi and Njυj.

  • Equation (28) for ΔG(3) uses the approximations Niνi and Njυj again.

In  Appendix B, we provide expressions for the accuracy of these approximations. Note that, without specifying the copy number of solvents Ns, we can only approximately compute changes in Gibbs free energy, not the instantaneous Gibbs free energy of the system. Typically, only the changes are of interest.

The definition of Gibbs free energy states that

ΔG=ΔHTΔS.
(30)

Comparing this expression to Eq. (28), we identify

ΔG0=ΔHTΔS0
(31)

and

TΔSmix=kBTlogQ,
(32)

where ΔH is the enthalpy of reaction, and ΔS0 and ΔSmix represent the changes in entropy due to the molecular conformations and the translational motions, respectively,

ΔS=ΔS0+ΔSmix.
(33)

It is instructive to realize that the discrepancies between the various expressions for ΔG are due to how the ΔSmix term in Eq. (33) is written.

When employing any of the above expressions that involve the logarithms of products of copy numbers that are raised to the power of other copy numbers, we recommend splitting the logarithm of products into a sum of logarithms as well as using log xy = y log x in order to prevent overflow resulting from computing very large numbers. The results of this section do not assume a system consisting of multiple, weakly interacting compartments, which we discuss next.

Now that we have treated the scenario of a reaction event occurring within a single compartment, we want to generalize to the case of a grid of compartments. For this, we develop an argument based on timescales that will allow us to track the copy numbers of the reactive solutes while ignoring those of the inert solvent. This is a necessary modeling feature due to the computational infeasibility of tracking the solvent copy numbers in each compartment for long times. After describing the thermodynamic framework for a grid of compartments, in Sec. II C, we derive the change in G resulting from diffusion of solutes between adjacent compartments.

In simulating a RDME, one commonly treats diffusion between adjacent compartments and chemical reactions within compartments using an augmented set of all species and reactions in the system that treats species as distinct if they belong to separate compartments. Thus, if there are |L| reacting species and |Ω| compartments, where L and Ω are the sets of solute species and compartments, respectively, then in the augmented set, there are |L∥Ω| species tracked. The number of reactions in the augmented system, including r chemical reactions per compartment and roughly z|L| diffusion events per compartments (where z is the assumed constant number of neighbors of each compartment, ignoring boundary compartments), is |Ω|(r + z|L|). This augmented set of species and reactions is then simulated using the Gillespie algorithm in which the reaction propensities are appropriately scaled according to the compartment volumes.4 

Crucial to the justification of this strategy to simulate a RDME is the assumption that, within each compartment, the reacting species can be considered homogeneously distributed so that one may use mass-action kinetics to determine the propensities. This assumption amounts to the condition that the timescale describing diffusion within compartments, τD, is much less than the timescale of chemical reactions, τC,

τDτC.
(34)

This comparison should be done for each diffusing and reacting species, and the timescale of the fastest reaction (taken as the inverse of the propensity) for each species should be used. If the condition holds, then the process of intra-compartmental diffusion will homogenize the solution faster than chemical reactions that occur, so the assumption of mass-action kinetics holds. Let the dimension of the space be d, the length of the (here assumed cubical) compartments be h, and the diffusion constant of a species be D. Then,

τDh22dD.
(35)

The Kuramoto length is given by

lK=2dDτC,
(36)

so one can see that the condition τDτC is equivalent to the condition lKh, and thus, one can enforce this condition by choosing a smaller compartment size h. For fixed total volume, there is a trade-off between h and |Ω|, which determines the size of the augmented system and therefore the computational efficiency. The timescale of intra-compartment diffusion τD is approximately equal to the timescale of inter-compartment diffusion (which can be given as the inverse of kD=Dh2), so a third way of describing this condition is that the frequency of jumps between adjacent compartments is much greater than the frequency of chemical reactions inside the compartments, which can be checked empirically in simulation.1,4 We note that, in the literature, these expressions may differ up to a constant coefficient depending on the reference.

In order to approximately describe the thermodynamics of this system, we distinguish between the inert solvent and the dilute, chemically reactive solutes. We assume here that the system is impermeable to the flow of either kind of species to the exterior. Within the system, all species are permeable (though local diffusion constants may be incorporated for each species1). We treat the solutes explicitly (at the level of their compartment copy numbers), whereas we model the solvent implicitly through an appropriate limit as in the steps leading to Eq. (21). To this end, we assume the existence of a laboratory timescale τl whose purpose is to define the temporal resolution of our measurements such that any changes occurring on a timescale longer than τl will be measured. We assume that τl is much longer than the timescale describing the local compartment fluctuations in the solvent copy numbers τs, yet shorter than but on the order of the timescale describing diffusion of the solutes τD. On this timescale τl, then in the time between the chemical reactions and diffusion events involving the solutes, the system is quasi-equilibrated with respect to fluctuations involving the fast process of solvent diffusion, and thus, the system may be assigned well-defined values of Gibbs free energy.44 This hierarchy of timescales can be written as

τsτlτDτC.
(37)

Typical ratios of the diffusion constants for solute to solvent are in the range of 1/10 to 1/100, placing τD/τs in the range of 10–100.45 

On the timescale τl, there is enough temporal resolution to track the diffusion and chemical reactions of the solutes while allowing averaging over the fluctuations of the solvent. To describe activity occurring over the large grid of compartments for extended systems, we introduce new timescales τxg, where x refers to any of the timescales defined above. If we hold the compartment size h and the chemical concentrations fixed and add more compartments to our system, the rates of solute diffusion events and chemical reactions occurring anywhere in the system will scale as |Ω|, the number of compartments, and thus, the timescales needed to describe them scale as |Ω|−1, i.e., τlgτl/|Ω|. One might be concerned that τlg will be less than τs for large systems, ostensibly violating our requirement that we will be able to average over the solvent fluctuations to define quasi-equilibrated states. However, the timescale of solvent fluctuations across the whole grid, τsg, will also scale inversely with |Ω|, so the condition for being able to average over solvent fluctuations expressed in Eq. (37) does not ultimately depend on the number of compartments. In other words, for systems with many compartments, as long as Eq. (37) holds for one compartment, we can be sure that our laboratory timescale that describes the whole grid, τlg, will be short enough to describe processes involving the solutes, while long enough to allow averaging over the fluctuations of the solvent occurring locally in each compartment.

We assume that the exterior of the system acts as a reservoir for the thermodynamic variables p and T. The volume V of the system also remains constant; however, under the assumption of the solvent being an incompressible liquid, the change in quantity pV is approximately zero, and for each reaction, the change in Gibbs free energy is equal to that of the Helmholtz free energy. Thus, it is inconsequential whether we consider p or V to be reservoir variables, and we choose p in order to speak of the Gibbs free energy of the system. We can write the Gibbs free energy of the system as G(N, p, T), where N={{Ni,A}iM}AΩ represents the set of copy numbers of all solute and solvent species in each compartment in the grid. We further assume the compartments to be only weakly interacting, that is, they can exchange energy and particles, but the interaction of the two subsystems does not contribute a term to the Gibbs free energy of the system. This is equivalent to assuming that the Gibbs free energy of the compartments GA is linearly additive,

G=AΩGA,
(38)

without any terms of the form GAB. To justify this, we first note that the interaction free energy between two adjacent compartments will primarily be due to the interaction of the solvent at the interface. This interfacial free energy will, even for mesoscopically sized compartments, be small compared to the free energy due to the bulk of the compartment. Furthermore, our main interest will be in changes in the Gibbs free energy of the system due to solute diffusion between compartments and chemical reactions, neither of which will significantly affect the interfacial free energy, so all terms of the form GAB will drop out of the expression for ΔG.

In our approach of averaging over fluctuations in the solvent amounts and taking the limit that this average is large compared to the copy numbers of the solutes, we are choosing to neglect the changes in Gibbs free energy of the system owing to the solvent fluctuations. We justify this with an argument that these fluctuations are small compared to those resulting from the activity of the solute molecules; it also arises out of necessity due to the computational expense of tracking the solvent fluctuations. On the laboratory timescale τl, each compartment has an average copy number of solvent molecules, Ns,A¯. The fluctuations in this quantity will have a standard deviation on the order of Ns,A¯1/2.44 As indicated previously [to arrive at Eq. (21)], we avoid specifying Ns,A¯ by taking the limit that it is much larger than the number of solute species, and thus, fluctuations in this quantity do not matter when computing nearly exact changes in the Gibbs free energy accompanying reactions involving the solutes. We need to establish that the change in Gibbs free energy of a compartment due to a fluctuation in the solvent copy number on the order of Ns,A¯1/2 is small compared to the change in Gibbs free energy accompanying chemical reactions and inter-compartment solute diffusion events and, thus, that the fluctuations in this quantity are not outweighing the changes that we are measuring. In  Appendix C, we show that if the concentrations of solutes are very different in the two compartments, then this Gibbs free energy change is on the order of εNs,A¯1/2kBT, where ε is the ratio of solutes to solvent. We show that this quantity is typically much less in magnitude than the change in Gibbs free energy accompanying a solute diffusion event. If the concentrations are nearly equal, then the Gibbs free energy is on the order of εkBT and is thus negligibly small. We make the modeling choice to ignore these changes in the Gibbs free energy resulting from solvent fluctuations because they tend to be small and they do not represent the processes we are interested in, which involve the activity of the solute molecules.

To summarize, we track the changes in the Gibbs free energy of the system by computing a value of ΔG whenever a chemical reaction occurs within a compartment or a solute diffusion event occurs between adjacent compartments. The timescales on which these events occur are much slower than that of solvent equilibration, allowing for well-defined values of G between the events. By the linear additivity of the compartments’ free energies, any change in the free energy of a single compartment is equal to the change in free energy of the whole system (i.e., ΔGA = ΔG). We assume that the diffusion of the solvent only contributes small, fluctuating, unbiased changes to the free energy that we ignore; we also assume that the amount of solvent is so large that one can take a limit and neglect the fluctuations in this quantity when computing the changes in free energy for processes involving the solute copy numbers.

To describe the change in Gibbs free energy accompanying a solute diffusion event between neighboring compartments, we use a similar approach to the one used above for chemical reactions. The key difference here is that, as opposed to multiple species being involved in a reaction taking place in a single compartment, we now have a single species involved in a diffusion event taking place between two compartments.

Consider species i diffusing from compartment A, where its initial copy number is Ni,A, to compartment B, where its initial copy number is Ni,B. The total copy numbers of molecules in compartments A and B are NA and NB, respectively. Assume there is just one spectator species constituting the solvent, labeled s with copy numbers Ns,A and Ns,B (we remove our uncertainty in the exact values of these numbers by taking a limit later). As a result of the diffusion event, we have the following changes in these quantities: Ni,ANi,A − 1, Ni,BNi,B + 1, Ns,ANs,A, Ns,BNs,B, NANA − 1, and NBNB + 1. The initial Gibbs free energy is

Ginitial=Ni,Aμi0̃+kBTlogNi,ANA+Ni,Bμi0̃+kBTlogNi,BNB+Ns,Aμs0,*̃+kBTlogNs,ANA+Ns,Bμs0,*̃+kBTlogNs,BNB,
(39)

and the final Gibbs free energy is

Gfinal=Ni,A1μi0̃+kBTlogNi,A1NA1+Ni,B+1μi0̃+kBTlogNi,B+1NB+1+Ns,Aμs0,*̃+kBTlogNs,ANA1+Ns,Bμs0,*̃+kBTlogNs,BNB+1.
(40)

The difference of these two expressions leads to an exact formula, similar to Eq. (16), which we omit here. Analogously to the approximation Nσ made in the context of chemical reactions, here we assume that NA, NB ≫ 1, which amounts to setting NA − 1 ≈ NA and NB + 1 ≈ NB in Eqs. (39) and (40). Using this approximation, we can express the change in Gibbs free energy as

ΔG=kBTlog(Ni,A1)(Ni,A1)Ni,ANi,A(Ni,B+1)(Ni,B+1)Ni,BNi,B.
(41)

We see that, analogously to making the approximation Niνi above, if we here take Ni,A, Ni,B ≫ 1, then Eq. (41) reduces to the common expression

ΔG=kBTlogNi,BNi,A.
(42)

The right-hand side of Eq. (42) will always be less than that of Eq. (41), which, although the difference is typically very slight, can lead to systematically biased calculations, as we describe in Sec. III.

The considerations leading to the main results of this paper [Eqs. (14), (21), and (41)] do not depend essentially on the assumption that the compartments comprising the simulation volume form a regular Cartesian grid (as illustrated in Fig. 1). A variety of strategies have been developed to discretize the simulation volume in a more sophisticated manner. These include unstructured meshes that use complex polygonal compartments useful for modeling curved surfaces, as well as adaptive meshes that use compartments with time-dependent boundaries that can increase computational efficiency.46–48 Hybrid, multiscale methods also exist, which combine a mesoscopic compartment-based description of part of the simulation volume with either a macroscopic or microscopic description elsewhere in the volume.49,50 Determining reaction and diffusion propensities in these frameworks is a somewhat complicated issue that we do not address here. Without undertaking a detailed analysis of how the above derivations of expressions for ΔG should be extended to apply in each of these computational frameworks, we observe that, given a method for calculating propensities leading to reaction and diffusion events, the expressions derived here can be applied directly to the case of an irregular grid of compartments, which encompasses several applications of interest. In an irregular grid, each compartment A has an arbitrary shape and volume VA, with VAVB, in general, for any pair A, B. It is important to ensure that for each of these compartments the assumptions regarding timescales described above still hold. If they do hold, then to apply the above results, it is only necessary to use compartment A’s local value of ΘA = NAvVA in place of Θ where it appears in the expressions for ΔG(1), ΔG(2), and ΔG(3) [Eqs. (21), (25), and (28), respectively].

Here, we perform stochastic simulations to illustrate the effects of approximating ΔG for reaction and diffusion events when copy numbers are considered small and discrete. We use MEDYAN, a simulation platform designed to study active networks at high resolution, which is equipped with a RDME simulation engine as described above.11 With this, we report on two different simulation setups to test the discrepancy between the nearly exact and approximate formulas for ΔG corresponding to reactions and to diffusion. We compare the nearly exact equations (21) and (41) for reactions and diffusion, respectively, with their approximate counterparts [Eqs. (28) and (42)]. We observe that only the nearly exact expressions result in sensible behavior, i.e., the rate of change of Gibbs free energy on average is 0 kBT/s when the system is at equilibrium.

To test the effect of approximation for ΔG of reactions, we consider the simple reaction scheme

A+Bkk+C,
(43)

where the rate constants to the right and left are k+ = 0.05 µM−1 s−1 and k = 0.01 s−1, respectively, giving an equilibrium constant of Keq = k/k+ = 0.2 µM. We perform stochastic simulations with the Next Reaction Method in MEDYAN51 using a single compartment of size 0.125 µm3 (i.e., in this example, there is no diffusion). To employ Eq. (21) in simulation, we first compute ΔG0σkBT log Θ − σkBT = 3.71 kBT for the forward reaction and −3.71 kBT for the reverse reaction (using ΔG0 = kBT log Keq), and then, when each reaction fires during simulation, the quantity Q̃(1) is computed from the instantaneous values of the copy numbers to determine ΔG(1). A similar approach is taken to employ Eq. (28). We begin with NA = 100, NB = NC = 50 and repeat a simulation of 100 s duration 3000 times to obtain averages of the trajectory of rates tΔG(t) resulting from the forward and reverse reactions [the notation tΔG(t) implies that the measured rates vary smoothly as a function of time; however, this quantity is calculated as the total change in Gibbs free energy resulting from discretely timed chemical events during 1 s-long windows]. Figure 2 displays the results of these simulations. Note how, while the two trajectories bear close similarity, the equilibrium value of tΔG(t) centers around 0 kBT/s for the nearly exact formulation [Eq. (21)], yet erroneously centers around ∼ −0.08 kBT/s for the approximate version [Eq. (28)].

FIG. 2.

Numerical results illustrating the difference between nearly exact and approximate formulations of ΔG for reactions and diffusion. (a) Averages over 3000 simulations of the chemical scheme A + BC, using the nearly exact Eq. (21). The blue curve represents the trajectory of tΔG resulting from the forward reaction, the red curve represents that from the reverse reaction, and the black curve represents their sum. Shaded regions represent the standard deviation over the 3000 repeated trials. The inset displays a blow-up of the black curve once the system has reached close to equilibrium. (b) The same as just described, but using the approximate Eq. (28). (c) A single trajectory of diffusion of 1000 molecules over a 1 µm3 cubic grid of eight compartments, beginning from a random initial spatial distribution. Values of ΔG are calculated using Eq. (41). (d) The same as just described; however, the values of ΔG are calculated using Eq. (42).

FIG. 2.

Numerical results illustrating the difference between nearly exact and approximate formulations of ΔG for reactions and diffusion. (a) Averages over 3000 simulations of the chemical scheme A + BC, using the nearly exact Eq. (21). The blue curve represents the trajectory of tΔG resulting from the forward reaction, the red curve represents that from the reverse reaction, and the black curve represents their sum. Shaded regions represent the standard deviation over the 3000 repeated trials. The inset displays a blow-up of the black curve once the system has reached close to equilibrium. (b) The same as just described, but using the approximate Eq. (28). (c) A single trajectory of diffusion of 1000 molecules over a 1 µm3 cubic grid of eight compartments, beginning from a random initial spatial distribution. Values of ΔG are calculated using Eq. (41). (d) The same as just described; however, the values of ΔG are calculated using Eq. (42).

Close modal

To test the effect of approximating ΔG for diffusion in a compartment-based reaction–diffusion scheme, we next employed MEDYAN to simulate diffusion of 1000 solute molecules with diffusion constant 20 µm2 s−1 in a 2 × 2 × 2 grid of compartments, each a cube with volume 0.125 µm3. The initial distribution of molecules is uniformly random over the compartment grid, and thus, the system begins near equilibrium and is then allowed to stochastically evolve for 100 s, i.e., the molecules hop randomly between adjacent compartments. For each diffusion event, the value of ΔG is determined using the nearly exact equation (28) for one run, and in another run, the approximate equation (42) is used. The difference in the trajectories of tΔG for these simulations is stark, as displayed in Fig. 2. While the trajectory centers around 0 kBT/s for the nearly exact formulation of ΔG, it erroneously centers around ∼ −1925 kBT/s for the approximate version. Diffusion events are very frequent in this system, occurring around 240 000 times/s (this number can be calculated from the parameters of the system and is also observed during simulations). Thus, since Eq. (42) is always less than the nearly exact quantity, even by a small amount on the order 0.05 kBT for this system, this systematic bias is amplified by the frequency of diffusion events to produce significant differences from the expected behavior, necessitating the use of a more exact formula for ΔG. Finally, we performed a simulation involving both reactions and diffusion across multiple compartments. We found again that only when the nearly exact formulas were used did the rate of change of Gibbs free energy center on 0 kBT/s at equilibrium. It is not additionally illuminating to show the data, so we do not display it here.

We have argued that when the copy numbers of the reactants and products are treated as small, discrete quantities, then certain approximations leading to the textbook formulas for ΔG of reaction and diffusion [Eqs. (2) and (3)] break down and lead to biased results. We emphasize that this is true only when the copy numbers and reaction occurrences are treated as discrete; when they are treated as continuous, one should use the textbook formulas. This can be shown by considering a continuous version of the chemical system described by Eq. (43). The time evolution of the concentrations of the chemical species is obtained by solving a system of ordinary differential equations that employ mass-action kinetics, and from this solution, the rates tΔG(t) resulting from the forward and reverse reactions are computed using the instantaneous values of the species’ copy numbers that enter into Eqs. (21) and (28). Figure 3 displays the results of these calculations. Here, the total rate of tΔG(t) only approaches 0 kBT/s, as it must at equilibrium, when Eq. (28) is used.

FIG. 3.

Analytical results illustrating that, when copy numbers are treated as continuous variables, the textbook formula for ΔG of reaction should be used. (a) Calculation of the rates tΔG(t) using Eq. (21) resulting from the forward (blue curve) and reverse (red curve) reactions, as well as their sum (black curve), as described in the main text. The inset shows a blow-up of these curves as the system approaches equilibrium. (b) The same as just described, but using Eq. (28).

FIG. 3.

Analytical results illustrating that, when copy numbers are treated as continuous variables, the textbook formula for ΔG of reaction should be used. (a) Calculation of the rates tΔG(t) using Eq. (21) resulting from the forward (blue curve) and reverse (red curve) reactions, as well as their sum (black curve), as described in the main text. The inset shows a blow-up of these curves as the system approaches equilibrium. (b) The same as just described, but using Eq. (28).

Close modal

When the copy numbers of the chemical species are treated as continuous, there is no notion of a single occurrence of a reaction; instead, the evolution of the system is parameterized by the continuous variable ξ that quantifies the extent of advancement of the reaction.52 In this framework, which is adopted in classical thermodynamics, the copy number of any species never jumps instantaneously from Ni to Ni + νi, and thus, the premise of the derivation presented above leading to Eq. (16) falls apart. This explains why the seemingly more accurate equation (21) is wrong when applied to chemical dynamics that are described using continuous variables to represent the copy number of species. When the chemical dynamics are modeled this way, the textbook expression [Eq. (28)] for the change in Gibbs free energy is valid. Another way to understand the difference in these expressions is to view Eq. (28) as giving the instantaneous slope of G with respect to the degree of advancement of the reaction when these quantities are viewed as continuous variables,38,42 whereas Eq. (16) effectively integrates the continuously varying quantity dG through a single discrete reaction occurrence.

If the copy numbers of chemical species are not large compared to the stoichiometric coefficients (i.e., Ni is not much greater than νi), it is best to describe the chemical dynamics by treating the copy numbers as discrete variables participating in stochastically timed chemical reactions. This is the philosophy adopted by several recent models of intracellular environments, where copy numbers of molecules of interest are sometimes quite small. In these cases, adoption of the more exact expression for ΔG can not only improve precision but also ensure correct behavior. If concentration gradients are expected to be a strong source of entropy production (e.g., diffusion of monomeric actin along the lengths of filopodia, see Refs. 53–55), then using the nearly exact formulas presented here can ensure that the resulting measurements of dissipation are not strongly biased.

The authors acknowledge A. Chandresekaran, H. Ni, and Q. Ni for improvement of this manuscript and helpful discussions. This work was supported by the National Science Foundation (Grant Nos. NSF 1632976, NSF DMR-1506969, and NSF CHE-1800418).

The Gibbs–Duhem equation from classical thermodynamics states

iMNidμi=SdT+Vdp=0,
(A1)

where M represents all chemical species in the system including the solvent, and the last equality holds at constant temperature and pressure.27 Applying the product rule56 to the expression for the Gibbs free energy, G = ∑iMNiμi, we get

dG=iMμidNi+iMNidμi=iMμidNi.
(A2)

Evaluating ∑iMμidNi using stoichiometric coefficients in place of dNi can be shown to lead to the expression ΔG(3).

If the copy numbers are small, we are better served using the discrete difference operator Δ and not the differential difference operator d. The product rule for the discrete difference operator applied to G gives

ΔG=iMμiΔNi+iMNiΔμi+iMΔμiΔNi.
(A3)

Here, we do not neglect the cross-term as we do in the differential product rule [Eq. (A2)]. We evaluate this expression term by term,

iMμiΔNi=ΔG(3),
(A4)
iMNiΔμi=kBTlogiRNiνiNiNijPNi+υjNjNjNN+σN,
(A5)
iMΔNiΔμi=kBTlogiRNiνiNiνijPNi+υjNjυjNN+σσ.
(A6)

Comparing these to the expressions for the accuracy of the various approximations given in  Appendix B, we have

iMNiΔμi+iMΔNiΔμi=ΔG(0)ΔG(3),
(A7)

and thus, combining all the terms in Eq. (A3), we recover our exact expression ΔG(0), and this approach can be seen as an alternate derivation of Eq. (14). To summarize, for small, discrete copy numbers, we obtain corrections to the expression for the change in Gibbs free energy accompanying chemical reactions that are not captured by the constraints imposed by the Gibbs–Duhem equation, which assumes that the copy numbers are continuous quantities.

We calculate the accuracy of the approximations ΔG(1), ΔG(2), and ΔG(3) by taking the difference of these quantities with ΔG(0). We have

(ΔG(1)ΔG(0))/kBT=logN+σNN+σσ,
(B1)
(ΔG(2)ΔG(0))/kBT=logN+σNN+σ+logiRNiNiνiNi×jPNjNj+υjNj,
(B2)
(ΔG(3)ΔG(0))/kBT=logN+σNN+σ+logiRNiNiνiNiνi×jPNjNj+υjNj+υj.
(B3)

The accuracy of ΔG(1) depends only on N for a given reaction, whereas the remaining approximations depend also on the values of Ni and Nj. These observations reflect the fact that, in order to arrive at the expression for ΔG(1), we leveraged the size of N compared to σ, and to arrive at the expressions for ΔG(2) and ΔG(3), we further successively leveraged the sizes of Ni, Nj compared to νi, νj.

Here, we first calculate an approximate expression for the change in Gibbs free energy of the system accompanying a fluctuation of n solvent molecules from compartment A to compartment B. We have

Ginitial=iLNi,Aμi0̃+kBTlogNi,ANA+Ns,Aμs0,*̃+kBTlogNs,ANA+iLNi,Bμi0̃+kBTlogNi,BNB+Ns,Bμs0,*̃+kBTlogNs,BNB,
(C1)

where L is the set of solute species. After the transfer of n solvent molecules from A to B, we have

Gfinal=iLNi,Aμi0̃+kBTlogNi,ANAn+Ns,Aμs0,*̃+kBTlogNs,AnNAn+iLNi,Bμi0̃+kBTlogNi,BNB+n+Ns,Bμs0,*̃+kBTlogNs,B+nNB+n.
(C2)

Taking the difference and simplifying, we arrive at the exact expression

ΔG=kBTlogNANA(NAn)NAn(Ns,An)Ns,AnNs,ANs,ANBNB(NA+n)NA+n×(Ns,B+n)Ns,B+nNs,BNs,B.
(C3)

To understand the magnitude of ΔG for typical values of n, Ns,A, and Ns,B compared to NA and NB, we first make the assumption that the two compartments initially have the same number of solvent molecules, i.e., Ns,A = Ns,BNs. Next, we assume that the solvent molecules dominate the proportion of total molecules, allowing us to write NANBN. These approximations will hold in the limit that the number of solute molecules is much less than the number of solvent molecules for each compartment. We next introduce the small parameters

εA=iLNi,AN
(C4)

and

εB=iLNi,BN,
(C5)

which capture the dilutions of the two compartments, and

ξ=nN,
(C6)

which represents the relative size of the fluctuation. For N = 109, we typically have ξ ∼ 10−4.5 (since nN1/2) and ε ∼ 10−6. Substituting these parameters into Eq. (C3), we have

ΔG=kBTlogN2N(N(1ξ))N(1ξ)(N(1+ξ))N(1+ξ)×(N(1εAξ))N(1εAξ)(N(1εB+ξ))N(1εB+ξ)(N(1εA))N(1εA)(N(1εB))N(1εB).
(C7)

We expand this expression to first order in εA and εB and then to second order in ξ. The result is

ΔG(εAεB)NξkBT+12(εA+εB)Nξ2kBT=(εAεB)nkBT+12(εA+εB)n2NkBT.
(C8)

The observed numerical agreement between Eqs. (C3) and (C8) is close for realistic values of the parameters: for NA = 109, NB = 5 × 108, εA = 10−6, εB = 3 × 10−6, and n=NA1/2, the prediction of Eq. (C3) is −0.0 632 401 kBT and the prediction of Eq. (C8) is −0.0 632 436 kBT. The first term in Eq. (C8) will dominate if the solute dilutions in the two compartments are very different from each other. In this case, we may compare the size of this change in Gibbs free energy to that accompanying the diffusion of a solute from compartment B to compartment A. This latter change in Gibbs free energy will be approximately kBTlogεAεB. If we now set εA = B, where a is of order 1 (typically, it will fall in the range [1/10, 10]), then ΔG for the solvent fluctuation will be kBT(a − 1)εBn and ΔG for the solute diffusion will be kBT log a. The product εBn will be typically on the order of ∼10−1.5, so one can see that for typical values of the parameters the change in Gibbs free energy from a solute diffusion event will be significantly greater in magnitude than that from a solvent fluctuation. If the dilutions are very similar, then εAεB ≈ 0, and the second term in Eq. (C8) dominates. This term is on the order of (εA + εB)kBT since n2N1. These changes in Gibbs free energy will typically be much smaller than those accompanying a chemical reaction or inter-compartment diffusion of the solute. Thus, we may neglect the activity of the solvent in tracking the Gibbs free energy of the system.

1.
R.
Grima
and
S.
Schnell
, “
Modelling reaction kinetics inside cells
,”
Essays Biochem.
45
,
41
56
(
2008
).
2.
N.
Tanaka
and
G. A.
Papoian
, “
Reverse-engineering of biochemical reaction networks from spatio-temporal correlations of fluorescence fluctuations
,”
J. Theor. Biol.
264
(
2
),
490
500
(
2010
).
3.
F.
Baras
and
M.
Malek Mansour
, “
Reaction-diffusion master equation: A comparison with microscopic simulations
,”
Phys. Rev. E
54
(
6
),
6139
(
1996
).
4.
D.
Bernstein
, “
Simulating mesoscopic reaction-diffusion systems using the Gillespie algorithm
,”
Phys. Rev. E
71
(
4
),
041103
(
2005
).
5.
D. T.
Gillespie
, “
Exact stochastic simulation of coupled chemical reactions
,”
J. Phys. Chem.
81
(
25
),
2340
2361
(
1977
).
6.
J.
Schaff
,
C. C.
Fink
,
B.
Slepchenko
,
J. H.
Carson
, and
L. M.
Loew
, “
A general computational framework for modeling cellular structure and function
,”
Biophys. J.
73
(
3
),
1135
1146
(
1997
).
7.
L. M.
Loew
and
C. S.
James
, “
The virtual cell: A software environment for computational cell biology
,”
Trends Biotechnol.
19
(
10
),
401
406
(
2001
).
8.
E.
Roberts
,
J. E.
Stone
, and
Z.
Luthey-Schulten
, “
Lattice microbes: High-performance stochastic simulation method for the reaction-diffusion master equation
,”
J. Comput. Chem.
34
(
3
),
245
255
(
2013
).
9.
M. J.
Hallock
and
Z.
Luthey-Schulten
, “
Improving reaction kernel performance in lattice microbes: Particle-wise propensities and run-time generated code
,” in
2016 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW)
(
IEEE
,
2016
), pp.
428
434
.
10.
J.
Hattne
,
D.
Fange
, and
J.
Elf
, “
Stochastic reaction-diffusion simulation with mesoRD
,”
Bioinformatics
21
(
12
),
2923
2924
(
2005
).
11.
K.
Popov
,
J.
Komianos
, and
G. A.
Papoian
, “
MEDYAN: Mechanochemical simulations of contraction and polarity alignment in actomyosin networks
,”
PLoS Comput. Biol.
12
(
4
),
e1004877
(
2016
).
12.
C.
Floyd
,
G. A.
Papoian
, and
C.
Jarzynski
, “
Quantifying dissipation in actomyosin networks
,”
Interface Focus
9
(
3
),
20180078
(
2019
).
13.
M.
Ander
,
B.
Pedro
,
B.
Di Ventura
,
J.
Ferkinghoff-Borg
,
M. A. F. M.
Foglierini
,
C.
Lemerle
,
I.
Tomas-Oliveira
, and
L.
Serrano
, “
Smartcell, a framework to simulate cellular processes that combines stochastic approximation with diffusion and localisation: Analysis of simple networks
,”
Syst. Biol.
1
(
1
),
129
138
(
2004
).
14.
B.
Drawert
,
S.
Engblom
, and
A.
Hellander
, “
URDME: A modular framework for stochastic simulation of reaction-transport processes in complex geometries
,”
BMC Syst. Biol.
6
(
1
),
76
(
2012
).
15.
S. S.
Andrews
,
J. A.
Nathan
,
R.
Brent
, and
P. A.
Adam
, “
Detailed simulations of cell biology with Smoldyn 2.1
,”
PLoS Comput. Biol.
6
(
3
),
e1000705
(
2010
).
16.
M.
Vigelius
,
A.
Lane
, and
B.
Meyer
, “
Accelerating reaction–diffusion simulations with general-purpose graphics processing units
,”
Bioinformatics
27
(
2
),
288
290
(
2010
).
17.
S.
Wils
and
E.
De Schutter
, “
Steps: Modeling and simulating complex reaction-diffusion systems with python
,”
Front. Neuroinf.
3
,
15
(
2009
).
18.
T. L.
Hill
,
Free Energy Transduction and Biochemical Cycle Kinetics
(
Courier Corporation
,
2004
).
19.
T. M.
Earnest
,
J. A.
Cole
, and
Z.
Luthey-Schulten
, “
Simulating biological processes: Stochastic physics from whole cells to colonies
,”
Rep. Prog. Phys.
81
(
5
),
052601
(
2018
).
20.
D. A.
Beard
and
H.
Qian
, “
Relationship between thermodynamic driving force and one-way fluxes in reversible processes
,”
PLoS One
2
(
1
),
e144
(
2007
).
21.
I.
Prigogine
,
Introduction to Thermodynamics of Irreversible Processes
, 3rd ed. (
Interscience
,
New York
,
1967
).
22.
R.
Rao
and
M.
Esposito
, “
Nonequilibrium thermodynamics of chemical reaction networks: Wisdom from stochastic thermodynamics
,”
Phys. Rev. X
6
(
4
),
041064
(
2016
).
23.
J. L.
England
, “
Dissipative adaptation in driven self-assembly
,”
Nat. Nanotechnol.
10
(
11
),
919
(
2015
).
24.
D. S.
Seara
,
V.
Yadav
,
I.
Linsmeier
,
A. P.
Tabatabai
,
W. O.
Patrick
,
S. M. A.
Tabei
,
S.
Banerjee
, and
M. P.
Murrell
, “
Entropy production rate is maximized in non-contractile actomyosin
,”
Nat. Commun.
9
(
1
),
4948
(
2018
).
25.
J. L.
England
, “
Statistical physics of self-replication
,”
J. Chem. Phys.
139
(
12
),
121923
(
2013
).
26.
I.
Prigogine
and
G.
Nicolis
, “
Biological order, structure and instabilities
,”
Q. Rev. Biophys.
4
(
2-3
),
107
148
(
1971
).
27.
H. B.
Callen
,
Thermodynamics and an Introduction to Thermostatistics
(
Wiley
,
1998
).
28.
P.
Attard
,
Thermodynamics and Statistical Mechanics
(
CUP Archive
,
2002
).
29.
G.
Nicolis
and
I.
Prigogine
,
Self-Organization in Nonequilibrium Systems: From Dissipative Structures to Order through Fluctuations
(
Wiley
,
1977
), pp.
339
426
.
30.
H.
Lodish
,
J. E.
Darnell
,
A.
Berk
,
C. A.
Kaiser
,
M.
Krieger
,
M. P.
Scott
,
A.
Bretscher
,
H.
Ploegh
,
P.
Matsudaira
 et al,
Molecular Cell Biology
(
Macmillan
,
2008
).
31.
R.
Phillips
,
J.
Theriot
,
J.
Kondev
, and
H.
Garcia
,
Physical Biology of the Cell
(
Garland Science
,
2012
).
32.
B.
Alberts
,
D.
Bray
,
J.
Lewis
,
M.
Raff
,
K.
Roberts
, and
J. D.
Watson
,
Molecular Biology of the Cell
, 4th ed. (
Garland
,
2002
).
33.
P. W.
Atkins
,
J.
De Paula
, and
J.
Keeler
,
Atkins’ Physical Chemistry
(
Oxford University Press
,
2018
).
34.
D. A.
Beard
and
H.
Qian
,
Chemical Biophysics: Quantitative Analysis of Cellular Systems
(
Cambridge University Press
,
2008
).
35.
D. A.
McQuarrie
and
J. D.
Simon
,
Physical Chemistry: A Molecular Approach
(
University Science Books
,
Sausalito, CA
,
1997
), Vol. 1.
36.
T. L.
Hill
, “
Perspective: Nanothermodynamics
,”
Nano Lett.
1
,
111
(
2001
).
37.
T. L.
Hill
,
Thermodynamics of Small Systems
(
Courier Corporation
,
1994
).
38.
J.
de Heer
, “
The free energy charge accompanying a chemical reaction and the Gibbs-Duhem equation
,”
J. Chem. Educ.
63
(
11
),
950
(
1986
).
39.
T.
Engel
,
P.
Reid
 et al,
Physical Chemistry
(
Pearson
,
2006
).
40.

The argument of the logarithm, NAΘA, has dimensions of concentration, however it is understood that such quantities are taken with reference to some standard concentration, typically 1M.

41.
D. M.
Anderson
,
J. D.
Benson
, and
A. J.
Kearsley
, “
Foundations of modeling in cryobiology—I: Concentration, Gibbs energy, and chemical potential relationships
,”
Cryobiology
69
(
3
),
349
360
(
2014
).
42.
G. M.
Barrow
, “
Free energy and equilibrium: The basis of G0 = -RT in K for reactions in solution
,”
J. Chem. Educ.
60
(
8
),
648
(
1983
).
43.

We have logN+σN(N+σ)=log1+σN(N+σ)σ for large N.

44.
Y.
Mishin
, “
Thermodynamic theory of equilibrium fluctuations
,”
Ann. Phys.
363
,
48
97
(
2015
).
45.
R.
Milo
and
R.
Phillips
,
Cell Biology by the Numbers
(
Garland Science
,
2015
).
46.
S.
Engblom
,
L.
Ferm
,
A.
Hellander
, and
P.
Lötstedt
, “
Simulation of stochastic reaction-diffusion processes on unstructured meshes
,”
SIAM J. Sci. Comput.
31
(
3
),
1774
1797
(
2009
).
47.
S. A.
Isaacson
and
C. S.
Peskin
, “
Incorporating diffusion in complex geometries into stochastic chemical kinetics simulations
,”
SIAM J. Sci. Comput.
28
(
1
),
47
74
(
2006
).
48.
B.
Bayati
,
P.
Chatelain
, and
P.
Koumoutsakos
, “
Adaptive mesh refinement for stochastic reaction–diffusion processes
,”
J. Comput. Phys.
230
(
1
),
13
26
(
2011
).
49.
L.
Ferm
,
A.
Hellander
, and
P.
Lötstedt
, “
An adaptive algorithm for simulation of stochastic reaction–diffusion processes
,”
J. Comput. Phys.
229
(
2
),
343
360
(
2010
).
50.
M. B.
Flegg
,
S. J.
Chapman
, and
R.
Erban
, “
The two-regime method for optimizing stochastic reaction–diffusion simulations
,”
J. R. Soc., Interface
9
(
70
),
859
868
(
2011
).
51.
M. A.
Gibson
and
J.
Bruck
, “
Efficient exact stochastic simulation of chemical systems with many species and many channels
,”
J. Phys. Chem. A
104
(
9
),
1876
1889
(
2000
).
52.
J.
Borge
, “
Reviewing some crucial concepts of Gibbs energy in chemical equilibrium using a computer-assisted, guided-problem-solving approach
,”
J. Chem. Educ.
92
(
2
),
296
304
(
2014
).
53.
R.
Erban
,
M. B.
Flegg
, and
G. A.
Papoian
, “
Multiscale stochastic reaction–diffusion modeling: Application to actin dynamics in filopodia
,”
Bull. Math. Biol.
76
(
4
),
799
818
(
2014
).
54.
D.
Ulrich
,
G. A.
Papoian
, and
R.
Erban
, “
Steric effects induce geometric remodeling of actin bundles in filopodia
,”
Biophys. J.
110
(
9
),
2066
2075
(
2016
).
55.
Y.
Lan
and
G. A.
Papoian
, “
The stochastic dynamics of filopodial growth
,”
Biophys. J.
94
(
10
),
3839
3852
(
2008
).
56.

The product rule can be derived as d(fg) = (f + df) (g + dg) − fg = fdg + gdf + dfdg = fdg + gdf since the term dfdg is assumed to be negligible. We give this reminder to emphasize the presence of the cross term, which we do not neglect in the discrete case.