
==== Front
ArXiv
ArXiv
arxiv
ArXiv
2331-8422
Cornell University

arXiv:2409.01529v1
2409.01529
1
preprint
Article
Algebraic and diagrammatic methods for the rule-based modeling of multi-particle complexes
Rousseau Rebecca J. 1*
Kinney Justin B. 2†
1 Department of Physics, California Institute of Technology, Pasadena, CA 91125
2 Simons Center for Quantitative Biology, Cold Spring Harbor Laboratory, Cold Spring Harbor, NY 11724
* rroussea@caltech.edu
† jkinney@cshl.edu
3 9 2024
arXiv:2409.01529v1https://creativecommons.org/licenses/by/4.0/ This work is licensed under a Creative Commons Attribution 4.0 International License, which allows reusers to distribute, remix, adapt, and build upon the material in any medium or format, so long as attribution is given to the creator. The license allows for commercial use.
nihpp-2409.01529v1.pdf
The formation, dissolution, and dynamics of multi-particle complexes is of fundamental interest in the study of stochastic chemical systems. In 1976, Masao Doi introduced a Fock space formalism for modeling classical particles. Doi’s formalism, however, does not support the assembly of multiple particles into complexes. Starting in the 2000’s, multiple groups developed rule-based methods for computationally simulating biochemical systems involving large macromolecular complexes. However, these methods are based on graph-rewriting rules and/or process algebras that are mathematically disconnected from the statistical physics methods generally used to analyze equilibrium and nonequilibrium systems. Here we bridge these two approaches by introducing an operator algebra for the rule-based modeling of multi-particle complexes. Our formalism is based on a Fock space that supports not only the creation and annihilation of classical particles, but also the assembly of multiple particles into complexes, as well as the disassembly of complexes into their components. Rules are specified by algebraic operators that act on particles through a manifestation of Wick’s theorem. We further describe diagrammatic methods that facilitate rule specification and analytic calculations. We demonstrate our formalism on systems in and out of thermal equilibrium, and for nonequilibrium systems we present a stochastic simulation algorithm based on our formalism. The results provide a unified approach to the mathematical and computational study of stochastic chemical systems in which multi-particle complexes play an important role.
==== Body
pmcI. INTRODUCTION

Large complexes of classically behaving particles play a central role in a variety of scientific disciplines. For example, many essential biological processes depend on large complexes formed by proteins, nucleic acids, and/or other macromolecules. A common theme in such systems is “combinatorial complexity” [1], i.e., that an immense (and often infinite) variety of molecular complexes can form from a relatively small number of interaction rules governing the assembly of a relatively small number of molecular components. Nevertheless, mathematical methods for analyzing stochastic chemical systems that exhibit such combinatorial complexity have yet to be developed.

In 1976, Masao Doi [2, 3] introduced a Fock space formalism for modeling many-body systems of classical particles. This approach was further developed by others [4–6], and has proven useful in the study of diffusion-limited aggregation [7–9] and other problems in statistical physics [10–12]. As in quantum field theory, Doi’s formalism supports the creation and annihilation of particles, but does not support the assembly of preexisting particles into complexes. Consequently, analyzing systems that involve multi-particle complexes using this formalism requires specifying one distinct field for every distinct species of complex. This makes Doi’s formalism unwieldy for analyzing systems that exhibit substantial combinatorial complexity.

Consider, for example, the homopolymer system illustrated in Fig. 1. This system comprises one type of component particle having two sites capable of forming an interaction. In thermal equilibrium the system’s behavior is governed by two quantities: the chemical potential of the particles and the energy of interaction [Fig. 1(a)]. Out of equilibrium the system is governed by four rate parameters describing the appearance and disappearance of particles, as well as their mutual binding and unbinding [Fig. 1(b)]. But despite how simple this system is to describe in words and pictures, modeling this system in Doi’s formalism is complicated because the above rules lead to an infinite number of possible polymeric complexes. To apply Doi’s formalism, one must define an infinite number of fields, one for every species of complex. For equilibrium systems, one must then manually specify the chemical potential of each species [Fig. 1(c)]. For systems out of equilibrium, one must manually specify the rate of reaction between all reacting sets of species [Fig. 1(d)]. And in doing so, one must take care that the chemical potentials and/or reaction rates written down are expressed correctly as functions of the underlying model parameters, as Doi’s formalism provides no means of computing these quantities.

Clearly something is missing. Ideally, the formalism one uses to describe systems of multi-particle complexes should allow one to mathematically derive the set of possible complexes, the chemical potentials of each complex, and the rates of reaction between complexes from the underlying rules and their associated parameters. This paper develops a formalism that does this.

The problem of combinatorial complexity has long been recognized in the field of computational systems biology. Starting in the 2000s, researchers studying biological signaling pathways began developing “rule-based” approaches for simulating chemical systems of multi-particle complexes [1, 13–19]. Some of these efforts have produced sophisticated software ecosystems, such as BioNetGen [14, 16, 18, 20] and Kappa [17, 21]. These simulation approaches, however, are based on formal representations that do not lend themselves to analytical calculations using the mathematical methods of statistical physics. For example, BioNetGen is based on a process algebra describing the formation and dynamics of port graphs (i.e., graphs with edges attached through ports) [22] while Kappa is based on the κ-calculus process algebra [17, 21]. As a result, work in this area has been confined to computational analyses rather than analytic calculations.

Here we bridge the divide between Doi’s mathematical formalism and rule-based methods for computationally modeling biochemical systems. Using a Fock space for classical particles reminiscent of but distinct from that of Doi, we develop an operator algebra that allows not only for the creation and annihilation of particles, but also for the assembly of particles into complexes. We show that this operator algebra allows one to mathematically analyze equilibrium and nonequilibrium systems that are defined in a rule-based manner, and can also be used as a basis for computational analyses using stochastic simulations.

After introducing how microstates and macrostates are represented, we apply the formalism to three systems in thermal equilibrium: a monomer system, a homodimer system, and a homopolymer system. Next we show how our formalism can be used to define and analyze nonequilibrium systems, both through the analytic derivation of master equations and through computational analysis carried out using a stochastic simulation algorithm. We end by illustrating the versatility of our formalism, showcasing the variety and complexity of system behavior that can arise from positing different sets of rules.

II. FOUNDATIONS

A. Microstates

Following Doi we define a set 𝒮 of microstates where each microstate s∈𝒮 corresponds to a unit vector |s⟩. The resulting set of pure states forms an orthonormal basis for the Fock space. The vector |ψ⟩ describing the system is then given by a probabilistic mixture of pure states, (1) ψ=∑s∈𝒮pss,

where ps represents the probability of the system being in state s. It is useful to define the sum of all possible states as the “sum vector” (2) ∣sum⟩=∑s∈𝒮|s⟩.

The expectation value of any operator O is then given by ⟨sum|O|ψ⟩, and probability normalization requires that ⟨sum∣ψ⟩=1.

The Fock space supports systems both in and out of thermal equilibrium. In equilibrium, the state vector |ψ⟩ for the system (which is taken to be in the grand canonical ensemble) can be expressed in terms of a Hamiltonian operator H that assigns a free energy to each microstate: (3) |ψ⟩=e-βHZ|sum⟩,whereZ=sume-βHsum

is the partition function, β=1/kBT where kB is Boltzmann’s constant and T is temperature, and H|s⟩=Hs|s⟩ where Hs denotes the free energy of state s. The dynamics of the system state |ψ⟩ out of equilibrium is described by (4) ddtψ=Wψ,

where W is a transition operator. We call this the “macrostate master equation.” In terms of the scalar transition rates Ws→t from microstate s to microstate t, the transition operator is (5) W=∑s,t∈𝒮Ws→t(|t⟩⟨s|-|s⟩⟨s|).

Note that the first term in the summand reflects the flow of probability into pt, while the second term reflects the flow of probability out of ps. We call these the “reaction” and “depletion” terms, respectively.

B. Macrostates

Unlike in Doi’s formalism, the microstates in our formalism represent not only the externally observable properties of particles, but also their unobservable internal states. These internal states are, in fact, what make particles in our formalism identifiable and thus allow complexes to be constructed from preexisting particles. We therefore distinguish between the microstates of a system (represented by the |s⟩ vectors) and the macrostates of the system.

In this work we focus on zero-dimensional (i.e., well-mixed) populations of particles and complexes. Assuming there are K possible observably distinct species of complex, each macrostate is characterized by a vector n→=n1,…,nK where each nk quantifies the number of complexes of species k. The corresponding macrostate vector is defined to be the sum of all microstate vectors consistent with the macrostate, i.e., (6) n→=∑s∣n→s.

Note that |sum⟩=∑n→|n→⟩, and that the probability of a the system being in a macrostate n→ given |ψ⟩ is (7) Pn→=n→ψ.

In equilibrium systems, each macrostate is an eigenstate of the Hamiltonian: (8) Hn→=-∑knkμkn→,

where μk is the (bare) chemical potential of species k. Consequently, the probability of the system having macrostate n→ is a product of species-specific Poisson distributions, i.e., (9) Pnk=1Zkeβμknknk!,whereZk=expeβμk.

In nonequilibrium systems, this probability becomes a function of time t and evolves according to the “macrostate master equation” (10) ddtPtn→=∑n→′Wn→n→′Ptn→′,

where Wn→n→′=n→|W|n→′/n→′|n→′ are the macrostate-specific transition rates. In later sections we compute these Wn→n→′ from the transition operator W, but instead of calculating each rate directly we find it simpler to calculate the vector (11) Jn→=W†n→.

We call |J(n→)⟩ the “flux projector” since taking the inner product of it with |ψ⟩ yields a vector of probability fluxes, i.e., P˙t(n→)=⟨J(n→)|ψ(t)⟩.

III. MONOMER IN EQUILIBRIUM

A. Microstates and macrostates

We now construct the Fock space on which our formalism is based, using a system of monomeric particles for concreteness. The particles are represented using a hard-core boson field, A, which is assumed to have N excitation modes. Each mode Ai is indexed by a number i∈𝒩={1,…,N} that represents the internal state of a particle. This index allows the formalism to track individual particles that are outwardly identical. In what follows we keep N finite for concreteness, but all physically meaningful calculations are performed in the N→∞ limit.

Each mode Ai can be in one of two orthonormal states: |1⟩i represents the presence of a particle with internal state i, while |0⟩i represents its absence. Microstates are given by tensor products over all modes. Specifically, a microstate representing K particles having indices ℐ⊆𝒩 is represented by (12) |ℐ⟩=⨂i∈𝒩|1⟩iifi∈ℐ,|0⟩iotherwise.

These states are orthonormal, i.e., ⟨ℐ∣𝒥⟩=δℐ𝒥. The resulting macrostates of the system are (13) n=∑ℐ:ℐ=nℐ,forn=0,…,N.

The vacuum state, |0⟩=|∅⟩, is both a microstate and a macrostate.

These definitions readily extend to systems defined by multiple fields. Consider a mixture of monomer species A and B. The microstate comprising A monomers having indices ℐ and B monomers having indices 𝒥 is given by (14) |ℐ,𝒥⟩A,B=|ℐ⟩A⊗|𝒥⟩B,

where the subscripts indicate the Fock space in which each state vector lives. The corresponding macrostates and sum states are given by analogous tensor products. Systems with three or more fields are defined similarly.

B. Mode and field operators

We now define four types of mode-specific operators: creation, annihilation, presence, and absence. The creation operator for mode i is defined to be Aˆi=|1⟩i0i. When applied to a microstate |ℐ⟩, this operator has the effect (15) Aˆi|ℐ⟩=0ifi∈ℐ,|ℐ∪{i}⟩otherwise.

The corresponding annihilation operator is defined to be Aˇi=Aˆi†, and has the effect (16) Aˇi|ℐ⟩=|ℐ\{i}⟩ifi∈ℐ,0otherwise.

The presence operator is defined as A‾i=AˆiAˇi. A‾i|ℐ⟩ is one if i∈ℐ and zero otherwise. The absence operator is defined to be A˜i=AˆiAˇi=1-A‾i. Note that A‾i and A˜i are self-adjoint.

We highlight several key algebraic properties of these operators. First, creation and annihilation operators are nilpotent, i.e., Aˆi2=Aˇi2=0. Second, the commutator (17) [Aˇi,Aˆj]=δij1-2A‾i

is very different than one finds in the harmonic oscillator algebra, and thus in other algebras used to model classical particles [2, 6]. Third, the commutator (18) [A‾i,Aˆj]=δijAˆj,

plays an important role later when constructing multi-particle complexes from component particles. Appendix A lists some additional useful properties.

We further define field-specific creation, annihilation, presence, and absence operators as sums over all of the corresponding mode operators, i.e., (19) Aˆ=∑iAˆi,

and similarly for Aˇ, A‾, and A˜. These field operators satisfy the useful commutation relations (20) [Aˇ,Aˆ]=N-2A‾,[A‾,Aˆ]=Aˆ.

Macrostates are given by (21) n=Aˆnn!0.

Here the combinatorial factor corrects for each set of modes being summed over n! times in the operator product Aˆn. Note that |n⟩=0 if n>N, since this would cause every term in Aˆn to contain at least one factor of Aˆi2. Applied to a macrostate, one finds that (22) Aˆn=n+1n+1,A‾|n⟩=n|n⟩,

(23) Aˇn=N-n+1n-1,A˜|n⟩=(N-n)|n⟩.

See Appendix A for a derivation of these results.

In the large N limit, the number of modes that are excited in any macrostate |n⟩ with substantial physical probability becomes negligible compared to N. In what follows we therefore approximate Eq. (23) as (24) Aˇn≈Nn-1,A˜n≈Nn.

By similar logic we can also approximate [Aˇ,Aˆ]≈N, etc.

Finally, it is useful to consider the coherent state, (25) z=∑n=0∞znn=ezA0.

This allows one to express the generating function for the distribution over macrostates as ⟨z∣ψ⟩. Note that setting z=1 recovers the sum state, i.e., |1⟩=|sum⟩.

C. Hamiltonian operator

In what follows we assume that each system of interest is contained within a volume V. For a gas of monomers, the relevant Hamiltonian is H=-μA‾, where μ denotes a bare chemical potential. We use the term “bare” to emphasize that μ determines the excitation probability for each independent mode Ai, whereas the concentration of monomers depends on μ, N, and V . The generating function for the equilibrium probability distribution over macrostates is found by: (26) ⟨z∣ψ⟩=1Z⟨0|ezAˇeβμA‾eAˆ|0⟩

(27) =1Z⟨0|ezAˇeλAˆ|0⟩(definingλ=eβμ)

(28) =1Z∏i⟨0iezAˇieλAˆi∣0⟩i

(29) =1Z∏i⟨0i(1+zAˇi)(1+λAˆi)∣0⟩i

(30) =1Z(1+zλ)N

(31) =1+zλ1+λN.

In the first step we used the fact that f(A‾)g(Aˆ)|0⟩=g(f(1)Aˆ)|0⟩ for any functions f and g. The resulting quantity λ is the per-mode fugacity corresponding to chemical potential μ. In the second step we used the fact that operators for different modes commute. In the third step we used the nilpotency of Aˆi to truncate the expansions of each exponential. Finally, we used the normalization requirement ⟨1∣ψ⟩=1 to determine the partition function Z=(1+λ)N. The result is the generating function for the binomial distribution corresponding to N modes with a per-mode excitation probability of λ/(1+λ). From this generating function we find that the expected concentration of monomers is (32) ⟨A‾⟩V=1Vddz⟨z∣ψ⟩z=1=1VNλ1+λ=λ′1+Vλ′/N,

where λ′=NVλ. Keeping ⟨A‾⟩/V constant while taking N→∞ requires holding λ′ approximately constant and thus rescaling λ∼VN. In this limit we get (33) zψ=ez-1Vλ′⇒ψ=e-Vλ′Vλ′.

The corresponding partition function is Z=eVλ′. Note that ⟨z|ψ⟩ is the generating function for a Poisson distribution with mean Vλ′. We thus see that (34) μ′=kBTlogλ′=μ+kBTlogNV

is the effective chemical potential, i.e., the chemical potential appropriately renormalized to account for the N modes available for excitation in volume V . It is therefore μ′, not μ, that reflects the physically measurable chemical potential.

IV. HOMODIMER IN EQUILIBRIUM

A. Composite operators

Multi-particle complexes are represented as products of mode operators for three kinds of fields: particle fields, interaction fields, and site fields. For example, we define the creation operator for a dimer of two A particles by the composite operator (35) Dˆij=IˆijaˆiaˆjAˆiAˆj,

which is the product of mode operators for a particle field A, an interaction field I, and a site field a. More specifically, Aˆi and Aˆj create the two component particles, Iˆij registers that these two particles interact with one another, and aˆi and aˆj respectively indicate that the Ai and Aj particles are each participating in an interaction and are therefore not free to interact with additional particles. Note that the index of the dimer creation operator is the pair of monomer indices, (i,j). Since the monomer is symmetric, we assume that Iˆij=Iˆji, and so Dˆij=Dˆji. Note also that Dˆii=0 because of the nilpotency of Aˆi and aˆi. The number of internal states for the dimer is therefore ND=N2≈N2/2. The dimer annihilation, presence, and absence operators are defined in terms of the creation operator in the same manner as for a single particle: (36) Dˇij=Dˆij†=IˇijaˇiaˇjAˇiAˇj,

(37) D‾ij=DˆijDˇij=I‾ija‾ia‾jA‾iA‾j,

(38) D˜ij=DˇijDˆij=I˜ija˜ia˜jA˜iA˜j.

The field operator Dˆ is given by Dˆ=12∑i,jDˆij, where the factor 1/2 compensates for double counting in the sum. The field operators Dˇ, D‾, and D˜ are defined similarly.

The homodimer system also comprises free monomers. We represent these by a separate composite field M defined by the mode operator Mˆi=Aˆia˜i. The corresponding number of internal states is NM=N, and the three related mode operators are Mˇi=Aˇia˜i, M‾i=A‾ia˜i, and M˜i=A˜ia˜i. The corresponding field operators are defined as sums over i. Note the inclusion of a˜i in Mˆi ensures that M‾ does not count A particles that are components of dimers, Mˇ does not annihilate such particles (which would leave dangling I and a modes), etc.

The macrostate comprising m monomers and d dimers is given by (39) m,d=Mˆmm!Dˆdd!0.

Using a˜i,aˆj=-δijaˆi and Aˆi2=0, one can readily verify that Mˆ and Dˆ commute. The corresponding coherent state is therefore (40) zM,zD=∑m=0∞∑d=0∞zMmzDdm,d=ezMMˆ+zDDˆ0.

B. Sectoring by species

Now consider a Hamiltonian in which each A particle has chemical potential μ and each interaction has Gibbs free energy ϵ: (41) H=-μ∑iA‾i+ϵ12∑i,jI‾ij.

To compute the equilibrium state of the system, we re-express the Hamiltonian as a sum of terms that operate separately on monomers and dimers. Using the identity 1=a˜i+a‾i, we split the Hamiltonian into two parts, H=HM+HD, where (42) HM=-μ∑iA‾ia˜i,HD=-μ∑iA‾ia‾i+ϵ2∑i,jI‾ij.

These operators satisfy the commutation relations (43) [HM,Mˆ]=-μMMˆ,[HD,Dˆ]=-μDD,ˆ[HM,Dˆ]=0,[HD,Mˆ]=0,

where μM=μ and μD=2μ-ϵ are the bare monomer and dimer chemical potentials. Next we compute the generating function: (44) zM,zD∣ψ=Z-1⟨0|ezMMˇ+zDDˇe-βHM+HDeMˆ+Dˆ|0⟩=Z-1⟨0|ezMMˇ+zDDˇeλMMˆ+λDDˆ|0⟩≈Z-1⟨0|ezMMˇeλMMˆezDDˇeλDDˆ|0⟩.

In the first step we used Eq. (43) and defined the fugacities λM=eβμM and λD=eβμD. The second step follows from the approximation (see Appendix B) (45) [Mˆ,Dˇ]=12∑i,j(Mˆi+Mˆj)Dˇij≈0.

This commutator is not exactly zero because annihilating a dimer frees up A modes that can be used to create two monomers. But by way of comparison, (46) [Mˆ,Mˇ]=NM-2M‾and[Dˆ,Dˇ]=ND-2D‾

have terms that scale as N and N2, respectively. The effect of the commutator [Mˆ,Dˇ] on a physical state is consequently negligible in the large N limit. This reflects the number of modes available to create a monomer not being limiting in the physically meaningful regime.

Next we insert a copy of the identity operator, (47) 1=∑m,d|m,d⟩⟨m,d|m!d!,

into the right-hand side of Eq. (44) and observe that only the |0⟩⟨0| term survives. Consequently, (48) zM,zD∣ψ≈Z-1⟨0|ezMMˇeλMMˆ|0⟩⟨0|ezDDˇeλDDˆ|0⟩=Z-1zM∣λMzD∣λD.

Setting ⟨1,1|ψ⟩=1 we find that Z≈ZMZD, where ZM and ZD are the respective partition functions for the monomer and dimer species. We thus obtain (49) |ψ⟩≈|ψ⟩M⊗|ψ⟩D,

where |ψ⟩M describes a monomer-only system, |ψ⟩D describes a dimer-only system, and both have the same Poisson form as in Eq. (33).

There are two important caveats to the result in Eq. (49). First, the sectors for distinct species are only independent in the N→∞ limit. For example, the right-hand side of Eq. (49) has nonzero |m⟩⊗|d⟩ terms for all values of m≤N and d≤N/2, whereas each |m,d⟩ term on the left-hand side is nonzero only if m+2d≤N. Second, this sectoring result holds only in equilibrium systems; indeed, the populations of particles in different sectors will generally be coupled out of equilibrium.

Finally we discuss the scaling behavior of the system with N and V. As in the previous section, requiring the concentration of monomers ⟨M‾⟩/V to be constant as N→∞ reveals an effective monomer chemical potential of μM′=μM-+kBTlogNMV. Similarly, requiring the concentration of homodimers ⟨D‾⟩/V to be constant as N→∞ reveals an effective monomer chemical potential of μD′=μD+kBTlogNDV. These relations are realized by renormalizing the parameters of the Hamiltonian so that (50) μ′=μ+kBTlogNVandϵ′=ϵ-kBTlogV

are held constant. In terms of these quantities, the effective dimer chemical potential is μD′=2μ′-ϵ′-kBTlog2, where the logarithmic term accounts for the symmetry of the molecule. This system is therefore exactly renormalizable. We thus see that, to maintain a fixed concentration of dimers, the bare interaction energy ϵ must become weaker as system volume increases. This makes sense: if V increases while monomer concentration stays fixed, the number of monomers available to bond to a given monomer will increase in proportion to V. To keep the probability of the given monomer forming a dimer constant, the bare interaction energy ϵ must weaken as V increases so that e-βϵ∝V-1. This implies that e-βϵ′=Ve-βϵ will be fixed.

C. Gallery operators

Consider more generally a system that realizes K distinct species of complex. Let Gˆk denote the creation operator for complex k and assume that (51) [G‾k,Gˆk′]=δkk′Gˆk.

We refer to the vector G→=(Gˆ1,…,Gˆk)⊤ as the “gallery,” as it exhibits creation operators for all possible complexes. The gallery allows us to define the coherent state (52) z→=expz→⊤G→0,

where z→=z1,…,zK⊤ is a vector of scalars. As in the monomer and homodimer systems, the generating function for a system |ψ⟩ is ⟨z→|ψ⟩, and the sum state is (53) sum=e∑kGˆk0=|1→⟩.

The macrostates of the system are given by (54) n1,…,nK=∏k=1KGˆknknk!0,

and the effects of the four field operators on macrostates are (55) Gˆkn1,…,nK=nk+1n1,…,nk+1,…,nK,G‾kn1,…,nK=nkn1,…,nK,Gˇkn1,…,nK≈Nkn1,…,nk-1,…,nK,G˜kn1,…,nK≈Nkn1,…,nK.

where Nk is the number of internal states for species k.

Because different complexes can share internal components, the Hamiltonian can be defined in a rule-based manner as in Eq. (41) instead of on a species-by-species basis. As we will see in later sections, such rule-based definitions can require far fewer than K terms. If there are no energetic interactions between separate complexes, the Hamiltonian can then be equivalently expressed as (56) H≃-∑kμkG‾k.

where μk is the bare chemical potential for species k. The generating function for the equilibrium state then factorizes, i.e. (57) ψ=⨂k|ψ⟩k,where|ψ⟩k=e-Vλk′Vλk′,

and where λk′=eβμk′ and μk′=μk+kBTlogNkV are the effective fugacity and chemical potential of species k. We note that it may or may not be possible to renormalize the parameters of the Hamiltonian so that all the effective chemical potentials are independent of V. For example, this is possible in the homodimer system, but not in the homopolymer system discussed in the next section.

D. Factory operators

Hamiltonians of the form in Eq. (56) describe systems of non-interacting particles and might understandably be viewed as trivial. They become less trivial, however, in systems comprising large (or infinite) numbers of distinct complexes, each complex having a chemical potential that is a function of the parameters used to define the Hamiltonian. In such systems, merely enumerating different species of complex and determining their chemical potentials can be nontrivial. It is therefore natural to instead define the set of possible complexes implicitly by specifying the rules for their construction, and to use these rules to then compute the different species of complex and their associated chemical potentials. We now show how our formalism enables this.

In the case of the homodimer, the sum state can be expressed as (58) sum=eF2eF10,

where (59) F1=∑iAˆia˜i,F2=12∑i,jIˆijaˆiaˆjA‾iA‾j.

Here, F1=Mˆ creates free monomers, while F2 joins two monomers into a dimer. Specifically, F2 tests for the presence of two particles, Ai and Aj, and if these already exist it joins them into a dimer Dij. Note that neither A particle can be part of an existing dimer due to the excitation of site fields ai and aj. For example, (60) F2F122|0⟩=14∑i,j,k,lIˆijaˆiaˆjA‾iA‾jAˆka˜iAˆla˜l|0⟩

(61) =14∑i,j,k,lIˆijaˆiaˆjAˆiAˆjδikδjl+δilδjk|0⟩

(62) =0,1,

where the first step uses the identities aˆia˜i=aˆi and [A‾i,Aˆj]=δijAˆi. We will soon show more generally that (63) F2pp!F1qq!|0⟩=|q-2p,p⟩ifq≥2p,0otherwise.

Summing this over all p and q establishes the |sum⟩ state in Eq. (58).

The sum of states for complexes generated in any system can thus be specified by a vector of operators F→=F1,…,FL⊤ via (64) sum=eFL…eF10.

We call this vector the “factory.” We emphasize that the order of the operators within the factory is important, as these operators, unlike gallery operators, do not generally commute.

Using the factory instead of the gallery to define the set of possible complexes in a system can have an important advantage: the factory often comprises far fewer operators than the gallery. This is not the case for the homodimer system, but it is so for the homopolymer system presented in Section V.

There is a disadvantage, however, to defining a system using the factory: one loses access to the generating function. One can, of course define a coherent state analogous to |z→⟩ via (65) x→=exLFL⋯ex1F10,

where x→=x1,…,xL⊤. It is questionable, however, how useful the corresponding generating function ⟨x→|ψ⟩ is for analysis. As we will see, each term in the expansion of Eq. (65) can yield multiple distinct mixtures of complexes. One thus generally cannot read off the macrostate distribution P(n→) from the expansion of ⟨x→|ψ⟩.

E. Wick’s theorem

The algebraic manipulations needed to show Eq. (63) become unwieldy as p and q become large. Wick’s theorem, a foundational result in quantum field theory, makes these calculations significantly more straightforward by providing a systematic procedure for reordering operators in a multi-operator product.

To see how Wick’s theorem can be applied to our formalism, define the compound operators A‾ia=aˆiA‾i and Aˆia=aˆiAˆi. Ignoring the interaction field I for the moment, each term in the expansion of the left-hand side of Eq. (63) has the form (66) A‾i1a⋯A‾i2qaMˆj1⋯Mˆjp0.

Note that all compound presence operators appear to the left of all creation operators. We refer to this as “productive ordering.” Given an operator product X1X2⋯Xn, we denote its productive ordering by 𝒫X1X2⋯Xn. Each term in the Taylor expansion of the factory representation is productive ordered because the instructions for assembling each complex are applied after the instructions for creating its components. In contrast, each term in the expansion of the gallery representation in Eq. (39) has the form (67) Mˆi1⋯AˆimMˆj1a⋯Aˆj2da0.

The key difference from Eq. (66) is that this term contains only creation operators, all of which commute.

Every term of the form in Eq. (66) is in fact equal to a sum of terms having the form in Eq. (67) with m=q-2p and d=p. To transform the former to the latter, we iteratively apply the exchange rule (68) A‾iaMˆj=MˆjA‾ia+δijAˆia

until no A‾ia operators appear to the left of any Mˆj operators. Each application of the exchange rule adds another term to the expansion. The result is a sum of operator products such that all A‾ia in each product appear to the right of all Mˆi and Aˆia. Such products are said to be “normally ordered.” More generally, an operator product is normally ordered if all presence operators appear to the right of all creation operators. The normally ordered form of an operator product X1X2⋯Xn is denoted by 𝒩X1X2⋯Xn. Normally ordered products are useful because any such products containing presence operators vanish when applied to the vacuum state.

Wick’s theorem provides an equality between productive ordered and normally ordered operator products. To state Wick’s theorem, we define a “contraction” between two operators Xi and Xj to be (69)

The contraction of two specific operators within a larger product X1X2⋯Xn removes these operators from the product and replaces them with their contraction, i.e., (70)

A key assumption of Wick’s theorem is that the contraction of any two operators in a product is “central,” i.e., it commutes with all other operators in the product. For the homodimer algebra, the only contractions needed to transform Eq. (66) to Eq. (67) are of the form (71)

These contractions are indeed central, i.e., [Aˆja,Mˆj]=[Aˆja,A‾ia]=0.

Applied to our context, Wick’s theorem states that any productive ordered operator product is equal to the sum of all possible normally ordered contractions: (72)

(73) =∑allcontractions𝒞𝒩𝒞X1X2…Xn

For example, applying Wick’s theorem to the right-hand side of Eq. (60) gives (74)

(75) =MˆkMˆkA‾iaA‾ja+δikAˆiaMˆkA‾ja+δilAˆiaMˆkA‾ja+δjkAˆjaMˆkA‾ia+δjlAˆjaMˆkA‾ia+δikδjlAˆiaAˆja+δilδjkAˆiaAˆja.

When applied to the vacuum state, only the last two terms in Eq. (75) survive, thus yielding the expression in Eq. (61), (76) A‾iaA‾jaMˆkMˆk0=δikδjl+δilδjkAˆkaAˆla0.

More generally, Wick’s Theorem allows us to transform terms in the expansion of the factory expression for the sum vector (Eq. (64)) into a sum of terms in the expansion of the gallery expression for the sum vector (Eq.(53)).

F. Formal diagrams

We now introduce diagrammatic methods that aid in computations involving Fock space operators. Each diagram indicates an operator product or sums of such products over internal states. Fig. 2 shows several examples. The indices of mode operators are written as index names inside open dots [Fig. 2(a)]. Mode operators are indicated by the decorated operator name written next to their respective dots. Multiple operator names written next to the same dot indicate that those operators share the same index. Modes that have two indices are written next to lines that connect the two dots representing these indices. A closed dot indicates summation over the corresponding index. Field operators are thus distinguished from mode operators through the use of closed rather than open dots. Symmetry factors are also kept explicit [Fig. 2(b)].

This diagrammatic notation is helpful in computations involving Wick contractions; we demonstrate this by deriving Eq. (63). The effect of each Wick contraction is illustrated in Fig. 2(c): contracting an A‾a operator (part of F2) with F1=Mˆ eliminates the F1 and replaces the A‾a in F2 with an Aˆa. To avoid unnecessary notation going forward, we represent this operation using the same diagrams but showing only the decorations on the A operators. Now consider the left-hand side of Eq. (63) with p=2 and q=5 [Fig. 3(a), line 1]. Because the two F2 operators are applied after the five F1 operators, this product is productive ordered [Fig. 3(a), line 2]. Next we use Wick’s theorem to convert this to a sum of normally ordered products [Fig. 3(a), lines 3–8]. We also evaluate the combinatorial coefficients that arise due to distinct contractions producing topologically identical products. Consider, for example, the coefficient for term (iii). For the first contraction, there are five choices of Mˆ and four choices of A‾a. For the second contraction, there are four remaining choices of Aˆ but only one possible choice for A‾a – that which is linked with the first A‾a through an I-field. Since interchanging the order in which the contractions are performed does not change the result, we divide the result by two. This yields a combinatorial coefficient of (5·4)×(4·1)/2!. The combinatorial coefficients for the other terms follow similarly.

Each normally ordered term in lines 3–8 yields one of the resulting operator products shown in lines 9–11, as indicated. We leave it to the reader to check that the combinatorial coefficients computed in lines 3–8 do in fact match those shown in lines 9–11, which are as expected based on symmetry considerations. Note in particular that all terms except term (vi) contain A‾a operators. Consequently, only term (vi) survives when applying this result to the vacuum state [Fig. 3(b)].

We are now in a position to evaluate the left-hand side of Eq. (63) for general values of p and q (Fig. 4). If q<2p, then F1q does not supply enough Mˆ operators to contract all the A‾a operators supplied by F2p. The expression therefore vanishes. If q≥2p, however, there are q!/(q-2p)! ways to contract all the A‾a with all the Mˆ, thereby leaving p copies of Dˆ and q copies of Mˆ. The resulting combinatorial factor replaces the 1/q! with 1/(q-2p)!, thus providing the factors needed to correct the redundancies in Dˆp and Mˆq. This completes the derivation of Eq. (63), and thus proof of the factory/gallery equivalence for the homodimer system.

V. HOMOPOLYMER IN EQUILIBRIUM

We now turn to a system that is far simpler to define in a rule-based manner than in a species-based manner. Consider a factory comprising two operators: (77) F1=∑iAˆia˜ib˜i,F2=∑i,jA‾iA‾jaˆibˆjJˆij.

This is similar to the homodimer factory, but F2 differs in that it forms an asymmetric (rather than symmetric) bond between two A particles. Specifically, the summand in F2 occupies a site ai on the Ai monomer, a site bj on the Aj monomer, and forms a bond Jij between them. F2 is not multiplied by a symmetry factor because Jij≠Jji. These factory operators are represented graphically in Fig. 5(a). We define a rule-based Hamiltonian for this system in the familiar way: (78) H=-μ∑iA‾i+ϵ∑i,jJ‾ij,

as represented in Fig. 5(b). The resulting gallery [Fig. 5(c)] is far more complex than that of the homodimer: it comprises creation operators for polymer chains and polymer rings of all lengths. Here, x-chains and x-rings are created by the operators (79) Cˆx=∑i1,…,ixAˆi1⋯Aˆixaˆi1⋯aˆix-1×bˆi2⋯bˆixJˆi1i2⋯Jˆix-1ix,

(80) Rˆx=1x∑i1,…,ixAˆi1⋯Aˆixaˆi1⋯aˆix×bˆi1⋯bˆixJˆi1i2⋯Jˆix-1ixJˆixi1.

These operators are more clearly expressed in diagrammatic notation [Fig. 5(d)]. Note the factor of 1/x in Eq. (80); this is needed to compensate for redundancy in the sum over internal indices that results from x-rings having rotational symmetry.

The equivalent species-based Hamiltonian has an infinite number of terms, each with its own bare chemical potential: (81) H≃-∑x=1∞μCxC‾x+μRxR‾x,

where μCx=xμ-(x-1)ϵ and μRx=xμ-xϵ are the bare chemical potentials for x-chains and x-rings. The corresponding number of complex-specific microstates are NCx=Nx and NRx=1xNx. Putting these together, we obtain the effective chemical potentials of each species: (82) μCx′=xμ-ϵ+ϵ-kTlogNxV,μRx′=xμ-ϵ-kTlogNxxV.

Unlike in the homopolymer system, it is not possible to renormalize μ and ϵ so that all μCx′ and μRx′ are independent of volume. Keeping μC1′ independent of N and V requires fixing the value of μ′=μ-kBTlogNV as in the monomer and homodimer systems. Keeping all other μCx′ independent of N and V then requires fixing ϵ′=ϵ+kBTlogV. The effective chemical potentials for all species thus become (83) μCx′=xμ′-ϵ′+ϵ′,μRx′=xμ′-ϵ′-kBTlogx-kBTlogV.

Constraining the concentrations of all Cx species independent of V therefore requires that the concentration of all Rx species scale as V-1. We note, however, that the converse is not possible, i.e., one cannot choose a definition for μ′ and ϵ′ so that the concentrations of all ring species are independent of V.

We therefore conclude that, for the homopolymer system to behave sensibly in the V→∞ limit, the concentrations of all chain polymers must remain fixed whereas the concentrations of all ring polymers must vanish. This makes sense: as V increases, the number of free ends with which the free end of an x-chain can interact increases in proportion to V. To preserve the concentrations of all x-chains, e-βϵ must scale as V-1. The probability of one free end of an x-chain interacting with the other free end of the same polymer will thus scale as V-1, and the concentration of all ring polymers will also scale as V-1.

Given Eq. (83), defining η=eβμ′-ϵ′ and computing (84) logZ=∑x=0∞eβμCx+∑x=0∞eβμRx,

we find that the log partition function density of the system is, for 0<η<1, (85) logZV=eβϵ′1-η-log1-ηV.

The left and right terms of the resulting expression are the respective contributions from chains and rings. The V-1 scaling of the second term reflects the vanishing of ring species as V→∞. The remaining chain contribution diverges as η→1 from below. In this limit, the concentration of A particles diverges as (86) ⟨A‾⟩V=1β∂∂μ′logZV≈eβϵ′δ2.

Since the concentration of each chain species Cx is given by eβμCx′=eβϵ′ηx, the distribution of chain lengths is distributed exponentially with decay rate logη, and flattens out as η→1 from below. Defining δ=1-η, the mean and variance of these chain lengths diverge as (87) x=η1-η≈1δ,varx=η(1-η)2≈1δ2.

To show the equivalence of the factory and gallery representations of the homopolymer system, we again invoke Wick’s theorem. In this case, however, the allowable contractions are more complex. Paralleling the analysis for the homodimer, we define the compound operators (88) Mˆi=Aˆia˜ib˜i,Aˆia=Aˆiaˆib˜i,

(89) Aˆib=Aˆia˜ibˆi,Aˆiab=Aˆiaˆibˆi,

(90) A‾ia=A‾iaˆib˜i,A‾ib=A‾ia˜ibˆi.

These six operators obey the contraction rules (91)

(92)

These rules are shown diagrammatically in Fig. 5(e). The reader may notice that the first two contraction products are not central, and thus violate an assumption of Wick’s theorem. We find, however, that Wick’s theorem still holds if we allow contraction products to participate in additional contractions. Fully contracted operator products therefore have all A‾ia and A‾ib operators participating in one contraction each, whereas each Aˆ operator may participate in zero, one, or two contractions.

We can use the above contraction rules to derive the complexes generated from any given term in ∣sum⟩=eF2eF1|0⟩. Fig. 5(f) shows the result for one such term. A generalized version of this computation is used in Appendix C to prove the factory/gallery equivalence for the homopolymer system.

VI. NONEQUILIBRIUM SYSTEMS

A. Species-based formalism

We now turn to the problem of determining the macrostate master equation [Eq. (10)] given a rule-based microstate master equation [Eq. (4)]. We start, however, by investigating how to specify and analyze a more standard species-based microstate master equation.

Species-based master equations are built from individual reactions, each of which annihilates a specified set of complexes and creates a new set in their place. Suppose there are H distinct species-specific reactions. The transition operator in Eq. (4) will then have the form (93) W=∑h=1HrhQh-Q`h,

where rh is the rate at which reaction h occurs, Qh is a “reaction operator” that effects this reaction when applied to a macrostate |n→⟩, and Q` is a corresponding “depletion operator.” Each reaction operator has the form (94) Q=∏k=1KGˆkpkpk!Gˇkqkqk!,

where q→=q1,…,qK is an abundance vector that describes the reactants and p→=p1,…,pK is a vector that describes the products. For the sake of simplicity we assume that p→ and q→ do not overlap (i.e., pk=0 and/or qk=0 for every k). Applying the conjugate of this operator to the macrostate, one finds that (95) Q†|n→⟩=Np→Ω(n→,q→,p→)|n→+q→-p→⟩,

where (96) Ω(n→,q→,p→)=∏knk+qk-pkqk1nk≥pk

is a coefficient that depends on n→, p→, and q→, but not otherwise on the details of the reaction, and (97) Np→=∏kNkpk≈∏kNkpkpk!

is the number of distinct product microstates created when Q is applied to a single reactant microstate.

In our formalism, the depletion operator O` corresponding to any reaction operator O is given in terms of the microstate-specific components via (98) O=∑ℐ,𝒥oℐ𝒥|𝒥⟩⟨ℐ|⇒O`=∑ℐ,𝒥oℐ𝒥|ℐ⟩⟨ℐ|,

where ℐ and 𝒥 index all microstates of the system and oℐ𝒥 is the rate at which |ℐ⟩ is transformed into |𝒥⟩. Note that, by this definition, all depletion operators are self-conjugate. In Appendix D we show that Eq. (98), together with the assumption of non-overlapping p→ and q→, leads to a depletion operator corresponding to Q of (99) Q`=∏kG˜kpkG‾kqk.

Applying Q`†=Q` to the macrostate then gives (100) Q`n→=Np→Ωn→,q→,q→n→,

Adding back the subscript h on p→ and q→ and using Eq. (11), we obtain an expression for the flux projector: (101) |J(n→)⟩=∑h=1HrhNp→hΩn→,q→h,p→hn→+q→h-p→h-Ωn→,q→h,q→h|n→⟩.

The difficulty with this species-based formulation is that, in systems that admit multi-particle complexes, the form of the transition operator in Eq. (93) does not reflect the underlying simplicity of the system. Rather, the rates rh and reaction operators Qh are derived quantities that follow from an (often much smaller) set of informally stated rules. Furthermore, even manually specifying the right-hand side of Eq. (93) can be tricky: a small number of rules can lead to a very large (or even infinite) number of reactions H, and both rh and Qh can depend in nontrivial ways on the elemental parameters that govern those rules. We now show how our formalism addresses this problem by enabling the rule-based definition of W.

B. Rule-based formalism

We specify the transition operator in a rule-based manner as follows. Suppose we have a system defined by L reaction rules. For each rule l, we specify a rate rl and a “reaction rule operator” Rl. The transition operator is then given by (102) W=∑l=1LrlRl-R`l.

where R`l is the depletion operator corresponding to Rl.

Suppose a rule operator R is able to drive M different species-specific reactions. For each reaction m, let q→m denote the number of reactant species and p→m the number of product species. We find that (103) R†|n→⟩=∑m=1MσmΩn→,q→m,p→m|n→+q→m-p→m⟩,

where σm is the number of distinct product microstates that can result from each reactant microstate in an m-type reaction. The depletion operator R` follows from the rule operator R using Eq. (98). Applying R`†=R` to |n→⟩, (104) R`|n→⟩=∑m=1MNp→mΩ(n→,q→m,q→m)|n→⟩.

Adding back the l indices, we obtain the flux projector (105) |J(n→)⟩=∑l=1L∑m=1MlrlσlmΩ(n→,q→lm,p→lm)n→+q→lm-p→lm-Ωn→,q→lm,q→lm|n→⟩.

C. Macroscopic master equation for the homopolymer

We now use our rule-based formalism to derive the macrostate master equation for the homopolymer system. Out of equilibrium, the dynamics of this system can be defined by L=4 reaction rule operators: (106) R1=∑iAˆia˜ib˜i,R3=∑i,jJˆijaˆibˆjA‾iA‾j,R2=∑iAˇia˜ib˜i,R4=∑i,jJˇijaˇibˇjA‾iA‾j.

Here, R1 creates a monomeric particle, R2=R1† destroys a monomeric particle, R3 creates an interaction between two particles, and R4=R3† destroys an interaction. Note that these four operators are given by the factory operators of Section V and their conjugates. The corresponding depletion operators are (107) R`1=∑iA˜ia˜ib˜i,R`3=∑i,jJ˜ija˜ib˜jA‾iA‾j,R`2=∑iA‾ia˜ib˜i,R`4=∑i,jJ‾ija‾ib‾jA‾iA‾j.

Diagrammatic representations of these rule operators and depletion operators are shown in Fig. 6(a).

In Section V we showed that the macrostates of the homopolymer system are given by (108) |c→,r→⟩=c1,r1,c2,r2,…=∏x=1∞Cˆxcxcx!Rˆxrxrx!0,

where cx and rx respectively indicate the number of chains and rings of length x. To compute the macroscopic master equation, we compute the flux projector (109) |J(c→,r→)⟩=∑l=14rl(Rl†-R`l)|c→,r→⟩.

We now evaluate Eq. (109) term-by-term. To ease notation, we show in the macrostate only the elements of c→ and r→ that change upon application of each operator and denote the unchanged elements by “...”. The l=1 term corresponds to monomer creation. As illustrated in Fig. 6(b), R1 maps a single microstate (corresponding to no reactants) to σ1≈N microstates (corresponding to all possible monomer states). By Eq. (103) and Eq. (104), (110) (R1†-R`1)|c→,r→⟩≈Nc1-1,…-….

R2 maps a single monomeric particle microstate to σ2=1 microstate (i.e., no products), and so (111) (R2†-R`2)|c→,r→⟩=c1+1c1+1,…-c1….

The effects of R3 and R4 are more complex. R3 can effect three different types of reactions depending on the reactants. First, R3 can join together an x-chain and y chain (x<y) to get an (x+y)-chain; this can be done in σ3=2 different ways [Fig. 6(c)]. Second, R3 can can join together two x-chains to get a 2x-chain; this can be done in σ3=2 ways [Fig. 6(d)]. Third, R3 can join together the ends of an x-chain to get an x-ring; this can be done in only σ3=1 way [Fig. 6(e)]. We therefore find that (112) (R3†-R`3)|c→,r→⟩=∑x<y2cx+1cy+1cx+1,cy+1,cx+y-1,…-cxcy|…⟩+∑x2cx+2cx+12cx+2,c2x-1,…-cxcx-12|…⟩+∑xcx+1cx+1,rx-1,…-cx|…⟩.

Similarly, the inverse operator R4 can separate an (x+y)-chain into an x-chain and y-chain (x<y) in σ4=2 ways, can separate a 2x-chain into two x-chains in only σ4=1 way, and can cut an x-ring to get an x-chain in σ4=x different ways. Consequently, (113) (R4†-R`4)|c→,r→⟩=∑x<y2cx+y+1cx-1,cy-1,cx+y+1,…-cx+y|…⟩+∑xc2x+1cx-2,c2x+1,…-c2x|…⟩+∑xxrx+1cx-1,rx+1,…⋅-rx|…⟩.

The depletion terms in these expressions can be simplified as follows (114) ∑x<y2cxcy+∑xcxcx-1+∑xcx=nchain2,∑x<y2cx+y+∑xc2x+∑xxrx=nlink,

where nchain=∑xcx is the total number of chains and nlink=∑x(x-1)cx+xrx is the number of links among all chains and rings.

Evaluating the inner product ⟨J(n→)|ψ(t)⟩, we thus obtain the macroscopic master equation for the homopolymer: (115) ddtP(…)=Nr1Pc1-1,…+r2c1+1Pc1+1,…+r3[∑x<y2cx+1cy+1Pcx+1,cy+1,cx+y-1,…+∑xcx+2cx+1Pcx+2,c2x-1,…+∑xcx+1Pcx+1,rx-1,…]+r4[∑x<y2cx+y+1Pcx-1,cy-1,cx+y+1,…+∑xc2x+1Pcx-2,c2x+1,…+∑xxrx+1Pcx-1,rx+1,…]-Nr1+r2c1+r3nchain2+r4nlinkP….

To verify this result we focus on the depletion term, i.e., the coefficient of P(c→,r→). The overall monomer creation rate is Nr1. This makes sense, as r1 is the per-mode rate of excitation of the A field. r1 must therefore scale as V/N. The r2 term reflects our assumption that only A particles engaging in no interactions are able to be annihilated; it does not scale with V or N. The r3 term reflects the fact that new interactions form between the a-end of one chain and the b-end of either another chain or the same chain r3 scales as V-1. The total number of available reactants in the system is therefore nchain2. The r4 term reflects the assumption that any link can be annihilated regardless of the complex in which it occurs r4 does not scale with V or N).

VII. SIMULATIONS

A. Deterministic simulations

Our formalism enables the numerical solution of the microstate master equation – at least when N is sufficiently small. One first defines a set of modes ℳ, as well as a set of mode-specific rules ℛ=Rli,rl, where l indexes the qualitatively different rules as in Eq. (102), and i indexes the internal states of the particles that each rule acts upon (so that Rl=∑iRli). The transition operator is then computed using (116) W=∑l=1Lrl∑iRli-R`li.

Given an initial state vector |ψ(0)⟩, the state at time t is computed using (117) ψt=exptWψ0.

Figs. 7(a)–7(c) show the results of this computation for three example systems: a monomer system (N=20), a homodimer system (N=4), and a homopolymer system (N=3).

The primary limitation of this approach is that it requires computations involving very large vectors and matrices: |ψ⟩ is a 2|ℳ|-dimensional vector, W is a 2|ℳ|×2|ℳ| matrix, and |ℳ| is polynomial in N, e.g., |ℳ|=N for the monomer, |ℳ|=2N+N2 for the homodimer, and |ℳ|=3N+N2 for the homopolymer. Even using sparse matrix methods, we have found the direct evaluation of Eq. (117) to be impractical for all but very small values of N. Nevertheless, these deterministic simulations provide a valuable way to check the accuracy of the stochastic simulations that we now turn to.

B. Stochastic simulations

Our formalism also enables stochastic simulations using the Gillespie algorithm [23, 24]. Importantly, these stochastic simulations can be carried out using much larger values of N. Algorithm 1 is one algorithm that does this. After explaining how the algorithm works, we illustrate its operator stepping through one iteration of the algorithm for a homodimer system. We then present computational results obtained using this algorithm and discuss the algorithm’s current limitations.

The input to Algorithm 1 consists of two strings, strans and sinit. The string strans specifies the set of rules and their corresponding rates, while sinit specifies the initial state of the system, s0. The output of the algorithm is a trajectory object 𝒯, the downstream parsing of which provides time traces for the abundances of all single particles and complexes. Post-hoc processing of 𝒯 then provides time traces for all species of complex. We emphasize that this algorithm only tracks the excitation states of modes. In particular, there is no need during the execution of the algorithm to enumerate the different possible species of complex or even to track which species occur. Rather, time traces for the abundances of different complexes are determined only during the post-processing of 𝒯.

After processing its inputs, the algorithm sets time t to zero and the trajectory object 𝒯 to be the empty set. Next, the function Initialize takes strans and sinit as inputs and outputs four sets of objects:

𝒪 is the set of all mode-specific operators (i.e., creation, annihilation, presence, and absence operators) for all fields. Every O∈𝒪 has the attributes O.mode, O.rules, and O.eligible. O.mode is a reference to the operator’s mode M∈ℳ. O.rules is a set of references to the rules R∈ℛ that include O in their operator product. O.eligible is a Boolean flag indicating whether O can be applied to st without killing it.

ℳ is the set of all modes for all fields. Every M∈ℳ has the attributes M.operators and M.excited. M.operators is a set of references to the four mode-specific operators (i.e., Mˆ, Mˇ, M‾, and M˜). M.excited is a Boolean value indicating whether the mode is excited. If TRUE, Mˆ and M˜ are ineligible, while Mˇ and M‾ are eligible. If FALSE, Mˆ and M˜ are eligible, whereas Mˇ and M‾ are ineligible.

ℛ is the set of all mode-specific rules, each defined as a product of mode-specific operators. Every rule R∈ℛ has the attributes R.operators and R.rate. R.operators is a set of references to the mode-specific operators O∈𝒪 that comprise R. R.rate is the rate at which the rule is applied when eligible.

𝒞 is the “state constructor,” i.e., the set of creation operators that, when applied to |0⟩, yield st, the system state at time t. The specific 𝒞 returned by Initialize corresponds to s0, which is specified by the string sinit. In particular, if the user specifies that s0=|0⟩, then 𝒞={}.

𝒪, ℳ, and ℛ are static, whereas 𝒞 evolves in time. Also note that (118) M=O.mode⇔O∈M.operators,

(119) R∈O.rules⇔O∈R.operators

for all O∈𝒪, M∈ℳ, and R∈ℛ.

After initialization, Algorithm 1 enters a while loop. Each execution of the loop applies one rule R to the system state st, then steps the system forward in time from t to t+Δt. The contents of this while loop are as follows:

Line 5 identifies the set of eligible rules and saves them in set ℛ*. This is the most computationally expensive part of Algorithm 1, since it must loop through every possible rule and test that rule for eligibility using the function IsRuleEligible. A rule R is eligible if and only if every operator in R.operators is eligible.

Line 6 uses the rates of the rules in ℛ* to randomly sample a time step Δt and a corresponding rule R via the Gillespie algorithm, which is implemented by GillespieStep. Time t is then incremented by Δt.

Lines 8–14 update the state constructor 𝒞 and the excitation state of modes affected by the rule R. To do this, a for loop is carried out over all operators O∈R.operators. If O is a creation operator, then O is added to the state constructor 𝒞. Alternatively, if O is an annihilation operator, then the corresponding creation operator O† is removed from the state constructor 𝒞. In either case, FlipModeExcitation is applied to the mode M=O.mode. This flips the excitation attribute of M, as well as the eligible attribute of all operators in M.operators.

Finally, the loop stores a tuple reporting the updated time t, the rule R, and the resulting state constructor 𝒞 in the trajectory object 𝒯 .

We now illustrate how this algorithm works by following it through initialization and one execution of the while loop. Assume that the string strans specifies the following kinetic rules for a homodimer system:

(120) R1=∑iAˆia˜i,R3=∑i<jA‾iA‾jaˆiaˆjIˆij,R2=∑iAˇia˜i,R4=∑i<jA‾iA‾jaˇiaˇjIˇij.

We further suppose that the string sinit specifies an initial state containing four monomers having indices 2, 3, 5, and 7, i.e., (121) s0=Aˆ2Aˆ3Aˆ5Aˆ70.

Taking strans and sinit as input, the function Initialize returns the following sets of objects: (122) ℳ=Aii∪aii∪Iiji<j,𝒪=Aˆi,Aˇi,A‾i,A˜ii∪aˆi,aˇi,a‾i,a˜ii∪Iˆij,Iˇij,I‾ij,I˜iji<j,ℛ=Aˆia˜ii∪{Aˇia˜i}i∪A‾iA‾jaˆiaˆjIˆiji<j∪A‾iA‾jaˇiaˇjIˇiji<j,𝒞=Aˆ2,Aˆ3,Aˆ5,Aˆ7.

Here the indices i and j are understood to run over 1,...,N, and the rates rl corresponding to each Rli are kept implicit. Initialize also sets the values of M.excited for all M∈ℳ, and of O.eligible for all O∈𝒪. M.excited is TRUE for modes A2, A3, A5, and A7, and is FALSE for all other modes (including all modes of the fields a and I). Consequently, the set 𝒪* of eligible operators is (123) 𝒪*=Aˆi,A˜ii∉{2,3,5,7}∪{Aˇi,A‾i}i∈{2,3,5,7}∪aˆi,a˜ii∪Iˆij,I˜iji<j.

Now consider the first execution of the while loop. In line 5 of Algorithm 1, IsRuleEligible is evaluated on every rule R∈ℛ. Based on the eligibility of each operator in R.operators, the eligible rules are found to be (124) ℛ*=Aˆi,a˜ii∉{2,3,5,7}∪{Aˇi,a‾i}i∈{2,3,5,7}∪A‾iA‾jaˆiaˆjIˆiji<j∈{2,3,5,7}.

Next, GillespieStep chooses a random eligible rule R∈ℛ* and a time increment Δt that is then added to t. Suppose (125) R=A‾3A‾7aˆ3aˆ7Iˆ37

is chosen. Applying R to the initial state vector then yields the updated state vector, (126) sΔt=Rs0=Aˆ2Aˆ3Aˆ5Aˆ7aˆ3aˆ7Iˆ370.

To register this change, the state constructor is updated to (127) 𝒞=Aˆ2,Aˆ3,Aˆ5,Aˆ7,aˆ3,aˆ7,Iˆ37.

Calls to FlipModeExcitation then set the excited attributes of modes a3,a7,I37 to TRUE and flip the eligible attribute of the four operators corresponding to each of these modes. The resulting set of eligible operators is (128) 𝒪*=Aˆi,A˜ii∉{2,3,5,7}∪{Aˇi,A‾i}i∈{2,3,5,7}∪aˆi,a˜ii∉{3,7}∪aˇi,a‾ii∈{3,7}∪Iˆ25,I˜25,Iˇ37,I‾37,.

Finally the updated time t, the chosen rule R, and the updated state constructor 𝒞 are added as a tuple to the trajectory 𝒯 .

In the next execution of the while loop, the set ℛ* of eligible rules is found by IsRuleEligible to be (129) ℛ*=Aˆia˜ii∉{2,3,5,7}∪{Aˇia˜i}i∈{2,5}∪A‾2A‾5aˆ2aˆ5Iˆ25,A‾3A‾7aˇ3aˇ7Iˇ37.

It is worth noting the changes to ℛ* vs. Eq. (124). The monomer annihilation rules Aˇ3a˜3 and Aˇ7a˜7 have been removed because the modes a3 and a7 are now excited, making the operators a˜3 and a˜7 ineligible. This prevents monomers joined by interactions from being annihilated, thereby leaving “dangling” interactions. In addition, all interaction creation rules except A‾2A‾5aˆ2aˆ5Iˆ25 have become ineligible. This prevents the A3 and A7 monomers, which are already interacting with each other, from participating in multiple interactions. Finally, the interaction annihilation rule A‾3A‾7aˇ3aˇ7Iˇ37 becomes eligible, allowing the newly-formed bond to dissociate.

Fig. 7 shows this stochastic algorithm applied to monomer, homodimer, and homopolymer systems. Panels (a-c) validate this algorithm by showing that the mean abundance of each species found across 500 simulations closely traces the mean abundance predicted by the deterministic algorithm. As noted above, however, these comparisons can only be carried out at small N due to limitations of the deterministic algorithm. This stochastic algorithm can be performed at much larger values of N, e.g., Figs. 7(d)–7(f) show such simulations performed using N=100.

The size of N is still a limitation for Algorithm 1. The bottleneck is line 5, which requires iterating through all mode-specific rules in ℛ and testing each one for eligibility. This step takes o(|ℛ|) time, and |ℛ| is polynomial in N:|ℛ|=2N for the monomer, ℛ=2N+2N2 for the homodimer, and |ℛ|=2N+2N2 for the homopolymer. That said, Algorithm 1 was developed only as proof-of-principle and has not been optimized for efficiency. Indeed, we expect the bottleneck can be eliminated by using more sophisticated methods for tracking which operators and rules are eligible given the state constructor, thereby enabling simulations using arbitrarily large values of N.

VIII. EXPRESSIVENESS

Our formalism provides a rule-based approach for defining and analyzing a diverse array of stochastic chemical systems in which multi-particle complexes can form. Here we illustrate the expressiveness of our formalism by briefly considering a variety of such systems, both in and out of equilibrium.

In equilibrium, systems are defined by a factory F→ and a Hamiltonian H. Putting the Hamiltonian aside for the moment, it is interesting to consider the qualitatively different sets of complexes that can arise from simple changes to the factory. This expressiveness is perhaps most apparent in polymer systems. Fig. 8 shows nine such polymeric systems derived from variations on the homodimer and homopolymer systems described in previous sections. Two such systems, the isotropic homopolymer [Fig. 8(b)] and branched homopolymer [Fig. 8(c)] are further analyzed below.

Out of equilibrium, systems are defined by a set of rules, with each rule Rl assigned a corresponding rate rl. Fig. 9 shows five such nonequilibrium systems. These rules again derive from variations on the homodimer and homopolymer systems. Fig. 9 also shows the results of stochastic simulations carried out using Algorithm 1 for specific choices of the rate parameters.

A. Isotropic homopolymer

We now analyze the isotropic homopolymer shown in Fig. 8(b). This system is defined by one species of monomeric subunit (A) with two sites (a and b) capable of participating in two classes of symmetric interaction (I and J). We aim to compute the partition function for the system assuming the Hamiltonian (130) H=-μ∑iA‾+ϵ2∑i,jI‾ij+J‾ij.

As with the homopolymer of Section V, the two factory operators for this system generate chains and rings. However, the use of two distinct types of symmetric interaction complicates these species. First, there are no self-interacting monomers. Second, the I and J interactions alternate in all multimers. For x-chains this means that, when x is even, there are two different species related by the exchange of I and J interactions. When x is odd, there is instead only one-species of x-chain, as exchanging I and J is equivalent to flipping the order of the subunit indices. Moreover, x-rings occur only when x is even. These rings have x/2-fold rotational symmetry, and for x≥4 they additionally have 2-fold mirror symmetry.

The log partition function density of the system is therefore given by (131) logZV=∑x=1∞logZ(2x-1)-chain+∑x=1∞logZ2x-chain+logZ2-ring+∑x=2∞logZ2x-ring.

In terms of the effective chemical potential μ′=μ+kBTlogNV and effective interaction energy ϵ′=ϵ-kBTlogV, as well as the control parameter η=eβμ′-ϵ′ (which represents the Boltzmann weight of a single particle with a dangling bond), the single-complex partition functions in Eq. (131) are (132) logZ(2x-1)-chain=eβϵ′η2x-1,logZ2x-chain=2eβϵ′η2x,logZ2-ring=η22V,logZ2x-ring=η2x4xV.

Summing the terms in Eq. (131) we find that, for 0<η<1, (133) logZV=eβϵ′η1-η+eβϵ′η21-η2+η24V-log1-η24V.

As with the homopolymer, the partition function diverges as η→1 from below. Defining δ=1-η, we find that the concentration of A particles diverges in this limit, scaling as (134) ⟨A‾⟩V≈3eβϵ′2δ2.

Again, this divergence is dominated by the x-chains even at finite V. The factor of 32 difference between this result and the homopolymer result in Eq. (86) reflects the fact that, while both even and odd chains contribute to Eq. (134), there are twice as many species of even chains (but the same number of odd chains) as in the homopolymer system.

B. Branched directed homopolymer

Consider now the branched directed homopolymer system shown in Fig. 8(c). This system is defined by one species of monomeric subunit (A) with three sites (a,b, and c) capable of participating in two classes of directed interaction (I and J). We assume that the system is in thermal equilibrium and is governed by the Hamiltonian (135) H=-μ∑iA‾+ϵ∑i,jI‾ij+J‾ij.

The two factory operators for this system generate two broad classes of complex: “trees,” which branch out from a single A in which site c is unoccupied, and “groves,” which consist of multiple trees branching out from a central closed ring. The log partition function density of this system is therefore given by (136) logZV=Ztree+∑x=1∞Zx-grove,

where Ztree is the partition function for all trees and Zx-grove is the partition function for all groves with trees extending off a central x-ring. To proceed, we let ξ represent the partition function for a tree with a dangling interaction extending off the root. In terms of this quantity, (137) Ztree=eβϵ′ξ,Zx-grove=2xηx(1+ξ)xVx,

where η=eβμ′-ϵ′ is again the Boltzmann weight of a single particle with a dangling bond. In Ztree, the factor of eβϵ′ removes the effect of the dangling bond. In Zx-grove, the factor of ηx accounts for the particles and bonds in each x-ring, the factor of 2x accounts for the fact that each bond in the x-ring can be of type I or J, and the (1+ξ)x factor accounts for the fact that each A within the ring can either be bare or have a tree attached. As with the homopolymer, the factor of 1x compensates for the rotational symmetry of the ring while 1V reflects the entropic cost of self-circularization.

We solve for ξ by noting that the self-similar structure of each tree complex yields the recursion relation (138) ξ=η1+2ξ+ξ2.

Solving this quadratic equation and using the limiting behavior ξ≈η as η→0, we derive (139) ξ=1-2η-1-4η2η.

Expressing Ztree and Zx-grove in terms of ξ and summing Zx-grove over all x, we find that for 0<η<14, (140) logZV=eβϵ′1-2η-1-4η2η-log1-4η2V.

As with the homopolymer system, the contribution from the circularized (i.e., x-grove) species vanishes in the V→∞ limit. We also find that the partition function is not defined for η>14. In particular, the mean concentration of A particles is found to diverge as δ=14-η→0 according to (141) ⟨A‾⟩V≈eβϵ′2δ+18Vδ,

with the first and second terms respectively corresponding to trees and groves. This result is qualitatively different from the corresponding results for the directed homopolymer in Section V and for the isotropic homopolymer analyzed above. In particular, the circularized species (the groves) dominate over the linear species (the trees) in the δ→0 limit when V is kept finite. The asymptotic behavior of the system therefore depends on which limit one takes first, η→14 or V→∞. Moreover, the divergence is more mild than in the non-branched systems, scaling as either δ-1/2 or δ-1 (depending on how one handles V) rather than δ-2.

IX. DISCUSSION

We have introduced an algebraic formalism for the rule-based modeling of multi-particle complexes in stochastic chemical systems. This algebra is based on a Fock space that allows not only the creation and annihilation of particles, but also the joining of particles into complexes based on specified rules. The Fock space comprises three types of hard-core boson fields; these represent particles, particle-particle interactions, and occupied binding sites. We have also described a formal diagrammatic approach that facilitates the use of this algebra.

For equilibrium systems, we showed that the set of all possible complexes can be rigorously specified by a “factory” and a Hamiltonian. The factory is an ordered set of operators that define how to construct complexes; the Hamiltonian is an operator that specifies rules for computing the Gibbs free energy of a complex based on its components. We showed for multiple systems how these rule-based definitions can be used to compute generating functions and partition functions, and to analyze scaling behavior near critical polymerization concentrations.

For nonequilibrium systems, we showed how to rigorously specify system dynamics in a rule-based manner. Specifically, we showed how a set of reaction rule operators and corresponding rates can be used to define the transition matrix of a “microstate” master equation. From this transition matrix one can then analytically compute the corresponding “macrostate” master equation, which governs the time evolution of observables. We also developed a Gillespie algorithm for simulating stochastic chemical systems based on these rule operators and corresponding rates.

The essential feature of our formalism, one that distinguishes it from previous approaches for modeling many-body systems of classical particles, is that it explicitly represents internal particle states. These internal states endow each particle with its own identity, thus allowing preexisting particles to join together into multi-particle complexes. Notably, our approach to modeling these internal states is consistent with the behavior of quantum systems in the decoherence limit. In this limit, the reduced density matrix for each particle becomes diagonal with elements along the diagonal quantifying the probability of each energy eigenstate [25, 26]. The orthonormal microstates in our formalism correspond to these diagonal positions in the reduced density matrix (i.e., the energy eigenstates), and the probabilities that multiply these microstates correspond to the values of the reduced density matrix at these positions. Classical particles, such as proteins, have many distinct energy eigenstates corresponding to different internal excitation modes, and the internal states of particles in our formalism stand in for these modes. Our formalism assumes a specific number N of such modes, but this choice does not affect results when N is sufficiently large.

We envision a variety of potential analytic applications for our formalism. As in the work of Doi [2, 3], it may be possible to use the algebra we have introduced to carry out diagrammatic perturbation theory calculations. As in the work of Peliti [6] and Goldenfeld [5], it may also be possible to identify a path integral formulation of this algebra. While our analysis focused on zero-dimensional (i.e., well-mixed) systems, we expect that it should be straightforward to apply our formalism to spatially extended systems in which diffusion plays an important role (as in [2, 3, 5, 6]).

We also envision a variety of computational applications. Serious computational applications will require reworking our proof-of-principle algorithm (Algorithm 1) so that its speed does not scale with N, but we expect this will be straightforward using more advanced bookkeeping methods. The resulting algorithm may provide advantages over existing rule-based modeling approaches, since the underlying objects in our formalism (hard core bosons) are simpler than those of existing algorithms (e.g., port graphs). Our formalism might also facilitate the development of qualitatively different computational strategies for rule-based modeling, e.g., the use of finite state projections [27] or tensor networks [28] to approximate the solution of the master equation.

In this paper we have focused on polymer systems, which best illustrate how large complexes can arise from simple interaction rules. Our initial motivation in pursuing this project, however, was in biological systems. We specifically sought to develop methods for modeling the biophysical mechanisms of gene regulation. Gene regulation is controlled by large protein-nucleic acid complexes [29, 30], with individual regulatory sequences able to nucleate the formation of large numbers of different macromolecular assemblies. A major goal in this field is to understand these complexes as well as their effects on gene expression using biophysical models, and substantial progress has been made using both equilibrium models [31–34] and nonequilibrium models [35–37]. In particular, biophysical modeling provides a principled approach to deciphering how gene regulatory programs are encoded in DNA and RNA sequences [38–42]. Our work provides a formal language in which such biophysical models can be expressed and then analyzed. This capability may allow researchers to systematically explore the space of biophysical models of gene regulation, thereby automating the construction and inference of biophysical models for gene regulatory codes.

ACKNOWLEDGMENTS

We thank Rob Phillips for his support and encouragement throughout this project, as well as Muir Morrison for his early work with JBK on this topic. RJR further thanks Sergei Gukov and Alexei Kitaev for helpful discussions. The work of JBK was supported by NIH grants GM133777 and HG011787. The work of RJR was supported by NIH grant GM118043. This research was performed in part at Aspen Center for Physics, which is supported by NSF grant PHY-2210452.

Appendix A: Algebra of mode and field operators

Since each mode represents a hard-core boson, mode-specific creation and annihilation operators are nilpotent, i.e., (A1) Aˆi2=Aˇi2=0.

When multiplied by presence and absence operators for the same mode, the creation and annihilation operators are readily seen to satisfy (A2) A‾iAˆi=AˆiA˜i=Aˆi,A˜iAˆi=AˆiA‾i=0,

(A3) A˜iAˇi=AˇiA‾i=Aˇi,A‾iAˇi=AˇiA˜i=0.

Presence and absence operators are both idempotent, i.e., (A4) A‾i2=A‾i,A˜i2=A˜i,

and mixed products of presence and absence operators for the same mode vanish: (A5) A‾iA˜i=A˜iA‾i=0.

From these properties and the fact that operators for distinct modes commute, we get the following commutation relations for mode operators: (A6) [Aˇi,Aˆj]=δij(A˜i-A‾i)=δij1-2A‾i,

(A7) [A‾i,Aˆj]=[Aˆi,A˜j]=δijAˆi,

(A8) [A˜i,Aˇj]=[Aˇi,A‾j]=δijAˇi.

Summing over indices gives the corresponding commutation relations for field operators: (A9) [Aˇ,Aˆ]=A˜-A‾=N-2A‾,

(A10) [A‾,Aˆ]=[Aˆ,A˜]=Aˆ,

(A11) [A˜,Aˇ]=[Aˇ,A‾]=Aˇ.

These commutation relations, together with A‾|0⟩=0, allow us to compute the impact of each field operator on the macrostate |n⟩, (A12) Aˆ|n⟩=AˆAˆnn!|0⟩=(n+1)Aˆn+1(n+1)!|0⟩=n+1n+1,

(A13) A‾|n〉=A‾A^nn!|0〉=1n![A‾,Aˆn]|0〉=1n!∑k=0n−1Aˆn−k−1[A‾,Aˆ]Aˆk|0〉=nAˆnn!|0〉=n|n〉,

(A14) Aˇ|n〉=AˇAˆnn!|0〉=[Aˇ,Aˆnn!]|0〉=1n!∑k=0n−1Aˆn−k−1[Aˇ,Aˆ]Aˆk|0〉=1n!∑k=0n−1Aˆn−k−1(N−2A‾)Aˆk|0〉=1n!∑k=0n−1(N−2k)Aˆn−1|0〉=1n!(nN−2n(n−1)2)Aˆn−1|0〉=(N−n+1)Aˆn−1(n−1)!|0〉=(N−n+1)|n−1〉,

(A15) A˜|n⟩=A˜Aˆnn!|0⟩=(N-A‾)Aˆnn!|0⟩=(N-n)Aˆnn!|0⟩=N-nn.

Appendix B: Non-commutation of monomer and dimer operators

Here we derive the expression for [Mˆ,Dˇ] in Eq. (45). We begin by evaluating the commutator on the individual composite mode operators: (B1) [Mˆk,Dˇij]=IˇijAˆka˜k,AˇiAˇjaˇiaˇj=δkiIˇijAˇjaˇjAˆka˜k,Aˇiaˇi+δkjIˇijAˇiaˇiAˆka˜k,Aˇjaˇj.

Considering the commutator in the k=i term and dropping the subscripts for brevity, we find that (B2) [Aˆa˜,Aˇaˇ]=Aˆ[a˜,Aˇ]aˇ+[Aˆ,Aˇ]a˜aˇ+AˇAˆ[a˜,aˇ]+Aˇ[Aˆ,a˜]aˇ=A‾aˇ=Aˆa˜Aˇaˇ.

The same holds for the k=j term. Substituting these back into Eq. (B1) and summing over k gives (B3) ∑k[Mˆk,Dˇij]=∑kδkiAˆia˜iIˇijAˇiaˇiAˇjaˇj+∑kδkjAˆja˜jIˇijAˇiaˇiAˇjaˇj=(Mˆi+Mˆj)Dˇij.

From this we recover Eq. (45).

Appendix C: Factory/gallery equivalence for the homopolymer

We now prove the equivalence of the factory and gallery representations for the homopolymer system in Section V. As with the homodimer, we do this by evaluating individual terms in the Taylor expansion of eF2eF1|0⟩, with factory operators defined as in Fig. 5(a). By inspection we see that all full contractions of F2pF1q must consist of a product of x-chain (Cˆx) and x-ring (Rˆx) operators. Wick’s theorem therefore gives (C1) F2pF1q0=∑cx,rx∣p,qΩcx,rxp,q∏x=1∞CˆxcxRˆxrx0,

where cx denotes the number of x-chains, rx denotes the number of x-rings, cx,rx∣p,q denotes all sets of these numbers that are consistent with p bonds and q particles, i.e., which satisfy (C2) q=∑x=1∞xcx+xrx,p=∑x=1∞[x-1]cx+xrx,

and Ωcx,rx(p,q) is a combinatorial coefficient that quantifies the number of distinct contractions that yield cx,rx∣p,q.

We now compute Ωcx,rx(p,q). The number of ways to partition q monomers among cx distinct x-chains and rx distinct x-rings is given by (C3) ωq=q!∏x=1∞(x!)cx+rx.

Similarly, the number of ways to partition p bonds among the x-chains and x-rings is (C4) θp=q!∏x=1∞([x-1]!)cx(x!)rx.

Since one can rearrange the x-chains among themselves and the x-rings among themselves without changing the result, the number of unique partitions is the product of ωq and θq divided by an exchange factor of πp,q=∏x=1∞cx!rx!. Moreover, there are σx=(x-1)!x! distinct ways of constructing each x-chain from a given set of x particles and x-1 bonds, and ρx=x!x! ways to construct each x-ring from a set of x particles and x bonds. Note that the circular symmetry, which contributes a factor of 1/x to this second quantity, is already accounted for in the definition of the x-ring and should not be doublecounted here. We therefore find that (C5) Ωcx,rx(p,q)=ωqθpπp,q∏x=1∞σxcxρxrx=p!q!∏x=1∞cx!rx!.

Consequently, (C6) eF2eF2|0⟩=∑p,q1p!q!∑cx,rx∣p,qp!q!∏xcx!rx!∏xCˆxcxRˆxrx|0⟩=∑c1,c2,…∑r1,r2,…∏xCˆxcxcx!Rˆxrxrx!|0⟩=e∑x(Cˆx+Rˆx)0.

This establishes the factory/gallery equivalence for the homopolymer.

Appendix D: Derivation for the species-specific depletion operator

Here we derive the species-specific depletion operator in Eq. (99). First we express the species-specific reaction operator as (D1) Q=∏k=1KPkQk,wherePk=Gˆkpkpk!,Qk=Gˇkqkqk!.

The assumption of non-overlapping p→ and q→ implies that Pk=1 and/or Qk=1 for all k. Because of this, (D2) Q`=∏k=1KP`kQ`k.

Note that this will generally not be true if any Pk and Qk are both non-unity, since the depletion version of a product of operators for a given field is generally not the product of the depletion version of each operator.

Next we express Pk and Qk in terms of mode-specific operators. Using the fact that the mode-specific creation and annihilation operators are nilpotent (and dropping k to ease notation), (D3) P=∑ℐ:|ℐ|=p∏i∈ℐGˆi,Q=∑ℐ:|ℐ|=q∏i∈ℐGˇi.

Transforming Gˆi→G˜i and Gˇi→G‾i, we obtain the corresponding depletion operators (D4) P`=∑ℐ:|ℐ|=p∏i∈ℐG˜i,Q`=∑ℐ:|ℐ|=q∏i∈ℐG‾i.

Applying each of these operators to a macrostate |n⟩, one finds that (D5) P`n=N-npn=G˜pn,

(D6) Q`n=nqn=G‾qn,

and therefore, (D7) P`=G˜p,Q`=G‾q.

Reintroducing the k subscripts, we obtain Eq. (99): (D8) Q`=∏k=1KP`kQ`k=∏k=1KG˜kpkG‾kqk.

Appendix E: Code availability

Python code implementing Algorithm 1 as well as the Jupyter Notebooks used to perform the simulations in Fig. 7 and Fig. 9 are available at github.com/Rebecca-J-Rousseau/RousseauKinney2024_algebra.

FIG. 1. Homopolymer in zero dimensions. The system comprises a single species of component particle, with each particle having two domains capable of forming heterotypic interactions. (a) In thermal equilibrium, system behavior is governed by the chemical potential (μ) and interaction energy (ϵ). (b) Out of equilibrium, the system is governed by the kinetic rates for monomer creation r1, monomer annihilation r2, interaction formation r3, and interaction dissolution (r4. (c) The rules in panel (a) generate an infinite number of possible complexes: x-chains and x-rings for all x=1,2,…, each complex with a distinct chemical potential. The log terms in the chemical potentials of the x-rings result from their rotational symmetry. (d) The rules in panel (b) result in an infinite web of reactions between different polymeric complexes, with reaction rates proportional to r3 and r4 in ways that depend nontrivially on the identities of the specific reactants and products.

FIG. 2. Diagrammatic notation for operator products and sums thereof over internal indices. (a) Examples of simple and compound mode operators. (b) Examples of simple and compound field operators. (c) Wick contractions relevant to the homodimer system.

FIG. 3. Diagrams facilitate algebraic calculations. (a) Evaluation of the p=2, q=5 term of Eq. (63) in terms of normally ordered operator products. (b) Result of the computation in panel (a) applied to the vacuum state.

FIG. 4. Diagrammatic proof of the factory/gallery equivalence for the homodimer system.

FIG. 5. Diagrammatic specification of the homopolymer system. (a) The factory. (b) The Hamiltonian. (c) The resulting gallery. (d) The Cˆx and Rˆx comprising the gallery, and (e) the four contraction rules for the system. For conciseness, panels (c-e) hide all operator names except for the decorators on the A operators. (f) An example algebraic calculation carried out using the commutation rules in panel (e).

FIG. 6. The eight types of species-specific reactions that occur in the nonequilibrium homopolymer system. (a) Diagrammatic representations of the four reaction rule operators and their corresponding depletion operators. (b) R1 maps the vacuum state to N different monomer microstates, while R2 maps each monomer microstate to a single vacuum state. (c-e) The three types of species-specific reactions effected by R3 and R4. (c) R3 can link an x-chain and y-chain together in two different ways, while R4 can split an (x+y)-chain into an x-chain and y-chain in two different ways. (d) R3 can link two x-chains together in two different ways, while R4 can split a 2x-chain into two x-chains in only one way.(e) R3 can circularize an x-chain in only one way, while R4 can linearize an x-ring in x different ways. Nodes are shown as open dots in panels (b-e) to indicate specific internal states; the numbers below each node indicate example values for the internal index i of each component particle.

FIG. 7. Simulations of nonequilibrium systems. (a-c) Deterministic and stochastic simulations for (a) a monomer system (N=20), (b) a homodimer system (N=4), and (c) a homopolymer system (N=3). (d-f) Stochastic simulations for the same systems as in panels (a-c) using N=100. Black dashed lines indicate mean abundances from deterministic simulations. Solid lines plot mean abundances from 500 stochastic simulations. Error bands indicate standard deviations in abundance across the stochastic simulations.

FIG. 8. Factories defining various polymer systems. Each panel shows the factory operators used to define system in equilibrium. The homopolymer system analyzed in Section V is shown in panel (a) for comparison.

FIG. 9. Stochastic simulations of various nonequilibrium models. Shown are results for (a) a heterodimer system, (b) an occlusive binding system, (c) a cooperative binding system, (d) a heteropolymer system, and (e) a branched homopolymer system. For each system, rl+ denotes the rate at which the forward rule, Rl, is applied, while rl- denotes the rate at which the reverse rule, Rl†, is applied. Solid lines indicate mean abundances from 500 stochastic simulations using N=100. Error bands indicate abundance standard deviations across the simulations.
==== Refs
[1] Hlavacek W. , Faeder J. , Blinov M. , Posner R. , Hucka M. , and Fontana W. , Rules for modeling signal-transduction systems, Sci. STKE 2006 , re6 (2006).16849649
[2] Doi M. , Second quantization representation for classical many-particle system, J. Phys. A: Math. Gen. 9 , 1465 (1976).
[3] Doi M. , Stochastic theory of diffusion-controlled reaction, J. Phys. A: Math. Gen. 9 , 1479 (1976).
[4] Grassberger P. and Scheunert M. , Fock space methods for identical classical objects, Fortschr Phys. 28 , 547 (1980).
[5] Goldenfeld N. , Kinetics of a model for nucleation-controlled polymer crystal growth, J. Phys. A 17 , 2807 (1984).
[6] Peliti L. , Path integral approach to birth-death processes on, J Phys-Paris 46 , 1469 (1985).
[7] Lee B. and Cardy J. , Renormalization group study of the A+B→∅ diffusion-limited reaction, J. Stat. Phys. 80 , 10.1007/BF02179861 (1995).
[8] van Wijland F. , Field theory for reaction-diffusion processes with hard-core particles, Phys. Rev. E 63 , 10.1103/PhysRevE.63.022101 (2001).
[9] Mattis D. C. and Glasser M. L. , The uses of quantum field theory in diffusion-limited reactions, Rev. Mod. Phys. 70 , 979 (1998).
[10] Walczak A. , Mugler A. , and Wiggins C. , A stochastic spectral analysis of transcriptional regulatory cascades, Proc. Natl. Acad. Sci. USA 106 , 6529 (2006).
[11] Mjolsness E. , Time-ordered product expansions for computational stochastic system biology, Phys. Biol. 10 , 10.1088/1478-3975/10/3/035009 (2013).
[12] Weber M. and Frey E. , Master equations and the theory of stochastic path integrals, Rep. Prog. Phys. 80 , 10.1088/1361-6633/aa5ae2 (2017).
[13] Regev A. , Silverman W. , and Shapiro E. , Representation and simulation of biochemical processes using the π-calculus process algebra, Pac. Symp. Biocomput. , 459 (2001).11262964
[14] Blinov M. , Faeder J. , Goldstein B. , and Hlavacek W. , BioNetGen: software for rule-based modeling of signal transduction based on the interactions of molecular domains, Bioinformatics 20 , 3289 (2004).15217809
[15] Baeten J. , A brief history of process algebra, Theoretical Computer Science 335 , 131 (2005).
[16] Faeder J. R. , Blinov M. L. , and Hlavacek W. S. , Rule-based modeling of biochemical systems with BioNetGen, in Methods in Molecular Biology, Systems Biology, edited by Moly I. V. (Humana Press, 2009) Chap. 5 , pp. 113–167.
[17] Feret J. , Danos V. , Krivine J. , and Fontana W. , Internal coarse-graining of molecular systems, Proc Natl Acad Sci USA 106 , 6453 (2009).19346467
[18] Sneddon M. , Faeder J. , and Emonet T. , Efficient modeling, simulation and coarse-graining of biological complexity with NFsim, Nat Methods 8 , 177 (2011).21186362
[19] Chylek L. , Harris L. , Tung C. , Faeder J. , Lopez C. , and Hlavacek W. , Rule-based modeling: a computational approach for studying biomolecular site dynamics in cell signaling systems, Wiley Interdiscip. Rev. Syst. Biol. Med. 6 , 13 (2014).24123887
[20] Harris L. , Hogg J. , Tapia J. , Sekar J. , Gupta S. , Korsunsky I. , Arora A. , Barua D. , Sheehan R. , and Faeder J. , BioNetGen 2.2: advances in rule-based modeling, Bioinformatics 32 , 3366 (2016).27402907
[21] Boutillier P. , Maasha M. , Li X. , Medina-Abarca H. , Krivine J. , Feret J. , Cristescu I. , Forbes A. , and Fontana W. , The Kappa platform for rule-based modeling, Bioinformatics 34 , i583 (2018).29950016
[22] Andrei O. and Kirchner H. , A rewriting calculus for multigraphs with ports, Electronic Notes in Theoretical Computer Science 219 , 67 (2008).
[23] Gillespie D. T. , A general method for numerically simulating the stochastic time evolution of coupled chemical reactions, Journal of Computational Physics 22 , 403 (1976).
[24] Gillespie D. T. , Exact stochastic simulation of coupled chemical reactions, The Journal of Physical Chemistry 81 , 2340 (1977).
[25] Zurek W. H. , Decoherence, einselection, and the quantum origins of the classical, Reviews of Modern Physics 75 , 715 (2003), quant-ph/0105127.
[26] Schlosshauer M. , Quantum decoherence, Physics Reports 831 , 1 (2019), 1911.06282.
[27] Munsky B. and Khammash M. , The finite state projection algorithm for the solution of the chemical master equation, The Journal of Chemical Physics 124 , 044104 (2006).16460146
[28] Nicholson S. B. and Gingrich T. R. , Quantifying Rare Events in Stochastic Reaction-Diffusion Dynamics Using Tensor Networks, Physical Review X 13 , 041006 (2023), 2301.03717.
[29] Bintu L. , Buchler N. E. , Garcia H. G. , Gerland U. , Hwa T. , Kondev J. , Kuhlman T. , and Phillips R. , Transcriptional regulation by the numbers: applications., Curr Opin Genet Dev 15 , 125 (2005).15797195
[30] Bintu L. , Buchler N. E. , Garcia H. G. , Gerland U. , Hwa T. , Kondev J. , and Phillips R. , Transcriptional regulation by the numbers: models., Curr Opin Genet Dev 15 , 116 (2005).15797194
[31] Shea M. A. and Ackers G. K. , The OR control system of bacteriophage lambda. A physical-chemical model for gene regulation., J Mol Biol 181 , 211 (1985).3157005
[32] Dodd I. B. , Shearwin K. E. , Perkins A. J. , Burr T. , Hochschild A. , and Egan J. B. , Cooperativity in long-range gene regulation by the lambda CI repressor., Genes & Development 18 , 344 (2004).14871931
[33] Kuhlman T. , Zhang Z. , Saier M. H. , and Hwa T. , Combinatorial transcriptional control of the lactose operon of Escherichia coli., Proc Natl Acad Sci USA 104 , 6043 (2007), q-bio/0703056.17376875
[34] Cui L. , Murchland I. , Shearwin K. E. , and Dodd I. B. , Enhancer-like long-range transcriptional activation by lambda CI-mediated DNA looping., Proc Natl Acad Sci USA 110 , 2922 (2013).23382214
[35] Estrada J. , Wong F. , DePace A. , and Gunawardena J. , Information integration and energy expenditure in gene regulation., Cell 166 , 234 (2016).27368104
[36] Scholes C. , DePace A. H. , and Sanchez A. , Combinatorial Gene Regulation through Kinetic Control of the Transcription Cycle., Cell Syst 4 , 97 (2017).28041762
[37] Wong F. and Gunawardena J. , Gene Regulation in and out of Equilibrium., Annual review of biophysics 49 , 199 (2020).
[38] Kinney J. B. , Murugan A. , Callan C. G. , and Cox E. C. , Using deep sequencing to characterize the biophysical mechanism of a transcriptional regulatory sequence, Proceedings of the National Academy of Sciences 107 , 9158 (2010).
[39] Segal E. and Widom J. , From DNA sequence to transcriptional behaviour: a quantitative approach., Nature reviews Genetics 10 , 443 (2009).
[40] Belliveau N. M. , Barnes S. L. , Ireland W. T. , Jones D. L. , Sweredoski M. J. , Moradian A. , Hess S. , Kinney J. B. , and Phillips R. , Systematic approach for dissecting the molecular mechanisms of transcriptional regulation in bacteria, Proceedings of the National Academy of Sciences 115 , 201722055 (2018).
[41] Ireland W. T. , Beeler S. M. , Flores-Bautista E. , McCarty N. S. , Röschinger T. , Belliveau N. M. , Sweredoski M. J. , Moradian A. , Kinney J. B. , and Phillips R. , Deciphering the regulatory genome of Escherichia coli, one hundred promoters at a time, eLife 9 , e55308 (2020).32955440
[42] Ishigami Y. , Wong M. S. , Martí-Gómez C. , Ayaz A. , Kooshkbaghi M. , Hanson S. M. , McCandlish D. M. , Krainer A. R. , and Kinney J. B. , Specificity, synergy, and mechanisms of splice-modifying drugs, Nature Communications 15 , 1880 (2024).
