A universal method for analyzing copolymer growth

Qureshi, Benjamin; Juritz, Jordan; Poulton, Jenny M.; Beersing-Vasquez, Adrian; Ouldridge, Thomas E.

doi:10.1063/5.0133489

Polymers consisting of more than one type of monomer, known as copolymers, are vital to both living and synthetic systems. Copolymerization has been studied theoretically in a number of contexts, often by considering a Markov process in which monomers are added or removed from the growing tip of a long copolymer. To date, the analysis of the most general models of this class has necessitated simulation. We present a general method for analyzing such processes without resorting to simulation. Our method can be applied to models with an arbitrary network of sub-steps prior to addition or removal of a monomer, including non-equilibrium kinetic proofreading cycles. Moreover, the approach allows for a dependency of addition and removal reactions on the neighboring site in the copolymer and thermodynamically self-consistent models in which all steps are assumed to be microscopically reversible. Using our approach, thermodynamic quantities such as chemical work; kinetic quantities such as time taken to grow; and statistical quantities such as the distribution of monomer types in the growing copolymer can be directly derived either analytically or numerically from the model definition.

I. INTRODUCTION

Copolymers are polymers consisting of more than one type of monomeric unit; the order of these monomers in the chain defines the copolymer sequence. In general, copolymerization mechanisms can be classified into two main categories: free copolymerization that does not rely on a template,¹ as shown in Fig. 1(a); and templated copolymerization, in which a template (usually another copolymer) is used to bias the distribution of sequences produced, as shown in Figs. 1(b) and 1(c). Polymers produced via both types of mechanism are of relevance to both biological and industrial systems. In living systems, O-glycans are sequences of monosaccharides that grow by free copolymerization from serine or threonine amino acids.² They play a key role as a physical protective barrier for cells from pathogens, as well as participating in other cellular processes.^2,3 Free copolymerization is also a common method for producing plastics and rubbers in commercial and industrial systems.^4,5 In addition, there have been recent experimental designs for free copolymerization systems that produce specific products utilizing DNA-nanotechnology-based reaction schemes.^6,7

FIG. 1.

View large Download slide

Comparison of the different types of copolymerization mechanism with three types of monomer (blue, red, and yellow). (a) shows free copolymerization, (b) templated self-assembly, and (c) templated copolymerization with autonomous separation. In (b) and (c) the template is shown with squares and the growing polymer with circles. (d) Example of a more detailed reaction scheme used to select the next monomer. In each of these sub-figures, different colors represent different monomer types, with bonds colored accordingly when their strength might depend on the monomer. In (d), different activation states of the monomer undergoing incorporation are represented by different shapes. The dashed bubble indicates how the arbitrary set of reactions in (d) may replace the simple reaction surrounded by a dashed bubble in (c) for a more complex model.

Templated copolymerization is the mechanism by which DNA, RNA, and polypeptides are produced in DNA replication, RNA transcription, and protein translation, respectively. These processes are at the heart of the central dogma of molecular biology⁸ and are the basis of the informational and biochemical complexity of life. In DNA replication, DNA templates the production of copies of itself; in transcription, DNA templates the production of RNA; and, in translation, messenger RNA (mRNA) is the template for the production of a polypeptide.⁹ Inspired by these biological templated copolymerization mechanisms, there has been recent interest in designing synthetic systems that can produce other sequence-controlled molecules via templated copolymerization.^10–16

Free polymerization can be modeled as a Markovian growth process under which monomers bind to the end of a growing polymer at a certain rate. Early free copolymerization models^17–19 built on this framework to allow for copolymerization via the incorporation of multiple types of monomeric unit, as shown in Fig. 1(a), albeit with irreversible polymerization reactions. In particular, Mayo and Lewis¹⁹ emphasized that in polymerization models, if the monomer binding events are irreversible and their rates are conditional on the terminal monomer type, then intra-sequence correlations are generated within the copolymer.

Although the use of models with irreversible transitions is reasonable in many contexts, thermodynamically self-consistent models require all transitions to be microscopically reversible.²⁰ In particular, if a transition from state A to state B is possible, then transitions from B to A must also be possible. Models with fully microscopically reversible polymerization reactions, as in Fig. 1(a), are more challenging to analyze but can be interpreted in a thermodynamic sense.^1,21,22

Templates can affect the rate at which monomers are added or removed from a growing copolymer, and hence, templated copolymerization models can be more complex than free copolymerization models. When the template consists of just one type of templating monomer (homopolymer), a templated copolymerization process can be mapped onto a free copolymerization model. Furthermore, if one assumes some symmetries regarding interactions between monomers in the growing copolymer and those in the template (such as all complementary bonds have equal strength and all non-complementary bonds have equal strength), models of sequence-bearing templates may be mapped onto models with homopolymeric templates and, hence, to models of free copolymerization.^23–27

Templated copolymerization models can be further divided into two main categories: templated self-assembly Fig. 1(b)^,^1,27–41 and autonomously separating mechanisms Fig. 1(c).^23,24,42 Templated self-assembly models are those in which all the monomers in the growing copolymer remain bound to the template. In autonomously separating models, the growing copolymer detaches as it extends.^23–25 There has been recent interest in explicitly modeling autonomous separation in templated growth in an attempt to understand models that give a better description of transcription or translation.^23–25 In autonomously separating models, the simultaneous growth and separation of the copolymer and template mean that the copy-template interactions are not permanent, and therefore, free energy released from such interactions cannot be part of the driving force of polymerization. In addition, since these copy-template bonds are temporary, they cannot stabilize the accurate copy directly in the long time limit. Furthermore, an ensemble of accurate polymers is a lower entropy state than an ensemble of random polymers. These conditions mean that non-equilibrium driving is required to generate accurate copies of the template if the copies are to spontaneously detach.⁴² Moreover, the separation of the lagging tail from the template as the copolymer grows naturally causes intra-sequence correlations within the product.²³

The models described above are maximally coarse-grained, in that they treat the binding of monomers to the growing tip of the copolymer as a simple, usually single-step, process. However, more generally, one may wish to study models in which polymerization occurs via a more detailed series of steps, as in Fig. 1(d). For instance, in order to explain the high accuracy observed in biological polymer copying systems, Hopfield⁴³ and Ninio⁴⁴ independently introduced the concept of kinetic proofreading: a reaction motif in which a monomer undergoes a free energy consuming activation reaction before it is polymerized into the copolymer. The introduction of kinetic proofreading reaction motifs presaged the investigation of more complex copolymerization mechanisms.^35,45

In summary, models that allow for multiple monomer types, intra-sequence correlations, reversible reactions, and general, multi-step monomer inclusion reactions represent a wide class of copolymerization processes. Previous techniques^{17–19,22,23,27–35,39–41,46–51,66–69} have not allowed analysis of thermodynamically self-consistent models of generalized free copolymerization processes in which monomer addition is given by an arbitrarily complex network of reversible reactions with rates that may depend on the terminal monomer type, and templated copolymerization models with a high symmetry that can be mapped to these free processes. Investigating the most general type of model in this class would require simulation.

In this paper, we present a universal method for studying this large class of copolymerization models. Drawing on the work of Gaspard and Andrieux⁵² for analyzing linear copolymerization processes, and Hill^53,54 for analyzing absorbing Markov processes, we present analytical methods for extracting: explicit expressions for the probability of inclusion of a given monomer; the growth rate of a copolymerization process; and the chemical work done by the process. Our method removes the need to extract the same features by simulation and often produces simple, analytic results.

In Sec. II A, we review and refine methods relating to absorbing Markov chains that are crucial to understanding our approach. In Sec. II B 1, we present our method. In Sec. III, we apply the method to a few example processes to demonstrate its use and power when considering models with certain features. First, we apply the method to models for which the rate of adding new monomers only depends on the monomer type being added. Next, we apply the method to templated copolymerization systems with autonomous-separation that do not have non-equilibrium kinetic proofreading cycles. Finally, we solve a generalized version of Hopfield’s kinetic proofreading model applied to a templated copolymerization system with an autonomously separating product.

II. METHODS

A. Absorbing Markov chains

We begin by reviewing and adapting some diagrammatic techniques introduced by Hill to analyze absorbing Markov chains.^53,54 An absorbing Markov chain is a Markov chain for which any trajectory through its state space with arbitrary initial conditions will reach an absorbing state in finite time almost surely.⁵⁵ We can decompose the state space of an absorbing Markov chain into absorbing states $A$ and transient states $X$ such that the state space is $V = A \cup X$ ⁠. Let us denote the rate function that describes the chain as $K : V \times V \to R^{+}$ such that K(x, y) is the rate of the transition from state x to state y. Then, we denote a Markov process as the tuple (V, K).

Throughout this section, we shall refer to the absorbing Markov chain given in Fig. 2(a), which possesses two absorbing states and non-trivial cycles, for illustrative purposes.

FIG. 2.

View large Download slide

Graphical representations of an absorbing Markov process to illustrate the methodology outlined in Sec. II A. (a) Example absorbing Markov process $(X \cup A, K)$ ⁠, with two absorbing states, $A = {A, B}$ ⁠, and four transient states $X = {1, 2, 3, 4}$ ⁠. (b) The closed process starting at state 1, $(X, K_{1})$ ⁠. (c) The cycle process $({c} \cup X / {1, 2, 3} = {c, 4}, K_{1, C})$ for the cycle C = 1 → 2 → 3 → 1 or C′ = 1 → 3 → 2 → 1.

1. Expectations of an absorbing process are steady state averages of a “closed process”

We will derive expressions for four main quantities: the probabilities of reaching certain absorbing states, the expected time taken to absorption, the expected net number of times traversing a given edge before absorption, and the expected number of times that a trajectory goes round a cycle before absorption. These quantities depend on the starting (transient) state $σ \in X$ and can be found in terms of the “closed” process.⁵⁴ The closed process is a modified version of an absorbing Markov chain in which transitions to the absorbing states are redirected to the starting state. Figure 2(b) shows the closed process starting at state 1 of our example absorbing chain of Fig. 2(a). The closed process for a Markov process $(X \cup A, K)$ starting at state σ is a new Markov process $(X, K_{σ})$ with a rate function given by

K_{σ} (x, σ) = K (x, σ) + \sum_{A \in A} K (x, A)

(1)

for $x \in X$ and agreeing with K on $X \times X / {σ}$ ⁠.

The closed process has a unique stationary distribution for the following reasons. From the definition of an absorbing Markov chain, there exists a path from any state to an absorbing state, taking finite time. Thus, in the closed process, there is a path from any state to the starting state, taking finite time. The set of states including the starting state and all those that may be reached from the starting state is, therefore, positive recurrent and furthermore, this set is the only recurrent set of states and will be reached from any other state. Since there is only one recurrent set of states, there is a unique stationary distribution.⁵⁵

Expected quantities of an absorbing Markov chain, such as the expected probability that a particular absorbing state is reached, can be found in terms of steady state quantities in the closed process. Whenever a trajectory of the original process reaches an absorbing state, in the closed process that same trajectory would have been reset back to the starting state. Hence, running the closed process for long times is equivalent to generating many independent trajectories to absorption for the original chain. Thus, averaging quantities in the steady state of the closed process is equivalent to taking expectations over independent trials of quantities in the absorbing chain. It is worth noting that the dependence of expected quantities on the starting state is encoded in the definition of the closed process. Finally, we can see that the definition of the closed process may permit self-transitions σ → σ, which, for continuous time Markov processes, have little meaning. However, for the purposes of calculating steady state probabilities of the closed process they may be ignored.

2. Steady state averages of the closed process are calculated using the Markov chain tree theorem

Given that we can turn the calculation of expectations of absorbing processes into steady state averages over closed processes, we can make use of tools developed for analyzing the steady state of Markov processes, such as the Markov chain tree theorem (MCTT).⁵⁶ The MCTT states that the steady state distribution of a Markov chain with a unique stationary distribution may be found by summing over rooted spanning trees of the process, where the transition rates are taken as weights on the edges of the graph. Explicitly, let $G$ be a directed weighted graph, with weight K(e) for an edge e of $G$ ⁠. A spanning tree of $G$ ⁠, rooted at a vertex ν, is a subgraph of $G$ with no cycles that connects all the vertices of $G$ and for which the out degree of every vertex, except ν, is one. The sets of spanning trees rooted at nodes 3 and 4 of the closed process given in Fig. 2(b), are shown in Fig. 3. Denote the set of all spanning trees rooted at x by $T (x)$ ⁠. The MCTT states that the steady state probability π(x) to be in state x is given by

π (x) = \frac{\sum_{T \in T (x)} \prod_{e \in T} K (e)}{\sum_{v \in X} \sum_{T \in T (v)} \prod_{e \in T} K (e)},

(2)

with e ∈ T representing the edges of the tree. The denominator here is simply a normalization constant.

FIG. 3.

View large Download slide

Spanning trees of the closed processes rooted at (a) node 3 and (b) node 4 derived from Fig. 2(b), with nodes labeled in the first spanning tree and all other trees following the same positioning. The spanning trees have been arranged in terms of the self-avoiding walk between nodes 1 and 3 for the trees rooted at node 3 and arranged in terms of the self-avoiding walks between nodes 1 and 4 for the trees rooted at node 4. More details on the relationship between self-avoiding walks and spanning trees are given in Appendix A.

We can define steady state currents from the steady state distribution of the closed process that corresponds to expected currents of the absorbing chain. Let a subscript σ denote quantities in the closed process starting at state σ. Then, π_σ is the steady state probability distribution and K_σ the rate function. The current along a given edge, e = x → y, is given by the probability to be in state x, π_σ(x), multiplied by the rate along said edge. We can, therefore, write the steady state current along all edges that originally led to absorbing states as

J_{Tot} (σ) = \sum_{A \in A} \sum_{x \in X} π_{σ} (x) K (x, A),

(3)

where, as before, $X$ is the set of transient states and $A$ ⁠, the set of absorbing states. J_Tot(σ) is the expected total current to absorbing states from state σ, and, therefore, its reciprocal is the expected time to absorption.

For the example process shown in Fig. 2, we present the spanning trees of the corresponding closed process rooted at node 3 and 4 in Fig. 3. Given the spanning trees, we can directly write down the total current to the absorbing states as

\begin{matrix} J_{Tot} (1) & = & \frac{1}{N} [k_{A} [r_{12} r_{24} r_{43} + r_{12} r_{23} (r_{42} + r_{43} + k_{B}))) \\ + ((r_{13} ((r_{43} + k_{B}) (r_{21} + r_{23} + r_{24}) + r_{42} (r_{21} + r_{23}))]) \\ + (k_{B} [r_{12} r_{23} r_{34} + r_{12} r_{24} (r_{31} + k_{A} + r_{32} + r_{34}))) \\ + ((r_{13} r_{34} (r_{21} + r_{23} + r_{24}) + r_{13} r_{32} r_{24}]], \end{matrix}

(4)

where $N$ is the normalization term, given in Appendix B. The terms multiplied by k_A are the partial current to absorbing state A, i.e., the current along transition 3 → A, coming from the trees rooted at node 3, and equivalently for k_B with state B, i.e., the current along transition 4 → B, coming from trees rooted at node 4.

3. Absorbing probabilities

Given a Markov chain with multiple absorbing states, we can ask for the probability of absorption in each absorbing state in the long time limit. The probability that a trajectory eventually ends in a specific absorbing state can be calculated from the closed process, by dividing the expected current along transitions that originally led to the absorbing state in question by J_Tot(σ) [Eq. (3)]. Therefore, the absorption probabilities can be written

P [σ \to A] = \frac{\sum_{x \in X} π_{σ} (x) K (x, A)}{\sum_{B \in A} \sum_{x \in X} π_{σ} (x) K (x, B)},

(5)

using the notation $P [σ \to A]$ to denote probability of being absorbed to A given that the trajectory started in state σ. It is worth noting here that given that this quantity is a ratio of currents, there is a factor of π_σ in both the denominator and the numerator of the expression. In practice, we see that the normalization factor from the MCTT [Eq. (2)] cancels out, which simplifies the quantities in the calculation.

For our example process shown in Fig. 2, we can use the partial currents to absorbing states A and B to write down the absorbing probabilities

\begin{matrix} P [1 \to A] & = & \frac{1}{N J_{Tot} (1)} k_{A} [r_{12} r_{24} r_{43} + r_{12} r_{23} (r_{42} + r_{43} + k_{B})) \\ + (r_{13} ((r_{43} + k_{B}) (r_{21} + r_{23} + r_{24}) + r_{42} (r_{21} + r_{23}))], \\ P [1 \to B] & = & \frac{1}{N J_{Tot} (1)} k_{B} [r_{13} r_{34} (r_{21} + r_{23} + r_{24}) + r_{13} r_{32} r_{24}) \\ + ((r_{12} r_{23} r_{34} + r_{12} r_{24} (r_{31} + k_{A} + r_{32} + r_{34})]] . \end{matrix}

(6)

The normalization factor, $N$ ⁠, propagated through from Eq. (2), conveniently cancels out with the $1 / N$ implicit in J_tot.

4. Counting edge and cycle transitions

In this subsection, we shall calculate the expected net number of times traversing a given edge of an absorbing Markov process before absorption. In addition, we shall calculate the expected number of times a non-recurrent cycle of an absorbing Markov process is traversed before absorption. Both of these will be of use later in defining a notion of chemical work.

To calculate the net number of times crossing a given edge, we find the expected current along the transition x ⇋ y, between states x and y of an absorbing process, $(X \cup A, K)$ ⁠, as in Sec. II A. The expected current through this edge, denoted J_x⇋y(σ), given starting in state $σ \in X$ ⁠, can be calculated from the closed process, $(X, K_{σ})$ ⁠, as in Eq. (1), as the difference between the steady state probability to be in state x multiplied by the rate from x → y and the steady state probability to be in state y multiplied by the rate from y → x,

J_{x ⇋ y} (σ) = π_{σ} (x) K (x, y) - π_{σ} (y) K (y, x) .

(7)

The net number of times traversing this edge (number of observed transitions x → y - number of observed transitions y → x) before absorption is, then, just the ratio between this current and total current to absorbing states,

N_{x ⇋ y} (σ) = \frac{J_{x ⇋ y} (σ)}{J_{Tot} (σ)} .

(8)

The current, given in Eq. (7), is intimately linked to the notion of cycles as pointed out by Wachtel et al.⁵⁷ and detailed in Appendix C. Thus, we also wish to find the expected number of times traversing a non-recurrent cycle. We define a non-recurrent cycle for a Markov chain to be a cycle of states, where each state, aside from the originating state, does not appear more than once in the cycle. For example, the cycle A → B → C → D → A is non-recurrent, but A → B → C → D → B → A is recurrent. Note that the originating state is arbitrary, and so A → B → C → D → A is equivalent to B → C → D → A → B. For a stationary process, the expected frequency with which a cycle is completed can be calculated from the one-way cycle current,^54,58 which is the probability current going around the cycle. For a chosen non-recurrent cycle, the one-way cycle current can be calculated diagrammatically from three terms. First, a cycle term given by the product of rates around the cycle in the chosen direction. Second, a spanning tree term that can be found by collapsing the nodes in the cycle into a single node [in Fig. 2(c), the cycle 1 → 2 → 3 → 1 has been collapsed in this way] and finding the sum of spanning trees of this new graph rooted at the collapsed cycle node. Finally, there is a normalization factor, which is the same normalization factor as for the current, $N$ ⁠. Explicitly, consider an absorbing Markov chain $G = (X \cup A, K)$ and its closed process starting at $σ \in X$ ⁠, $G_{σ} = (X, K_{σ})$ ⁠. Let C denote both the set of edges and set of nodes of a cycle in the closed process. To calculate the spanning tree term for the one-way cycle current, construct a new Markov chain, the cycle process $G_{σ, C} = ({c} \cup (X / C), K_{σ, C})$ ⁠, where ${c} \cup (X / C)$ is the set of transient states of the original process with the states in the cycle replaced by the single node c and K_σ,C is given by

\begin{aligned} K_{σ, C} (x, c) = \sum_{i \in C} K_{σ} (x, i) \\ K_{σ, C} (c, x) = \sum_{i \in C} K_{σ} (i, x) \\ K_{σ, C} (c, c) = 0 \end{aligned}

(9)

for $x \in X / C$ and agreeing with K_σ elsewhere. The cycle process for the cycle C = 1231 (or C′ = 1321) of the example system in Fig. 2(a) is shown in Fig. 2(c). Let $T_{σ} (x), T_{C} (x)$ be the sets of spanning trees rooted at x of the closed process, $G_{σ}$ ⁠, and cycle process, $G_{σ, C}$ ⁠, respectively. Then, the cycle current is given by⁵⁸

J_{Cyc} (σ, C) = \frac{\overset{Cycle}{\overset{︷}{(\prod_{e \in C} K (e))}} \overset{Spanning Trees}{\overset{︷}{\sum_{T \in T_{C}} \prod_{e \in T} K_{σ, C} (e)}}}{\underset{normalization}{\underset{︸}{\sum_{x \in X} \sum_{T \in T_{σ} (x)} \prod_{e \in T} K_{σ} (e)}}} .

(10)

Note that, for the cycle term, the edges are taken from the original process rather than the closed process. Given the cycle current for the closed process, the expected number of circulations of the cycle before absorption, N_Cyc(σ, C), is the ratio of the cycle current to the total current to absorbing states

N_{Cyc} (σ, C) = \frac{J_{Cyc} (σ, C)}{J_{Tot (σ)}} .

(11)

For our example process, the expected number of circulations of C = 1231 is

N_{Cyc} (1, C) = \frac{(r_{12} r_{23} r_{31}) (r_{42} + r_{43} + k_{B})}{N J_{Tot} (1)},

(12)

with the same implicit cancelation of normalization as before, since $J_{Tot} \propto \frac{1}{N}$ ⁠.

For an absorbing process starting at a given state, we may divide the cycles into internal and external cycles. External cycles are those which appear in the closed process and involve edges which were absorbing edges in the original process. The set of all cycles, sorted into internal and external, for the example process Fig. 2, is shown in Appendix D. The external cycles correspond to the pathways from the starting state to an absorbing state. Therefore, the expected number of times traversing an external cycle before absorption will be at most one and corresponds to the probability of following a given path to absorption. Furthermore, the sum of Eq. (11) over all external cycles will be one.

B. Copolymer methods

1. Philosophy of coarse-graining complex underlying copolymerization reaction networks

Armed with the techniques for solving absorbing Markov chains, here, we set out the method for the analysis of copolymerization processes. Gaspard and Andrieux⁵² presented a method to analyze Markov polymerization processes in which each monomer is added in a single step [i.e., if the internal reaction network shown in Fig. 1(d) were trivial], assuming long polymers. We shall present a method for mapping more complex models for the individual polymerization step onto coarse-grained descriptions that can be analyzed using this framework and, subsequently, show how to get the behavior of the full model from the results.

Consider a growing copolymer with M monomer types, which are assumed to be present in the environment at fixed concentrations. At a coarse grained level, we can define a state space of finite length sequences {x₁x₂…x_l|x_i ∈ {1, 2, …, M}, l ≥ 0}, where l is the length of the sequence. Let us refer to the coarse-grained states in this state space as completed states. On this coarse-grained level, a sequence of length l may increase in length by one unit by polymerizing one of M units at the growing tip (x₁x₂…x_l → x₁x₂…x_lx_l+1) or it may decrease in length by one unit (x₁x₂…x_l → x₁x₂…x_l−1). Such a coarse grained model is depicted in Figs. 1(a)–1(c) for free polymerization, templated self-assembly, and templated polymerization with simultaneous separation.

In general, copolymerization processes may be best described by models in which the underlying copolymerization reaction networks are complex, featuring multiple sub-steps in arbitrarily complex networks connecting the completed states, as suggested in Fig. 1(d). Hence, overall, we could consider a copolymerization process as having a tree-like structure with networks of reactions connecting completed states, as in Fig. 4. Such a class of models is wide-reaching, with many examples from the literature included in this class.^{17–19,22,23,27,28,30–35,39–41,45–51}

FIG. 4.

View large Download slide

(a) The step-wise process for an arbitrary model with three monomer types. This step-wise process is for a copolymer &xy. The edges colored red are the completion edges. The flower like structure of the step-wise process can be seen with four petals, each connected at the starting state &xy. (b) One of the petals of the step-wise process, which we use to define Λ^±(z, y). Λ⁺(z, y) is defined as the sum of spanning trees rooted at the rightmost state &xyz and Λ⁻(z, y) is defined as the sum of the trees rooted at the leftmost state &xy. (c) One of the petals [connecting &xy to &xyz, as in (b)] which has been linked back to the starting state. This graph is used to define Q(z, y) as the sum of spanning trees rooted at the leftmost state &xy.

We will define a Markov process at the level of the coarse-grained completed states that, by construction, preserves the probabilities of the transitions between the completed states of the fine-grained process and, therefore, preserves the statistics of the sequences produced. The coarse-grained Markov process does not preserve the distribution of the transition times between the completed states implied by the fine-grained model, which will, in general, be non-Markovian. Moreover, it does not provide fine-grained information on trajectories between the coarse-grained completed states. However, the temporal details and information about the fine-grained dynamics can be added back in at a later stage, once the statistics have been analyzed at the coarse-grained level.

2. Identifying propensities in the coarse-grained model

We find the transition rates of the coarse-grained model (hereafter labeled as propensities to avoid confusion with the underlying rates of the fine-grained process) by considering first passage problems between completed states. From a given completed state, there are M + 1 completed states that may be reached, corresponding to the M possible additions of a monomer and the removal of the monomer currently at the tip of the copolymer. For a first passage problem, we can convert each of these reachable completed states into an absorbing state by removing the transitions out of said states, as in Fig. 4(a), in the same vein as Cady and Qian.⁵⁹ Let us refer to this absorbing Markov process as the step-wise process and define step to mean the addition/removal of a monomer.

We shall work with the assumption that the transition rates depend on the two monomers at the growing tip of the copolymer following.^{1,17–19,21–23,36–38,52,60} There will therefore be M² flavors of this process corresponding to the combinations of the two terminal monomers of the copolymer, the central state &xy (here & represents an arbitrary sequence). We wish to find the absorbing probabilities, $P [& x y \to & x y z], z \in {1, \dots M}, P [& x y \to & x]$ given an initial condition of the central state &xy. As outlined in Sec. II B 1, Eq. (5), we can find these probabilities by constructing the closed process and finding sums over spanning trees rooted at different states. The step-wise process has M + 1 petal-like graphs each connected to the central state but disconnected from each other. Due to this structure, any sums over spanning trees of the full process will factorize into a product of sums over spanning trees of the petals. Thus, we find that the absorbing probabilities take the following form:

\begin{aligned} P [& x y \to & x y z] = \frac{1}{N} Λ^{+} (z, y) [\prod_{z^{'} \neq z} Q (z^{'}, y)] Q (y, x), \\ P [& x y \to & x] = \frac{1}{N} Λ^{-} (y, x) \prod_{z} Q (z, y) . \end{aligned}

(13)

Here, z ∈ {1, …, M}, $N$ is the normalization factor from Eq. (2); Λ⁺(z, y) is the sum over spanning trees of the petal connecting states monomers &xy and &xyz, rooted at the forward completed state &xyz; Λ⁻(y, x) is the sum over spanning trees of the petal connecting states monomers &x and &xy, rooted at the backward completed state &x, Fig. 4(b); and Q(y, x) is the sum over spanning trees of the petal connecting states &x and &xy, linked back to the central state and rooted at the central state, i.e., with edges redirected to the starting state as in the closed process, as in Fig. 4(c). Since Q is a sum over spanning trees rooted at the node to which edges have been redirected, the sum takes the same form for both the forward and backward petals, only depending on which two completed states it is connecting.

From these probabilities, we see that choosing propensities ω_±yx for the transitions $& x \overset{ω_{+ y x}}{\to} & x y$ and $& x y \overset{ω_{- y x}}{\to} & x$ such that

ω_{\pm y x} = \frac{Λ^{\pm} (y, x)}{Q (y, x)}

(14)

not only preserves the ratios of probabilities of transitions to completed states, but also ensures that ω_±yx only depends on monomers x and y.

We note here that this coarse graining process is different from lumping,^55,61 in which the state space is reduced while attempting to retain trajectory dynamics. In our approach, the coarse-grained process does not reproduce the dynamics of the fine-grained process, only the statistics of the completed states that are visited. However, dynamic quantities may be extracted exactly from the step-wise process, as we show in Sec. II B 4.

3. Solving the coarse-grained model

We now use the methods developed by Gaspard and Andrieux⁵² to solve the coarse-grained Markov model over the completed states with propensities ω_±yx. Gaspard and Andrieux’s approach considers a frame of reference that is comoving with the tip of the growing polymer and assumes that the state of the tip and nearby monomers reaches a stationary distribution to derive quantities at this steady state, such as the set of tip incorporation velocities v_x (the rates of adding monomers to a copolymer &x), the tip probabilities μ(x) (the probability at a given time that the growing polymer is in state &x), and the pair tip probabilities μ(x, y) (the probability of being in state &xy). However, we note that the time-dependent information is not physical at this stage due to the coarse-graining process. The above quantities are found from solving the following equations:⁵²

v_{x} = \sum_{y = 1}^{M} \frac{ω_{+ y x} v_{y}}{ω_{- y x} + v_{y}},

(15)

μ (x) = \sum_{y = 1}^{M} \frac{ω_{+ x y}}{ω_{- x y} + v_{x}} μ (y),

(16)

μ (x, y) = \frac{ω_{+ y x}}{ω_{- y x} + v_{y}} μ (x) .

(17)

Using μ and v, we can calculate the statistics of the copolymer sequence far behind the growing tip.⁵² We note that the distribution of monomers at the tip μ(x) is different from the distribution of monomers at sites behind the tip; we assume that this distribution reaches some limit far behind the growing tip in the bulk of the copolymer. This limiting distribution describes the probability that a monomer in the bulk of the copolymer takes a value x. Using ɛ(x) to denote the frequency of monomer x in the bulk of the copolymer,⁵²

ε (x) = \frac{μ (x) v_{x}}{\sum_{y} μ (y) v_{y}} .

(18)

We may similarly define ɛ(y|x) as the probability that in the bulk of the copolymer, a monomer y is observed given a monomer x behind it. ɛ(x) and ɛ(y|x) fully characterize the statistics of the bulk copolymer since under our assumptions, transitions only depend on the two monomers at the tip; the completed copolymer sequence is itself a Markov chain.²³

4. Extracting properties of the fine-grained model from the solution of the coarse-grained model

The easiest quantities to extract are the frequencies of monomers in the bulk of the copolymer. These quantities are identical in the coarse-grained and fine-grained models, since the coarse-graining preserves the statistical distribution of the sequences produced. Therefore ɛ(x), as defined in Eq. (18), and ɛ(y|x) apply directly to the fine-grained process.

The tip probabilities μ above give the fraction of the time spent in each tip state in the coarse-grained model. However, the coarse-grained model will not reproduce the time series of the fine-grained model, but only the sequences of completed states visited. We, therefore, quotient out the lifetime τ(x, y) of the tip state (x, y) to obtain the frequency with which the tip states are visited in the coarse-grained model,

ξ (x, y) = \frac{1}{\sum_{x^{'}, y^{'} = 1}^{M} \frac{μ (x^{'}, y^{'})}{τ (x^{'}, y^{'})}} \frac{μ (x, y)}{τ (x, y)},

(19)

τ (x, y) = \frac{1}{ω_{- y x} + \sum_{z = 1}^{M} ω_{+ z y}} .

(20)

This frequency defines a new tip distribution ξ. ξ(x, y) is the frequency that a given pair of monomers x, y is observed at the tip of the growing copolymer in the sequence of transitions. This distribution, ξ(x, y), applies to both the coarse-grained model and the sequence of completed states visited in the full fine-grained model. It can, therefore, be used to find averages of the key dynamic properties.

For example, we can calculate the probability P that a growing copolymer increases in length at each step of the step-wise process. P is calculated by averaging the probability of adding a monomer over the possible states &xy:

P = \sum_{x, y = 1}^{M} ξ (x, y) \frac{\sum_{z = 1}^{M} ω_{+ z y}}{ω_{- y x} + \sum_{z = 1}^{M} ω_{+ z y}} .

(21)

Upon averaging out the sequence information, we may treat the growth of a polymer as a random walk with probability P of stepping forwards and (1 − P) of stepping back. We can find the expected number of monomer inclusion/removal steps per net forward step as 1/(2P − 1) (for proof see Appendix E). A number of quantities scale with the total number of steps rather than the net number of steps, making the number of steps per net forward step a necessary quantity. For example, in order to find the expected time taken per net forward step, one can find the expected time to absorption for the step-wise process T(x, y) for a copolymer in state &xy by calculating 1/J_Tot(&xy) for the step-wise process using Eq. (3). Then, the expected time per net forward step is

τ_{step} = \frac{1}{2 P - 1} \sum_{x, y = 1}^{M} ξ (x, y) T (x, y) .

(22)

1/τ_step is, therefore, the physical average growth rate of the copolymer in the fine-grained model.

We may also calculate the chemical work done by the system in producing the copolymer. In a purely chemical system, with no time-varying externally applied protocols, the entropy increase of the universe is given by the decrease in the generalized free energy of the chemical system, including any coupled reservoirs of fuel molecules.²⁰ Since the total free energy must decrease, any increase in one contribution must be paid for by a decrease of at least the same magnitude in another contribution. It is common to describe the latter subsystem as doing work on the former.

For the polymerization systems analyzed here, the generalized free energy can be split into a term corresponding to the chemical free energy of the system averaged over the uncertain state of the system and a term related to the entropy arising due to the uncertainty of the state occupied.⁶²

G = \sum_{a} p (a) G_{chem} (a) + \sum_{a} p (a) \ln p (a),

(23)

where we use natural units such that k_BT = 1. Here, a is a chemical state of the system as a whole, G_chem(a) is the chemical free energy of state a, and p(a) is the probability that the system occupies the state a. G_chem(a) = −ln Z_a, where Z_a is the partition function of the system (explicitly including any large chemical buffers) restricted to the chemical state a. Z_a represents the contribution of concentrations and bond strength to the favorability of a molecular state. The principle of detailed balance²⁰ states that the chemical free energy change associated with a transition from a to b is given by

G_{chem} (b) - G_{chem} (a) = - \ln (\frac{K (a, b)}{K (b, a)}) .

(24)

The second term in Eq. (23) is information theoretic in character; it is equal to the negative of the Shannon entropy associated with the distribution over chemical states. For the systems studied here, in which we consider infinitely long copolymers that have reached steady state growth, the only relevant contribution to this term is the increase in Shannon entropy of the copolymer sequence produced as the polymer gets longer. Since the copolymer sequence is itself a discrete time Markov chain²³ the additional entropy per net forward step (the entropy rate) can be readily calculated,⁶³

H = - \sum_{x, y = 1}^{M} ε (x) ε (y | x) \ln ε (y | x),

(25)

with x, y representing the monomer types. Since the purpose of a copolymerization system is often to produce a low entropy (or “accurate”) sequence, it is reasonable to think of the chemical free-energy decrease per net forward step as the chemical work done to reduce the information entropy of Eq. (25) below that of a uniform, random polymer. Extending the definition provided by Poulton et al.,²³ we may define the efficiency of copolymerization as

η = \frac{\ln M - H}{\ln M + W_{chem}} \leq 1,

(26)

where ln M is the entropy per monomer (or entropy rate) of a uniform, random copolymer with M monomer types and $W_{chem}$ is the average decrease in chemical free energy per net forward step. This efficiency is, then, the ratio between the entropy drop due to the accuracy of the copolymer compared to that of a random one (ln M − H) and the chemical work used to drive the system $(W_{chem})$ above that required to make a random copolymer in equilibrium is (−ln M).⁶⁴

The expected work done during a transition adding or removing a monomer given starting in completed state &xy can be calculated by summing the contribution from Eq. (24) multiplied by the expected net current along the edge a ⇋ b prior to absorption over all the edges in the step-wise process,

\begin{align} w_{chem} (x, y) = \\ - Δ G_{chem} (x, y) = \sum_{b > a} \ln (\frac{K (a, b)}{K (b, a)}) N_{a ⇋ b} (& x y), \end{align}

(27)

where N_a⇋b(&xy) is the expected net number of times traversing edge a ⇋ b before absorption given started in the central state of the step-wise process &xy, as in Eq. (8). This sum will also require contributions from edges that lead to absorbing states. For such edges, the rate for the reverse transition in the logarithm of Eq. (27) is the rate from the full process.

Equivalently, however, as outlined in Appendix C, we may find this chemical work by considering the non-recurrent cycles of the process.⁵⁷ For a given internal cycle, C, we may define the affinity²⁰

A (C) = \ln \frac{\prod_{e \in C} K (e)}{\prod_{e \in C^{'}} K (e)},

(28)

where the products are over the edges e composing the cycle and C′, which is the cycle with the edges in a reversed direction. For external cycles, we may define the affinity in the same way, inferring the rate for the reversed edge of the transition to absorbing states from the full process. The expected work done before absorption of the cycle C having started in the state &xy is

\begin{align} w_{chem} (x, y) = \\ - Δ G_{chem} (x, y) = \sum_{C} A (C) \frac{J_{Cyc} (σ, C) - J_{Cyc} (σ, C^{'})}{J_{Tot (& x y)}} . \end{align}

(29)

Averaging w_chem(x, y) with ξ and multiplying by the expected number of steps per net forward step gives the expected chemical work done per net forward step,

W_{chem} = \frac{1}{2 P - 1} \sum_{x, y = 1}^{M} ξ (x, y) w_{chem} (x, y) .

(30)

Furthermore, the forms of Eqs. (22) and (30) may be applied to an arbitrary quantity for which one can find the expected value in the step-wise process starting in state &xy. Let this arbitrary quantity be A(x, y). One can then average this quantity using the distribution ξ to obtain the expected value of the quantity per step. Then, if appropriate, multiplying by 1/(2P − 1) gives the expected value of the quantity per net forward step. In practice, as shall be seen in Sec. III C, since the quantities we wish to calculate may be written in terms of sums over spanning trees, the quantities for the step-wise process may be written as a sum over the terms per petal, with the quantity for a given petal factorizing into some quantity that depends on the petal multiplied by Q’s for the other petals.

5. Stalled growth

Explicit simulation of copolymer growth is particularly challenging in regimes where P ≳ 0.5, since many backward and forward steps are taken per net forward step. At P = 0.5, then, the process will not reliably produce copolymers; for P < 0.5 polymers will tend to shrink. In general, for P = 0.5, we can say the model has stalled. Our approach is particularly beneficial in this case; indeed, it is possible to check whether a model is at the stall point by considering an M × M dimensional matrix of the ratios of forward to backward propensities,⁵² $Z_{y x} = (\frac{ω_{+ y x}}{ω_{- y x}}) = \frac{Λ^{+} (y, x)}{Λ^{-} (y, x)}$ ⁠. The model is at the stall point if and only if

\det (Z - 1_{M}) = 0,

(31)

where $1_{M}$ is the M × M identity matrix. The polymer will shrink if $\det (Z - 1_{M}) < 0,$ Since Z gives the ratios of adding a monomer to removing one, this condition essentially says that models will stall if the total rate of adding a monomer is equal to the total rate of removing one.

In a typical model, there exists at least one parameter that controls the driving. Often this parameter is related to the backbone strength of the polymer produced: e.g., the free energy drop associated with the formation of a generic backbone bond ΔG_pol. This parameter will be present in the rates of each external cycle so that by tuning it, the model can be moved all the way from stalling to irreversible growth, whereby monomers cannot be removed once polymerized. If such a parameter exists, we may rephrase the stall condition, given in Eq. (31), in terms of this parameter. For example, for the case of the parameter being ΔG_pol, we may find some threshold value Γ such that the model will stall for ΔG_pol = Γ.

6. Limiting behavior

We shall note two limits for which we may give analytic expressions for the frequency of monomer types in the copolymer bulk for all models. First, consider the case that the system is at the stall point [Eq. (31)]. In general, entropy production can still occur within cycles in the step-wise process; therefore, these frequencies cannot be determined from equilibrium arguments and are non-trivial. Nonetheless, at the stall point, we may express the monomer frequencies in the bulk relatively simply. The frequency ɛ_stall(x) of monomer x is proportional (up to normalization) to the cofactor of the diagonal element (corresponding to monomer x) of the matrix $(1_{M} - Z)$ ⁠, as proven in Appendix F. For example, for M = 2,

\begin{aligned} ε_{stall} (1) & \propto 1 - \frac{ω_{+ 22}}{ω_{- 22}}, \\ ε_{stall} (2) & \propto 1 - \frac{ω_{+ 11}}{ω_{- 11}}, \end{aligned}

(32)

and for M = 3, we have

\begin{aligned} ε_{stall} (1) & \propto (1 - \frac{ω_{+ 22}}{ω_{- 22}}) (1 - \frac{ω_{+ 33}}{ω_{- 33}}) - \frac{ω_{+ 23}}{ω_{- 23}} \frac{ω_{+ 32}}{ω_{- 32}}, \\ ε_{stall} (2) & \propto (1 - \frac{ω_{+ 11}}{ω_{- 11}}) (1 - \frac{ω_{+ 33}}{ω_{- 33}}) - \frac{ω_{+ 13}}{ω_{- 13}} \frac{ω_{+ 31}}{ω_{- 31}}, \\ ε_{stall} (3) & \propto (1 - \frac{ω_{+ 11}}{ω_{- 11}}) (1 - \frac{ω_{+ 22}}{ω_{- 22}}) - \frac{ω_{+ 12}}{ω_{- 12}} \frac{ω_{+ 21}}{ω_{- 21}} . \end{aligned}

(33)

On the other end of the spectrum, we can also solve for the monomer bulk frequencies in the irreversible limit, where ω_−yx = 0 for all x, y. Intuitively, we could consider the Markov process on the state space {1, …, M} representing copolymers with a given monomer at its tip and the transitions between those states with rates K_irrev(x → y) = ω_+yx. The steady state of this process will give the time dependent frequencies of having a given monomer at the tip of the copolymer. Therefore, dividing by the time spent in each state will give the bulk frequencies. A nice way to write out these frequencies in the style of the methods described thus far is as a sum over the spanning trees on the complete graph on M vertices with rate functions K_irrev(x, y) = ω_+yx. Explicitly, we may write these frequencies (up to normalization) as

ε_{irrev} (x) \propto (\sum_{T \in T (x)} \prod_{e \in T} K_{irrev} (e)) \sum_{y = 1}^{M} ω_{+ y x},

(34)

where $T (x)$ is the set of spanning trees of the complete graph on M vertices. This expression is derived formally in Appendix G. For example, with M = 2,

\begin{aligned} ε_{irrev} (1) & = \frac{ω_{+ 12} (ω_{+ 11} + ω_{+ 21})}{ω_{+ 12} (ω_{+ 11} + ω_{+ 21}) + ω_{+ 21} (ω_{+ 12} + ω_{+ 22})}, \\ ε_{irrev} (2) & = \frac{ω_{+ 21} (ω_{+ 12} + ω_{+ 22})}{ω_{+ 12} (ω_{+ 11} + ω_{+ 21}) + ω_{+ 21} (ω_{+ 12} + ω_{+ 22})} . \end{aligned}

(35)

7. Simplification for factorizable propensities

The presented method applies to arbitrary complex copolymerization models obeying the structure of Fig. 4. However, if we make some further common assumptions, much of the analysis simplifies. For example, consider the case in which the ratios of propensities may be factored,

\frac{ω_{+ y x}}{ω_{- y x}} = \frac{Λ^{+} (y, x)}{Λ^{-} (y, x)} = Y (y) X (x),

(36)

where Y is a function of monomer y only and X is a function of monomer x only. Intuitively, such a condition holds in the cases where there are no direct, type-dependent interactions between monomers in the growing polymer, such as when monomers only interact with a template.^{1,23,24,27–38} Under such an assumption, multiple calculations simplify, see Appendix H. For example, the stall condition becomes simply that the model will stall at

\sum_{x} X (x) Y (x) = 1 .

(37)

Bulk frequencies at stall are just,

ε_{stall} (x) = X (x) Y (x) .

(38)

III. EXAMPLE APPLICATIONS

We shall now consider some exemplar classes of models to: provide examples of how to utilize the methods; validate their accuracy; and to show the types of quantities and information that may be extracted.

A useful initial classification of these models is into those which we shall call balanced. We shall refer to a model as being balanced if its petals [see Fig. 4(b)] are detailed balanced. Such models are useful baseline checks as their cycles all have zero affinity, meaning no chemical work is done internally and, hence, the only contributions to chemical work are from external cycles. Furthermore, these models exhibit a proper equilibrium at the stall point and allow for equilibrium arguments to validate the method at this point. It is worth noting that although related to the notion of detailed balanced, the full model with its infinite state space is not detailed balanced.

A. Stalling behavior in a polymerization model with no neighbor–neighbor interactions

We shall start with the simplest case, where the propensities in the coarse-grained model only depend on the monomer type being added/removed: ω_±yx = ω_±y, such as in a simple model for templated self-assembly, shown in Fig. 1(b). Assume that there exists a backbone free energy ΔG_pol controlling the driving as described in Sec. II B 5. Any spanning tree in Λ^± must involve at least one incidence of ΔG_pol, since it appears in every external cycle. Therefore, we can split the ratio of propensities as follows:

\frac{ω_{+ y}}{ω_{- y}} = e^{Δ G_{y}} e^{Δ G_{pol}},

(39)

where ΔG_y encompasses the rest of the details about the models. We note, in general, that ΔG_y may be a function of ΔG_pol, however, in many cases, it is not. These cases include when there is only one completion reaction [highlighted in red in Fig. 4(a)] that contains the dependence on ΔG_pol or if the model is balanced. We may then interpret −ΔG_y as an effective binding free energy of monomer y. If we think of ΔG_pol as the free energy drive of the model away from stalling, we look for a threshold value ΔG_pol = Γ above which the model will not stall. Using Eq. (37), we see that

Γ = - \ln (\sum_{y} e^{Δ G_{y}}) = - \ln Z,

(40)

where $Z$ is the partition function for a system with one state for each monomer type, each state labeled by y and with free energy −ΔG_y. Furthermore, using Eq. (38), the bulk frequencies at the stall point may be written as

ε_{stall} (y) = \frac{e^{Δ G_{y}}}{\sum_{x} e^{Δ G_{x}}} = \frac{1}{Z} e^{Δ G_{y}},

(41)

which is the probability of selecting a state y with free energy −ΔG_y as predicted by equilibrium statistical mechanics. In these results, −ΔG_y looks like the equilibrium contribution to free energy, and the results follow fairly directly in equilibrium. However, these results hold even if the process involves fuel-consuming cycles: entropy may still be being produced at stall. In such cases, the effect of breaking the equilibrium will be to change the effective free energies of selecting a given monomer type.

B. Balanced models of templated polymerization with autonomous separation

Next, we shall consider a class of models where the ratio of propensities may be written as

\frac{ω_{+ y x}}{ω_{- y x}} = e^{Δ G_{y}} e^{- Δ G_{x}} e^{Δ G_{pol}} .

(42)

As before, ΔG_pol, coming from the polymerization reactions represents the driving of this process. Such a class of models includes, most notably, balanced models of templated polymerization with autonomous separation.²³ In these cases, the breaking of the previous copy-template bond every time a new bond is formed enforces the structure in Eq. (42). We shall assume, as in Ref. 23, that ΔG_y is independent of ΔG_pol.

Using Eqs. (37) and (38), we find the stall point to be ΔG_pol = Γ = −ln M and the bulk frequencies at stall $ε_{stall} (y) = \frac{1}{M}$ ⁠, where M the number of monomer types. Physically, we can understand these results by considering balanced models of templated polymerization with autonomous separation. For such models, by definition, there is no entropy production in internal cycles and therefore, the stall point must be thermodynamic equilibrium. In such models, the only driving comes from the polymerization ΔG_pol and the entropic effect having M monomers to choose. These two effect balance at equilibrium.⁶⁴

Next, let us consider the limit that the completion reactions highlighted in red in Fig. 4(a) are much slower than the other reactions. Explicitly, let k be some rate constant at the same order of magnitude of the rates of the process that are not the rates for the completion transitions indicated in red in Fig. 4(a). Write the completion rates as $k_{com} R_{com}^{+} (y, x)$ ⁠, where k_com ≪ k is a rate constant controlling the overall speed of the completion reactions and $R_{com}^{+} (y, x)$ provides any sequence dependence. Similarly, the reverse transitions along the completion edges have the rate $k_{com} R_{com}^{-} (y, x)$ ⁠. Furthermore, let there be n_com such completion reactions in a given petal of the step-wise process (we shall assume this number is the same for all pairs of monomers x, y).

Assume for simplicity that all completion reactions, $R_{com}^{\pm} (y, x)$ ⁠, take the same form in a given petal. Then, we can write the sum over spanning trees Q(y, x) as

Q (y, x) = \frac{1}{n_{com}} \frac{Λ^{-} (y, x)}{k_{com} R_{com}^{-} (y, x)} + O (\frac{k_{com}}{k}),

(43)

since Λ⁻(y, x) has first order terms in k_com/k. This fact can be seen from noting that the leading order terms in Q(y, x) are the trees with no completion reactions and the leading order terms in Λ⁻(y, x) are those same leading order trees of Q(y, x), except with one completion reaction added in. There are n_com such completion reactions and each adds the same leading order term to Λ⁻(y, x). With Q(y, x) taking this form, and remembering Eq. (42), the propensities take the following form:

\begin{aligned} ω_{+ y x} = n_{com} k_{com} R_{com}^{-} (y, x) e^{Δ G_{y} - Δ G_{x} + Δ G_{pol}} + O {(\frac{k_{com}}{k})}^{2}, \\ ω_{- y x} = n_{com} k_{com} R_{com}^{-} (y, x) + O {(\frac{k_{com}}{k})}^{2} . \end{aligned}

(44)

The n_comk_com term cancels in ratios of ω_±yx variables and, therefore, does not affect the sequence statistics. Thus, in the slow completion limit, such models are only affected by the binding free energy differences (ΔG_y − ΔG_x), the driving (ΔG_pol), and the nature of the final completion step $(R_{com}^{-})$ ⁠. Therefore, the finer details do not affect the statistics of the polymers.

Assuming that all completion edges are associated with the same free energy change −ΔG_pol, so that $R_{com}^{-} (y, x) = e^{- Δ G_{pol}}$ ⁠, we may solve for the statistics explicitly. For the case of two monomer types, M = 2, we find the bulk frequency to be ( Appendix I):

\begin{matrix} ε (1) & = & (1 - \frac{1}{2} (e^{- Δ G_{pol}} - 1) (e^{- D G} - 1)) \\ + {(\frac{1}{2} \sqrt{{(e^{- Δ G_{pol}} - 1)}^{2} {(e^{- D G} - 1)}^{2} + 4 e^{- D G}})}^{- 1}, \end{matrix}

(45)

where DG = ΔG₁ − ΔG₂. This expression is plotted in Fig. 5 for DG = 4. From this expression, we can confirm explicitly by substituting in the stall driving, ΔG_pol = −ln 2, that the bulk frequency, indeed, becomes $ε (1) = \frac{1}{2}$ ⁠. Furthermore, taking the irreversible limit, ΔG_pol → ∞, we find that the bulk frequency becomes

ε (1) = \frac{e^{Δ G_{1}}}{e^{Δ G_{1}} + e^{Δ G_{2}}},

(46)

which is the equilibrium statistical mechanics probability of choosing state 1 with free energy −ΔG₁, given state 2 has free energy −ΔG₂. Since the completion reactions are slow and irreversible, in this limit, the process of selecting the monomers is allowed to equilibrate. Therefore, copolymerization is simply sampling from the equilibrium distribution of this process and, hence, tends to the result predicted by equilibrium statistical mechanics.

FIG. 5.

View large Download slide

Plots of the frequency of the less stably-bound (incorrect) monomer with smallest binding free energy, labeled 2 for on- and off-rate discrimination balanced models with k_com = 100 and k_com → 0. The binding free-energy difference for these models is DG = ΔG₁ − ΔG₂ = 4. The models are topologically the Hopfield model as in Fig. 6(a), with ΔG_act = 0 and M_in = M_act = 1. However, for the on-rate discrimination, the free-energy terms are in the binding reactions instead of the unbinding ones. The specific models are given in Appendix J.

Equation (44) shows that in the slow completion limit, the finer details of the reaction network leading to the selection of a specific monomer becomes unimportant and the models collapse onto a single accuracy curve determined by DG, ΔG_pol, and $R_{com}^{-}$ ⁠. Conversely, if we fix all parameters except k_com, we seem to see that the bulk frequencies will tend monotonically to their limits as k_com/k → 0, either from above or below.

We can use this fact to compare bulk frequencies for certain types of models. For example, we may compare on-rate discrimination,²⁶ where incorrect monomers bind more slowly, to off-rate discrimination,²⁶ where incorrect monomers unbind more quickly. An example model comparing on-rate and off-rate discrimination is plotted in Fig. 5 for a model defined in Appendix J. Consider the bulk frequency of an incorrect monomer. On-rate discrimination benefits from fast polymerization and, therefore, tends to its slow polymerization limit from below, whereas off-rate discrimination benefits from allowing the process selecting monomers to equilibrate and, hence, tends to its slow copolymerization limit from above. This fact sets up a hierarchy for a given set of parameters and moderate or strong driving; for the bulk frequency of incorrect monomers, off-rate discrimination > slow copolymerization > on-rate discrimination. This observation is consistent with the results of Sartori and Pigolotti²⁶ and Poulton et al.²³ for kinetic (on-rate) and energetic (off-rate) discrimination.

C. Hopfield’s kinetic proofreading in a model of templated copying with autonomous separation

For our final example, we shall consider an explicit model of copolymerization, with Hopfield’s kinetic proofreading mechanism incorporated into a templated copolymerization system with autonomously separating product in a thermodynamically valid way. From this setup, we can provide a fully worked example of an explicit model and demonstrate the power of the method for analyzing sequences of models with recursive structures as we look at a generalized version of Hopfield’s proofreading incorporated into a model of templated polymerization with autonomous separation.

Explicitly, we first consider the one-loop model of kinetic proofreading shown in Fig. 6(a). There are two monomer types, the right ones x = r and wrong ones x = w. Note that we have already transformed the model so that the sequence of the copy is defined relative to that of the template.³⁵ These monomer types exist in inactive and active states with concentrations M_in and M_act, respectively, relative to some reference concentration, with each monomer type having the same concentration. As previously, we shall assume the environment is sufficiently large such that these concentrations remain constant.

FIG. 6.

View large Download slide

Reaction rates of the (a) 1-loop and (b) N-loop Hopfield kinetic proofreading models implemented in a templated polymerization model with autonomous separation. Each of these subfigures represents a single petal of the step-wise process as in Fig. 4, going from completed state &x → &xy. In both cases, the template is represented by red squares. In (a) the inactive monomer is represented by a white circle and the activated monomers by a dark blue circle. In (b), different levels of activation are represented by increasingly dark shades of blue circles. Furthermore, in (b) the numbers by the states represent the activation level of the monomer. In each case, the desired pathway is highlighted with red arrows.

The monomers may bind to the template either in an active or inactive state with binding free energies −ΔG_x for monomer type x. Inactive monomers may be activated on the template with a free-energy change of ΔG_act. Finally, active monomers may be polymerized into the copolymer chain, with a free-energy change −ΔG_pol. Subsequently, the penultimate monomer of the copolymer unbinds from the template. Each of these reactions is assigned a forward and reverse reaction rate consistent with the thermodynamic model; the full model is illustrated in Fig. 6(a). Conceptually, the proofreading motif functions by providing two opportunities to reject the unwanted monomer w: first, when the unactivated monomer binds and second, after it has been activated. To be effective, a non-zero affinity is required to drive the system around the cycle of states in the correct order: unbound template site → unactivated monomer bound → activated monomer bound.^20,43 We emphasize that this model differs from Hopfield’s original description in two important ways: first, we consider a full, microscopically reversible polymerization process, rather than a single incorporation step with irreversible polymerization; and second, we embed the proofreading motif into a non-trivial polymerization process involving autonomous detachment from the template.

Given the model as described in Fig. 6(a), we first identify the propensities ω_xy connecting completed states. Due to the petal-like structure, we can follow Eq. (13) and simply consider spanning trees of the petal sub-processes illustrated in Fig. 6; Λ⁻(y, x) rooted at &x, Λ⁺(y, x) rooted at &xy, and Q(y, x) for a petal connecting &x and &xy. Explicitly writing out the sums of spanning trees, we obtain

Λ_{1}^{+} (y, x) = [k_{1} k_{act} M_{in} + k_{K P} M_{act} (k_{1} e^{- Δ G_{y}} + k_{act})] k_{pol} e^{- Δ G_{x}},

(47)

\begin{align} Λ_{1}^{-} (y, x) = & [k_{1} k_{act} e^{Δ G_{act} - Δ G_{y}} + k_{K P} e^{- Δ G_{y}} (k_{1} e^{- Δ G_{y}} + k_{act})] \\ \times k_{pol} e^{- Δ G_{pol}}, \end{align}

(48)

\begin{align} Q_{1} (y, x) = & [k_{1} k_{act} e^{Δ G_{act} - Δ G_{y}} + (k_{K P} e^{- Δ G_{y}} + k_{pol} e^{- Δ G_{x}})) \\ (\times (k_{1} e^{- Δ G_{y}} + k_{act})] . \end{align}

(49)

Here, we add a subscript 1 to denote these quantities as applying to the simple, “1-loop,” Hopfield model, which we shall extend to allow more loops later. We note that the ratio Λ⁺(y, x)/Λ⁻(y, x) factorizes as in Eq. (36), and so, we can easily write down the stall condition as ΔG_pol = Γ with

\begin{matrix} Γ & = & - \ln (\frac{k_{1} k_{KP} M_{act} e^{- Δ G_{r}} + k_{act} (k_{KP} M_{act} + k_{1} M_{in})}{k_{1} k_{KP} e^{- Δ G_{r}} + k_{act} (k_{1} e^{Δ G_{act}} + k_{KP})}) \\ + (\frac{k_{1} k_{KP} M_{act} e^{- Δ G_{w}} + k_{act} (k_{KP} M_{act} + k_{1} M_{in})}{k_{1} k_{KP} e^{- Δ G_{w}} + k_{act} (k_{1} e^{Δ G_{act}} + k_{KP})}) . \end{matrix}

(50)

Note that by setting M_in = M_act = 1, ΔG_act = 0 in Eq. (50), Γ collapses to −ln 2 as these conditions reduce the system to a balanced one with a stall point at equilibrium, as in Sec. III B.

The frequency of the right and wrong monomers ɛ(x = r, w) may be calculated from Eq. (18) (the calculation is implemented in the supplementary material). We plot copying error, as represented by ɛ(w), in Fig. 7(a), and demonstrate that it agrees well with the results found from a Gillespie simulation⁶⁵ of the same model. We also compare to a “0-loop” version of the model, in which the inactivated monomers and the inactivated monomer bound state are omitted. As can be seen, the proofreading motif generally improves accuracy when driven above its stall point ΔG_pol = Γ. Indeed, we may write down expressions for the bulk frequency in the irreversible limit (ΔG_pol → ∞) using Eq. (34). In this irreversible limit, we recover Hopfield’s classic argument by taking some further limits consistent with his analysis. Namely, letting M_act, k_act, k_pol → 0, we find $ε (w) / ε (r) = e^{2 (Δ G_{w} - Δ G_{r})}$ ⁠. In this limit, the ratio of incorrect monomers to correct ones involves the square of the binding free energy difference, reflecting the fact that two steps of discrimination have occurred.

FIG. 7.

View large Download slide

Analytical method applied to a 1-loop proofreading model [Fig. 6(a)], compared to Gillespie simulation of the same model and a simpler 0-loop model. For these data, the following parameters were used: ΔG_r = 2, ΔG_w = −2, ΔG_act = −1, M_in = 1, M_act = 0.01, and k₁ = k_act = k_KP = 1. The stall point Γ is marked on each of the plots. The Gillespie simulations used a template of length 2000 and were run until completion with the first monomer being chosen as either r or w with probability 0.5. The statistics were averaged over 2000 copolymers per data point. The chemical work was calculated from the simulation as (Inactive monomers)*(ΔG_act − ln(M_act/M_in)) + L*(ΔG_pol + ln M_act), where “Inactive monomers” is the number of inactive monomers taken out of the environment and L is the length of the template.

We may also write down expressions for the expected chemical work done per net step of the process. This quantity will involve the total current to absorbing states of the step-wise process for starting with a copolymer &xy, which we may write as

\begin{matrix} J_{Tot} (y, x) = \frac{1}{N (y, x)} (Λ_{1}^{+} (r, y) Q (w, y) Q (y, x)) \\ (+ Λ_{1}^{+} (w, y) Q (r, y) Q (y, x) + Λ_{1}^{-} (y, x) Q (r, y) Q (w, y)), \end{matrix}

(51)

where $N$ is a normalization factor that will cancel out of calculations. In order to track each of the terms here, we shall break down the contributions to the chemical work done into three parts, one for each of the petals present in the step-wise process. These three petals correspond to adding a monomer type r, adding a monomer type w, or removing a monomer type y. Let us label each of these contributions to the chemical work with a subscript, $G_{r} (y, x)$ for the transition &xy → &xyr, $G_{w} (y, x)$ for the transition &xy → &xyw, and $G_{q} (y, x)$ for the transition &xy → &x. From the r petal, we have

\begin{matrix} G_{r} (y, x) & = & [(- Δ G_{act} + \ln \frac{M_{i n}}{M_{act}}) k_{1} k_{act} k_{K P} e^{- Δ G_{r}} (M_{i n} + M_{act} e^{Δ G_{act}})) \\ + ((Δ G_{pol} + Δ G_{r} - Δ G_{y} + \ln M_{i n} - Δ G_{a}) (k_{1} k_{act} k_{pol} M_{i n} e^{- Δ G_{y}})) \\ + ((Δ G_{pol} + Δ G_{r} - Δ G_{y} + \ln M_{act}) k_{K P} k_{pol} M_{act} e^{- Δ G_{y}}) \\ \times ((k_{1} e^{- Δ G_{r}} + k_{act})] \frac{Q (w, y) Q (y, x)}{N (y, x) J_{Tot} (y, x)} . \end{matrix}

(52)

The first line of Eq. (52) corresponds to the chemical work associated with the internal cycle (inactive monomer binds, gets activated, and the activated monomer unbinds). The second line corresponds to an external cycle: an inactive monomer binds to the template, is activated, and is polymerized into the chain with the previous monomer y, detaching from the template. The third line corresponds to the alternative external cycle: an active monomer binds to the template and is polymerized with monomer y, unbinding from the template. We may similarly write down $G_{w} (y, x)$ as Eq. (52), except swapping r and w. Finally, the contribution to the chemical work from the petal for removing monomer y may be written as

\begin{matrix} G_{q} (y, x) & = & [- (Δ G_{pol} + Δ G_{y} - Δ G_{x} + \ln M_{i n} - Δ G_{a})) \\ \times (k_{1} k_{act} k_{pol} e^{Δ G_{act} - Δ G_{y} - Δ G_{pol}}) - ((Δ G_{pol} + Δ G_{y})) \\ ((- Δ G_{x} + \ln M_{act}) k_{K P} k_{pol} e^{- Δ G_{pol} - Δ G_{y}} (k_{1} e^{- Δ G_{r}} + k_{act})] \\ \times \frac{Q (r, y) Q (w, y)}{N (y, x) J_{Tot} (y, x)} . \end{matrix}

(53)

Here, only external cycles are possible. The first line corresponds to monomer x rebinding to the template, monomer y being depolymerized, this monomer being deactivated and an inactive monomer y unbinding from the template; and the second line to x rebinding, y being depolymerized, and an active monomer y unbinding from the template. The distribution ξ(y, x) may be calculated from Eq. (20) and P from Eq. (21) (both demonstrated in the supplementary material), letting the chemical work done per net step of the 1-loop model be written as

Δ G = \frac{1}{2 P - 1} \sum_{x, y \in {r, w}} ξ (y, x) (G_{r} (y, x) + G_{w} (y, x) + G_{q} (y, x)) .

(54)

This chemical work done is plotted for a certain set of parameters in Fig. 7(b) and is also compared to both the results of direct simulation and the simpler “0-loop” model, which has chemical work ΔG_pol. The free-energy cost of the proofreading mechanism diverges as ΔG_pol → Γ since there will be a finite chemical work done per monomer addition/removal step due to the proofreading internal cycle and the number of addition/removal steps per net step diverges. Furthermore, for large ΔG_pol, the work tends to be dominated by ΔG_pol, albeit very slowly, as shown by the orange line gradually approaching ΔG_pol (the blue line) in Fig. 7(b).

In addition, we can find an expression for the time taken per net step forwards, Eq. (22). For this quantity, we need the explicit expression for the normalization, $N$ ⁠. Similar to the chemical work, we can split this term into contributions from the petal adding an r, $N_{r} (y, x)$ ⁠; from the petal adding a w, $N_{w} (y, x)$ ⁠; from the petal removing monomer y, $N_{q} (y, x)$ and a contribution from the central node. These normalization terms come from the sums of the spanning trees directed to the individual nodes in the closed step-wise process. We see that

\begin{matrix} N_{r} (y, x) & = & [k_{1} M_{in} (k_{act} e^{Δ G_{act}} + k_{K P} e^{- Δ G_{r}} + k_{pol} e^{- Δ G_{y}})) \\ + k_{K P} k_{act} M_{act} (e^{Δ G_{act}} + k_{1} k_{act} M_{in} + k_{K P} k_{act} M_{act}) \\ (+ k_{1} k_{K P} e^{- Δ G_{r}}] Q (y, x) Q (w, y), \end{matrix}

(55)

with a similar result for $N_{w} (y, x)$ except swapping r and w. Finally, for the monomer removal petal, we have

N_{q} (y, x) = k_{pol} e^{- Δ G_{pol}} (k_{1} e^{- Δ G_{y}} + k_{act} + k_{act} e^{Δ G_{act}}) Q (r, y) Q (w, y) .

(56)

The total normalization is, then,

N (y, x) = N_{r} (y, x) + N_{w} (y, x) + N_{q} (y, x) + Q (y, x) Q (r, y) Q (w, y),

(57)

with the last term being the contribution from the starting, central node. This normalization can be used in Eq. (51) to give the current to absorbing states, which can be used in Eq. (22) to find the expected time per net step. This time is plotted in Fig. 7(c), alongside a simulation of the same model and the simplified 0-loop model for comparison. Like the chemical work in Fig. 7(b), the time per net step diverges as ΔG_pol → Γ, since each monomer addition/removal step will take finite time, but the number of such steps required for a net forward step diverges. Unsurprisingly, the time taken for a given driving for the Hopfield model is longer than that of the simple model, due to the proofreading cycle.

Hopfield’s model for proofreading may be naturally extended to include N activation stages instead of just one.^48,50 We shall call these extensions the N-loop Hopfield models. These models can be solved recursively to write down expressions for the sums over spanning trees $Λ_{N}^{\pm} (y, x), Q_{N} (y, x)$ as a function of the number of loops N. We shall consider the model as in Fig. 6(b). A detailed derivation of the sums over the spanning trees is given in Appendix K. From these sums over spanning trees, we calculate the bulk frequencies, the time taken per net step, and the chemical work done per net step using recursive relations (see Appendix K).

For simplicity, we shall discuss the case where the monomer binding free energy is only dependent on monomer type, not on the activation stage; each activation stage is associated with a free energy change of ΔG_act; each active monomer is present in the environment at a concentration M_act except for the inactive monomers at concentration M_in; and the overall rate constants are k₁ for binding of inactive monomer, k_KP for binding of active monomers, and k_act for activation of monomers. Under these assumptions, the corresponding rates are given in Appendix K.

To reduce the frequency of incorrect monomers in the product, we wish to have a low concentration M_act of active monomers in solution to force the system into utilizing the proofreading cycles. Indeed, the bulk error probability in the irreversible limit [calculated using Eq. (34) and plotted in Fig. 8(a)] shows a strong improvement with loop number for low M_act, but larger values of M_act lead to a much worse performance and limited (or negative) returns to increasing the number of loops.

FIG. 8.

View large Download slide

Many-loop models have a limited efficacy for finite M_act for the proofreading model introduced in Fig. 6. Plots of the error in the irreversible limit ɛ_irrev(w) and the stall point driving Γ, for different values of the active monomer concentrations M_act and other parameters: ΔG_r = 2, ΔG_w = −2, M_in = 1, ΔG_act = −1, and ks = 1. The N = 0 error is not shown for clarity, but is 0.5 for all M_act.

However, for finite driving strength ΔG_pol, we cannot allow this concentration to be arbitrarily small. To see why, consider the stall point Γ(N) derived in Appendix K and plotted for a certain set of parameters in Fig. 8(b). It is observed that the stall point driving increases monotonically with N and that this increase is faster and tends to a higher limit for smaller M_act. We find that the limiting Γ scales approximately linearly with −ln(M_act). Intuitively, introducing more monomer states at a low concentration in the environment destabilizes the polymer. For a small M_act and driving ΔG_pol, the depolymerization of the polymer into these activated states competes with its tendency to grow by binding to and activating the inactive monomers.

One drawback of proofreading with a large number of loops is, therefore, that the tendency to disassemble the growing polymer increases. A second effect is a tendency to introduce errors by alternate pathways if M_act is non-zero. In particular, for M_act ≠ 0, we observe in Fig. 8(a) a minimum in ɛ_irrev(w) for a relatively small value of N. This minimum can be explained by splitting the pathways by which a monomer can go from solution to being incorporated into the polymer into two; either starting from a fully inactive monomer or from a partially activated one. The pathway starting with an inactive monomer will have the highest discrimination between the right and wrong monomers and will improve exponentially with more loops, as demonstrated by the exponential decrease in error for M_act = 0. However, the probability that a monomer taking this path will reach polymerization falls exponentially with loop number at the same time. On the other hand, the pathway from partially active monomers will give an error that reaches some non-zero limit as the number of loops N increases. Furthermore, the rate with which activated monomers bind to an available template site and subsequently get incorporated into the polymer will also tend to a constant. As such, the error will initially decrease exponentially with N, but for non-zero M_act, will eventually become dominated by the less discriminating, partially active monomer pathways through which monomers are more likely to be incorporated into the polymer.

Having calculated the error probability ɛ(w) at finite driving, plotted in Fig. 9(a); having used ɛ(x, y) to calculate the entropy rate; and having calculated $Δ G$ ⁠, we can evaluate the efficiency η, as in Eq. (26) (see the supplementary material for demonstrations). This efficiency is plotted in Fig. 9 for N = 0, 1, 5, 10 and a certain set of parameters. Although accuracy is generally increased above the stall point, we see that in this particular model, kinetic proofreading requires much more work than the minimum required to generate information and as such is inefficient. In addition, the gradient of the efficiency at minimum driving Γ_N is zero for N > 0, reflecting how at minimum driving, the number of monomer addition/removal steps diverges, but the chemical work done per such step remains finite.

FIG. 9.

View large Download slide

Plots of (a) the error and (b) efficiency of the N-loop proofreading model (Fig. 6) for a range of N, with the same parameters as in the one-loop Hopfield case, shown in Fig. 7, as a function of driving ΔG_pol. Proofreading is observed to generally increase the accuracy above its stall point, but in a thermodynamically inefficient way. The enhanced plot in the second graph shows the efficiencies near the stall point for each of the loop numbers on a non-logarithmic scale to emphasize the decreasing gradient at stall.

IV. CONCLUSION

We have presented a method for analyzing copolymerization models with complex networks of reactions leading to the incorporation or removal of monomers. By coarse graining, a model may be transformed into a simpler model that may be solved and, thereafter, information from the fine-grained process may be put back into the model to extract thermodynamic or kinetic quantities such as chemical work done, molecule exchange, or time taken. The approach allows for complex incorporation motifs to be considered alongside the nearest neighbor interactions in a thermodynamically well-defined model of polymerization with microscopic reversibility. We note that all of these features were present in the kinetic proofreading example in Sec. III C. Moreover, phenomena such as the shift in stall point with loop number and the non-monotonicity of the error rate with the loop number rely on these features being present in the model.

In general, this method provides a way to extract model predictions numerically quickly and without the need for simulations. Doing so is particularly useful when simulating polymer growth is slow, either due to the details of the incorporation process or because the polymer is near its stall point. In addition, the approach makes screening of a large parameter space for a given model topology feasible.

In addition to the numerical performance, this approach allows for analytic results in simpler models or those with helpful symmetries and in certain limits for more complex models. The process of summing over the spanning trees is particularly well suited to identifying the structure of the process and providing simplified results.

Moving forward, it is an open question as to whether components of the techniques developed here can be applied outside of the context of infinitely long polymers whose tips have reached a steady state. An obvious goal would be a simplified way to analyze finite-length “oligomers.”²⁴ More generally, we believe that the key equation of this paper, Eq. (14), may be applied more generally for the coarse graining of Markov processes and, in particular, that if a set of states are enclosed between two boundary states, in the sense that any path from one of the trapped states to outside must pass through one of the boundary states, then this set of states may be replaced by a pair of edges analogously to Eq. (14) that shall preserve the steady state properties of the Markov process.

This framework could be applied to explore models of copolymerization processes such as those presented in Refs. 17–19, 22, 23, 27, 28, 30–35, 39–41, and 45–51, more straightforwardly or more thoroughly. Alternatively, the method would allow for more complex reaction steps to be included in such models. The framework presented here is particularly useful when backward steps are relevant, either when the system is weakly driven and, thus, operating near the stall point or when thermodynamics is of importance or interest. We also predict that it will be useful to guide design principles for synthetic copolymerization systems, which are often particularly well-described by the class of models studied here.

SUPPLEMENTARY MATERIAL

The supplementary material contains a C++ script implementing the Gillespie algorithm that reproduces the data for the 1-loop Hopfield kinetic proofreading model presented in Fig. 7, and a MATLAB script for numerically calculating quantities of the 1-Loop and N-Loop Hopfield kinetic proofreading models presented in Sec. III C and shown in the solid lines of Fig. 7 and the points of Figs. 8 and 9.

ACKNOWLEDGMENTS

This work was part of a project that has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Program (Grant Agreement No. 851910). T.E.O. was supported by a Royal Society University Fellowship. J.J. was supported by a Royal Society Ph.D. studentship.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

All authors conceived of the project. B.Q. produced the methodology and analysis and wrote the initial draft. All authors interpreted results and reviewed and edited this paper.

Benjamin Qureshi: Conceptualization (equal); Formal analysis (lead); Methodology (equal); Writing – original draft (lead). Jordan Juritz: Conceptualization (equal); Visualization (equal); Writing – review & editing (equal). Jenny M. Poulton: Conceptualization (equal); Writing – review & editing (equal). Adrian Beersing-Vasquez: Conceptualization (equal); Writing – review & editing (equal). Thomas E. Ouldridge: Conceptualization (equal); Funding acquisition (lead); Supervision (lead); Writing – review & editing (equal).

DATA AVAILABILITY

The data that support the findings of this study are openly available in Zenodo at https://doi.org/10.5281/zenodo.7271702.

APPENDIX A: FACTORIZING SUMS OF SPANNING TREES

We note here that sums of spanning trees can be factorized in terms of self-avoiding walks (SAWs), a result which is both useful for generating sets of spanning trees and allows us to make statements about ratios of propensities of balanced models. For a given process, $G = (X, K)$ ⁠, for which we wish to find the sum of spanning trees rooted at $x_{1} \in X$ ⁠, we may factorize this sum in terms of self-avoiding walks (SAWs) between two vertices in the graph. Select some other arbitrary vertex $x_{2} \in X / {x_{1}}$ and let $S (x_{2}, x_{1})$ be the set of SAWs from x₂ to x₁. For each $S \in S (x_{2}, x_{1})$ ⁠, we can construct $G_{S} = ({s} \cup (X / S), K_{S})$ analogously to Eq. (9), whereby we collapse the nodes S in the SAW into the single node s. The sum over spanning trees rooted at x₁ may then be written as

\sum_{T \in T (x_{1})} \prod_{e \in T} K (e) = \sum_{S \in S (x_{2}, x_{1})} \underset{SAW term}{\underset{︸}{[\prod_{e \in S} K (e)]}} \underset{Spanning tree term}{\underset{︸}{[\sum_{T \in T_{S} (s)} \prod_{e \in T} K_{S} (e)]}},

(A1)

where $T (x), T_{S} (x)$ are the sets of spanning trees directed to x for the original process $G$ and the new process $G_{S}$ ⁠. For example, in Fig. 3(a), the spanning trees are arranged in terms of SAWs from node 1 to node 3, with the first row for SAW: 1 → 2 → 4 → 3; the second row for 1 → 2 → 3; and the last three rows for 1 → 3. Similarly, for Fig. 3(b), the trees are arranged in terms of SAWs from node 1 to node 4 with row one for 1 → 2 → 3 → 4; row two for 1 → 2 → 4; row three for 1 → 3 → 4; and row four for 1 → 3 → 2 → 4.

APPENDIX B: NORMALIZATION CONSTANT FOR EXAMPLE ABSORBING MARKOV PROCESS

The normalization constant for the closed example process, given in Fig. 2(b), can be found by considering the spanning trees rooted at each of the nodes. Factorizing these in terms of SAWs, we write

\begin{matrix} N & = & [r_{34} r_{42} r_{21} + r_{32} r_{24} k_{B} + r_{32} r_{21} (r_{43} + r_{42} + k_{B}) + r_{34} k_{B} (r_{24} + r_{23} + r_{21})) \\ + ((r_{31} + k_{A}) (r_{21} r_{43} + r_{21} k_{B} + r_{42} r_{21} + r_{23} r_{43} + r_{23} k_{B} + r_{42} r_{23} + r_{24} r_{43} + r_{24 k_{B}})] \\ + [r_{13} r_{34} r_{42} + r_{13} r_{32} (r_{43} + r_{42} + k_{B})) \\ + (r_{12} (r_{34} r_{42} + r_{43} r_{32} + r_{32} k_{B} + r_{32} r_{42} + r_{43} (r_{31} + k_{A}) + k_{B} (r_{31} + k_{A}) + r_{42} (r_{31} + k_{A}) + r_{34} k_{B})] \\ + [r_{12} r_{24} r_{43} + r_{12} r_{23} (r_{42} + r_{43} + k_{B}) + r_{13} ((r_{43} + k_{B}) (r_{21} + r_{23} + r_{24}) + r_{42} (r_{21} + r_{23}))] \\ + [r_{13} r_{34} (r_{21} + r_{23} + r_{24}) + r_{13} r_{32} r_{24} + r_{12} r_{23} r_{34} + r_{12} r_{24} (r_{31} + k_{A} + r_{32} + r_{34})] . \end{matrix}

(B1)

The first square bracket corresponds to the trees rooted at node 1, organized by SAWs from node 3; the second, to trees rooted at 2 organized by SAWs from 1; the third, to trees rooted at 3 organized by SAWs from 1; and the fourth, to trees rooted at 4 organized by SAWs from 1.

APPENDIX C: EQUIVALENCE BETWEEN CHEMICAL WORK CALCULATED FROM EDGES AND CYCLES

Here, we shall show the equivalence of chemical work for a process calculated by summing over edges vs summing over cycles. For this comparison, consider a process $(X, K)$ without any absorbing states (for simplicity) and such that every edge is microscopically reversible, and let π(x) be the steady state probability to be in state x. For an edge x ⇋ y, as described in Sec. II A 4, the net current through this edge is

J_{x ⇋ y} = π (x) K (x, y) - π (y) K (y, x) .

(C1)

We can write π(x) in terms of the spanning trees by MCTT, and by Appendix A, we may expand the sum over the spanning trees by SAWs from y to x. For π(y), we may expand by SAWs from x to y such that the spanning tree terms of both expansions are the same and only the direction of the edges in the SAW terms is flipped. The net current may then be written as

\begin{align} J_{x ⇋ y} = & \frac{1}{N} \sum_{S \in S (y, x)} [K (x, y) \prod_{e \in S} K (e) - K (y, x) \prod_{e^{'} \in S} K (e^{'})] \\ \times [\sum_{T \in T_{S} (s)} \prod_{e \in T} K_{S} (e)], \end{align}

(C2)

where $S (x, y)$ is the set of SAWs from node x to node y; $N$ is the normalization as in Eq. (2), and e′ is the edge in the opposite direction, i.e., if e = x → y, e′ = y → x; and the last bracketed term is the spanning tree part for SAW, S, as in Eq. (A1). One of the SAWs from y to x will simply be the single transition x → y, however, this term will cancel out from the sum leaving just the non-trivial SAWs. Taking a non-trivial SAW from y to x and multiplying by the rate K(x, y) gives a cycle containing the edge x → y. Therefore, the current may be written as a sum over cycle currents, as in Sec. II A 4, of cycles which contain the edge x → y minus those which contain y → x. Each of the edges contains a contribution to the chemical work $\ln (\frac{K (x, y)}{K (y, x)})$ ⁠. The total chemical work before absorption is the sum over all edges of these contributions,

W_{chem} = \sum_{x ⇋ y} \ln (\frac{K (x, y)}{K (y, x)}) \frac{J_{x ⇋ y}}{J_{Tot}} .

(C3)

Since J_x⇋y may be split up as a sum over cycles in this sum, we may collect the parts of this sum corresponding to given cycles and convert the sum over edges into a sum over cycles. By performing these operations, we find the contribution to the chemical work from cycle C to be $\ln (\frac{A (C)}{A (C^{'})})$ ⁠, i.e., the affinities as we might expect. Hence, the sum over cycles is equivalent to the sum over edges.

APPENDIX D: CYCLES OF THE EXAMPLE ABSORBING PROCESS

We divide the cycles of the example process, shown in Fig. 2(a), into internal cycles, external cycles to absorbing state A, and external cycles to absorbing state B. First, the internal cycles are

where the cycle is written out below in the clockwise direction. Similarly, we find the external cycles to state A,

Finally, the external cycles to absorbing state B are

APPENDIX E: NUMBER OF STEPS PER NET FORWARD STEP OF A RANDOM WALK

Here, we shall derive the number of steps per net forward step of a random walk. Let us set up a random walk as follows: let the state space be the nodes {0, 1, …L}, where L is the length of the walk (polymer). Let the transition 0 → 1 have probability 1, i → i + 1 for i = 1, …L − 1 have probability p, i → i − 1 for i = 1, …L − 1 have probability q = 1 − p, and let state L be an absorbing state as in Fig. 10.

FIG. 10.

View large Download slide

Graphical representation of the random walk process considered.

Then, we wish to find the expected number of steps to absorption given we start from state 0, for which we can utilize the spanning tree methods with Eq. (22). Since the total rate out of any state sums to one, the expected number of steps equals the expected time to absorption. Thus, we can form the closed process starting at 0. Let f(n) be the sum over spanning trees rooted at node n for the closed process. f(n) is given by

f (n) = \{\begin{cases} \sum_{i = 0}^{L - 1} p^{i} q^{L - 1 - i} & for n = 0 \\ p^{n - 1} \sum_{i = 0}^{L - 1 - n} p^{i} q^{L - 1 - n - i} & for n = 1, \dots, L - 1 . \end{cases})

(E1)

From this equation, the expected number of steps before absorption is

E [steps] = \frac{\sum_{n = 0}^{L - 1} f (n)}{p f (L - 1)} .

(E2)

By utilizing the formulas for finite geometric series, we can find the expected number of steps to be

E [steps] = \frac{1}{2 p - 1} (L - 1 - \frac{q}{p^{L}} (\frac{p^{L} - q^{L}}{p - q}) + \frac{q^{L}}{p^{L}} + \frac{p^{L} - q^{L}}{p^{L - 1}}) .

(E3)

Most of this expression is sub-linear in L and as such,

\lim_{L \to \infty} \frac{E [steps]}{L} = \frac{1}{2 p - 1},

(E4)

which is the net number of steps per net forward step.

APPENDIX F: THE FREQUENCY AT STALL IS GIVEN BY THE DIAGONAL COFACTORS OF A MATRIX

We wish to show that, at the stall point, the frequency with which a monomer appears in the bulk of the copolymer is proportional to the cofactor of the corresponding diagonal element of a matrix,

ε (x) \propto A_{x x},

(F1)

where A_ij is the cofactor of element i, j of the matrix $1 - Z$ ⁠. To show this relation, we will rely on the relationship between the cofactors and vectors of the nullspace of a matrix. Let M be an arbitrary matrix with a one dimensional nullspace, and let A be its matrix of cofactors. Recall that

M A^{T} = \det (M) 1 = 0 .

(F2)

Thus, any column of A^T is in the nullspace of M. In anticipation, let $\vec{μ}$ be a vector in the nullspace of M and $\vec{v}$ be a vector in the nullspace of M^T. Since M has a one dimensional nullspace, then,

\frac{μ_{x}}{μ_{y}} = \frac{A_{i x}}{A_{i y}},

(F3)

for some arbitrary i. Similarly,

\frac{v_{x}}{v_{y}} = \frac{A_{x j}}{A_{y j}},

(F4)

for arbitrary j.

Looking at Eqs. (15) and (16), noting that near the stall point, v_z ≪ ω_±y,x, we see that the tip probabilities μ(x) form a vector in the nullspace of $1_{M} - Z$ and the tip velocities v_x form a vector in the nullspace of $1_{M} - Z^{T}$ ⁠. Hence, we have that

\frac{μ (y) v_{y}}{μ (x) v_{x}} = \frac{A_{j y} A_{y i}}{A_{j x} A_{x i}},

(F5)

for arbitrary i, j. Thus, we may choose j = y and i = x leading to cancellation such that

\frac{μ (x) v_{x}}{μ (y) v_{y}} = \frac{A_{x x}}{A_{y y}} .

(F6)

Since

ε (x) = \frac{μ (x) v_{x}}{\sum_{y} μ (y) v_{y}} = \frac{A_{x x}}{\sum_{y} A_{y y}},

(F7)

we get the required result.

APPENDIX G: THE FREQUENCIES IN THE IRREVERSIBLE LIMIT ARE GIVEN BY THE STEADY STATE OF A PROCESS OF THE COMPLETE GRAPH

We wish to find an expression for the frequency with which monomer x appears in the bulk of the copolymer in the irreversible limit. This limit is such that the backward propensities ω_−yx = 0. With this assumption, from Eq. (15), we have

v_{x} = \sum_{y} ω_{+ y x} .

(G1)

With this form for the velocities, we may manipulate Eq. (16),

\begin{aligned} μ (x) = \sum_{y} \frac{ω_{+ y x}}{\sum_{z} ω_{+ z x}} μ (y), \\ \sum_{z \neq x} ω_{+ z x} μ (x) + ω_{+ x x} μ (x) = \sum_{y \neq x} ω_{+ x y} μ (y) + ω_{+ x x} μ (x) . \end{aligned}

(G2)

The last line is the equation for the steady state of a Markov process with probability μ(x) to be in state x and rate ω_+yx of transition from state x to state y. Thus, set μ(x) to be the steady state probability distribution of the Markov process on M states with transition rates from state x to y given by ω_+yx and v_x = ∑_yω_+yx. Then, calculating

ε (x) \propto μ (x) v_{x},

(G3)

gives the required result. Finding the distribution μ(x) in terms of the spanning trees of the complete graph on M elements gives Eq. (36).

APPENDIX H: SIMPLIFICATION OF RESULTS FOR FACTORIZABLE RATIOS OF PROPENSITIES

We shall show that, if the ratio of propensities factorizes as in Eq. (38), then we may simplify the stall condition and frequency of monomers at stall. Thinking of the functions X and Y as column vectors, since they have a discrete domain, the matrix Z may be written,

Z = \vec{X} {\vec{Y}}^{T} .

(H1)

By a well-known result,⁷⁰

\det (1_{M} - \vec{X} {\vec{Y}}^{T}) = 1 - {\vec{Y}}^{T} \vec{X} = 1 - \sum_{x} X (x) Y (x) .

(H2)

Here, rearranging gives Eq. (39). As shown, the frequency of monomer x in the bulk of the copolymer is given by the cofactor of the diagonal elements of $1_{M} - \vec{X} {\vec{Y}}^{T}$ ⁠. The cofactor A_xx may be written as

A_{x x} = \det (1_{M - 1} - {\vec{X}}_{[x]} {\vec{Y}}_{[x]}^{T}) = 1 - \sum_{y \neq x} X (y) Y (y) = X (x) Y (x),

(H3)

using the stall condition, where ${\vec{X}}_{[x]}$ is the vector $\vec{X}$ missing element X(x), i.e., ${\vec{X}}_{[x]} = {(X (1), \dots X (x - 1), X (x + 1), \dots X (M))}^{T}$ ⁠. In addition, because of the stall condition ∑_xX(x)Y(x) = 1,

ε_{stall} (x) = X (x) Y (x)

(H4)

is already normalized.

APPENDIX I: FREQUENCY FOR A BALANCED MODEL WITH TWO MONOMER TYPES IN THE SLOW POLYMERIZATION LIMIT

We shall derive the frequency of monomer x in the bulk of the copolymer with propensities given by Eq. (46), canceling n_comk_com, with $R_{com}^{-} = e^{- Δ G_{pol}}$ ⁠, and with M = 2. With these propensities, Eq. (15) becomes

v_{1} = \frac{v_{1}}{e^{- Δ G_{pol}} + v_{1}} + \frac{e^{- D G} v_{2}}{e^{- Δ G_{pol}} + v_{2}}

(I1)

v_{2} = \frac{e^{D G} v_{1}}{e^{- Δ G_{pol}} + v_{1}} + \frac{v_{2}}{e^{- Δ G_{pol}} + v_{2}},

(I2)

where DG = ΔG₁ − ΔG₂. These equations may be solved by the following form the velocities:

v_{x} = e^{Δ G_{y} - Δ G_{x}} v_{y} .

(I3)

Doing so reduces Eq. (15) to a quadratic equation,

0 = e^{D G} v_{1}^{2} + (1 + e^{D G}) (e^{- Δ G_{pol}} - 1) v_{1} + e^{- Δ G_{pol}} (e^{- Δ G_{pol}} - 2)

(I4)

with one positive root,

\begin{align} v_{1} = & \frac{1}{2} ((1 - e^{- Δ G_{pol}}) (e^{- D G} + 1)) \\ (+ \sqrt{{(e^{- Δ G_{pol}} - 1)}^{2} {(e^{- D G} - 1)}^{2} + 4 e^{- D G}}), \end{align}

(I5)

when the system is not stalling. v₂ can be found from in terms of v₁ as v₂ = e^DGv₁. Furthermore, a quick check confirms v₁ = 0 if ΔG_pol = −ln 2. Furthermore, with v_y known, Eq. (16) is a simple linear equation, and μ can be found as the eigenvector of the matrix

(\begin{matrix} \frac{1}{e^{- Δ G_{pol}} + v_{1}} & \frac{e^{D G}}{e^{- Δ G_{pol}} + v_{1}} \\ \frac{e^{- D G}}{e^{- Δ G_{pol}} + v_{2}} & \frac{1}{e^{- Δ G_{pol}} + v_{2}} \end{matrix}),

(I6)

with eigenvalue 1 and normalized to sum to 1. Combining the solutions for μ and v, using Eq. (18), gives Eq. (47).

APPENDIX J: MODEL USED FOR BALANCED ON-RATE VS OFF-RATE DISCRIMINATION COMPARISONS

Figure 11 shows the on-rate and off-rate discrimination models used to produce Fig. 5.

FIG. 11.

View large Download slide

Reaction rates of the (a) off-rate and (b) on-rate discrimination models used to produce the results of Fig. 5. These reactions represent a single petal of the step-wise process (Fig. 4) between the completed states &x and &xy. For the results in Fig. 5 for the off-rate and on-rate curves, the following parameters were used: ΔG₁ = 2, ΔG₂ = −2, k₁ = k_KP = k_act = 1, and k_com = 100.

APPENDIX K: EQUATIONS FOR N-LOOP HOPFIELD MODEL

The sum over the spanning trees of the N-loop model can be written in terms of sums over spanning trees of the lower loop number models. We label the reaction rates for the N-loop process as shown in Fig. 6. The N-loop model has one more node and two more edges than the N − 1-loop model. Let a subscript N denote the sums over the spanning trees for the N-loop models. Tracking the spanning trees, we see,

Λ_{N}^{+} = R_{act}^{+} (N) Λ_{N - 1}^{+} + R_{K P}^{+} (N) \frac{R_{pol}^{+}}{R_{pol}^{-}} \sum_{i = 0}^{N} [\prod_{j = 0}^{i - 1} R_{act}^{+} (N - j)] Λ_{N - 1 - i}^{-},

(K1)

Λ_{N}^{-} = R_{act}^{-} (N) Λ_{N - 1}^{-} + R_{K P}^{-} (N) \sum_{i = 0}^{N} [\prod_{j = 0}^{i - 1} R_{act}^{+} (N - j)] Λ_{N - 1 - i}^{-},

(K2)

Q_{N} = \frac{1}{R_{pol}^{-}} (Λ_{N}^{-} + R_{pol}^{+} \sum_{i = 0}^{N} [\prod_{j = 0}^{i - 1} R_{act}^{+} (N - j)] Λ_{N - 1 - i}^{-}),

(K3)

with the initial conditions

\begin{aligned} Λ_{- 1}^{\pm} = R_{pol}^{\pm}, \\ Λ_{0}^{\pm} = R_{pol}^{\pm} R_{in}^{\pm}, \\ Q_{0} = R_{pol}^{+} + R_{in}^{-} . \end{aligned}

(K4)

The sum-product can be eliminated by subtracting terms proportional to $Λ_{N - 1}^{\pm}, Q_{N - 1}$ ⁠, leaving just

\begin{matrix} Λ_{N}^{+} & = & (R_{act}^{+} (N) + \frac{R_{K P}^{+} R_{act}^{+} (N)}{R_{K P}^{+} (N - 1)}) Λ_{N - 1}^{+} \\ - \frac{R_{K P}^{+} (N) R_{act}^{+} (N) R_{act}^{+} (N - 1)}{R_{K P}^{+} (N - 1)} \times Λ_{N - 2}^{+} + R_{K P}^{+} (N) \frac{R_{pol}^{+}}{R_{pol}^{-}} Λ_{N - 1}^{-}, \end{matrix}

(K5)

\begin{matrix} Λ_{N}^{-} & = & (R_{act}^{-} (N) + R_{K P}^{-} (N) + \frac{R_{K P}^{-} (N) R_{act}^{+} (N)}{R_{K P}^{-} (N)}) Λ_{N - 1}^{-} \\ \times - \frac{R_{K P}^{-} (N) R_{act}^{+} (N) R_{act}^{-} (N - 1)}{R_{K P}^{-} (N - 1)} Λ_{N - 2}^{-}, \end{matrix}

(K6)

Q_{N} = R_{act}^{+} (N) Q_{N - 1} + \frac{Λ_{N}^{-}}{R_{pol}^{-}} + (R_{pol}^{+} - R_{act}^{+} (N)) \frac{Λ_{N - 1}^{-}}{R_{pol}^{-}},

(K7)

with the same initial conditions as above. This system of recursion relations may be used to generate the terms of the spanning tree sums quickly.

In certain simple cases, Eqs. (K5)–(K7) can be solved as a function of N. For example, when the reaction rates are not a function of N, such as

\begin{aligned} R_{in}^{+} = k_{1} M_{in}, \\ R_{in}^{-} = k_{1} e^{- Δ G_{y}}, \\ R_{act}^{+} (n) = k_{act}, \\ R_{act}^{-} (n) = k_{act} e^{Δ G_{act}}, \\ R_{K P}^{+} (n) = k_{K P} M_{act}, \\ R_{K P}^{-} (n) = k_{K P} e^{- Δ G_{y}}, \end{aligned}

(K8)

for n ∈ {1, …N}, for the step-wise process with monomer tip &xy, where k_i are some overall rates, M_in and M_act represent the concentrations of inactive or active monomers, respectively, and ΔG_act represents the chemical work upon moving a monomer up one activation stage. The rates in Eq. (K8) are used for the numeric results in Figs. 8 and 9. In this case, the sums over spanning trees are

\begin{matrix} Λ_{N}^{+} (y, x) & = & k_{pol} e^{- Δ G_{x}} (k_{act}^{N} (k_{1} M_{in} - k_{1} M_{act} + k_{act} e^{Δ G_{y}} M_{act} (e^{Δ G_{act}} - 1))) \\ + (\frac{k_{K P} M_{act} e^{- Δ G_{y}}}{Δ} [\frac{(k_{1} λ_{+} + k_{act} (k_{K P} - k_{1}))}{{(λ_{+} - k_{act})}^{2}} λ_{+}^{N + 1})) \\ ((- \frac{(k_{1} λ_{-} + k_{act} (k_{K P} - k_{1}))}{{(λ_{-} - k_{act})}^{2}} λ_{-}^{N + 1}]), \end{matrix}

(K9)

\begin{matrix} Λ_{N}^{-} (y, x) & = & \frac{k_{pol} e^{- Δ G_{pol}} e^{- Δ G_{y}}}{Δ} [(k_{1} λ_{+} + k_{act} (k_{K P} - k_{1})) λ_{+}^{N}) \\ (- (k_{1} λ_{-} + k_{act} (k_{K P} - k_{1})) λ_{-}^{N}], \end{matrix}

(K10)

\begin{matrix} Q_{N} (y, x) & = & \frac{k_{pol} e^{- Δ G_{x}}}{Δ} [(k_{1} e^{- Δ G_{y}} + k_{act} - λ_{-}) λ_{+}^{N} - (k_{1} e^{- Δ G_{y}})) \\ ((+ k_{act} - λ_{+}) λ_{-}^{N}] + \frac{e^{- Δ G_{y}}}{Δ} [(k_{1} λ_{+} + k_{act} (k_{K P} - k_{1})) λ_{+}^{N}) \\ (- (k_{1} λ_{-} + k_{act} (k_{K P} - k_{1})) λ_{-}^{N}], \end{matrix}

(K11)

where

λ_{\pm} = \frac{1}{2} (k_{act} + k_{K P} e^{- Δ G_{y}} + k_{act} e^{Δ G_{act}} \pm Δ),

(K12)

Δ = \sqrt{{(k_{act} + k_{K P} e^{- Δ G_{y}} + k_{act} e^{Δ G_{act}})}^{2} - 4 {(k_{act})}^{2} e^{Δ G_{act}}} .

(K13)

From Eqs. (K9) and (K10), we may write the stall condition. Noting that, Λ⁺(y, x) is independent of ΔG_pol and Λ⁻(y, x) is proportional to $e^{- Δ G_{pol}}$ ⁠, we may write the stall point as

Γ (N) = - \ln (\frac{Λ_{N}^{+} (r, r)}{e^{Δ G_{pol}} Λ_{N}^{-} (r, r)} + \frac{Λ_{N}^{+} (w, w)}{e^{Δ G_{pol}} Λ_{N}^{-} (w, w)}),

(K14)

such that the dependence on ΔG_pol in the logarithm is canceled out.

REFERENCES

1.

P.

Gaspard

,

Philos. Trans. R. Soc., A

374

,

20160147

(

2016

).

https://doi.org/10.1098/rsta.2016.0147

Google Scholar

Crossref

2.

A. P.

Corfield

,

Biochim. Biophys. Acta, Gen. Subj.

1850

,

236

(

2015

).

https://doi.org/10.1016/j.bbagen.2014.05.003

Google Scholar

Crossref

3.

S.

Pinzón Martín

,

P. H.

Seeberger

, and

D.

Varón Silva

,

Front. Chem.

7

,

710

(

2019

).

https://doi.org/10.3389/fchem.2019.00710

Google Scholar

Crossref

PubMed

4.

M.

Chanda

,

Introduction to Polymer Science and Chemistry: A Problem-Solving Approach

(

CRC Press

,

Boca Raton, FL

,

2013

).

Google Scholar

Crossref

5.

C. G.

Overberger

,

J. Polym. Sci., Polym. Symp.

72

,

67

(

1985

).

https://doi.org/10.1002/polc.5070720110

Google Scholar

Crossref

6.

W.

Meng

,

R. A.

Muscat

,

M. L.

McKee

,

P. J.

Milnes

,

A. H.

El-Sagheer

,

J.

Bath

,

B. G.

Davis

,

T.

Brown

,

R. K.

O’Reilly

, and

A. J.

Turberfield

,

Nat. Chem.

8

,

542

(

2016

).

https://doi.org/10.1038/nchem.2495

Google Scholar

Crossref

PubMed

7.

H.

Zhang

,

Y.

Wang

,

H.

Zhang

,

X.

Liu

,

A.

Lee

,

Q.

Huang

,

F.

Wang

,

J.

Chao

,

H.

Liu

,

J.

Li

et al,

Nat. Commun.

10

(

1

),

1006

(

2019

).

https://doi.org/10.1038/s41467-019-09004-4

Google Scholar

Crossref

PubMed

8.

F.

Crick

,

Nature

227

,

561

(

1970

).

https://doi.org/10.1038/227561a0

Google Scholar

Crossref

PubMed

9.

B.

Alberts

,

A.

Johnson

,

J.

Lewis

,

D.

Morgan

, and

M.

Raff

,

Molecular Biology of the Cell

(

Garland Science, Taylor & Francis Group

,

New York

,

2014

).

Google Scholar

10.

J.

Niu

,

R.

Hili

, and

D. R.

Liu

,

Nat. Chem.

5

,

282

(

2013

).

https://doi.org/10.1038/nchem.1577

Google Scholar

Crossref

PubMed

11.

D.

Kong

,

W.

Yeung

, and

R.

Hili

,

ACS Comb. Sci.

18

,

355

(

2016

).

https://doi.org/10.1021/acscombsci.6b00059

Google Scholar

Crossref

PubMed

12.

A. E.

Stross

,

G.

Iadevaia

,

D.

Núñez-Villanueva

, and

C. A.

Hunter

,

J. Am. Chem. Soc.

139

,

12655

(

2017

).

https://doi.org/10.1021/jacs.7b06619

Google Scholar

Crossref

PubMed

13.

J.-F.

Lutz

,

Sequence-Controlled Polymers

(

Wiley-VCH

,

Weinheim, Germany

,

2018

).

Google Scholar

Crossref

14.

D.

Núñez-Villanueva

,

M.

Ciaccia

,

G.

Iadevaia

,

E.

Sanna

, and

C. A.

Hunter

,

Chem. Sci.

10

,

5258

(

2019

).

https://doi.org/10.1039/c9sc01460h

Google Scholar

Crossref

PubMed

15.

D.

Núñez-Villanueva

and

C. A.

Hunter

,

Acc. Chem. Res.

54

,

1298

(

2021

).

https://doi.org/10.1021/acs.accounts.0c00852

Google Scholar

Crossref

PubMed

16.

J.

Cabello-Garcia

,

W.

Bae

,

G.-B. V.

Stan

, and

T. E.

Ouldridge

,

ACS Nano

15

,

3272

(

2021

).

https://doi.org/10.1021/acsnano.0c10068

Google Scholar

Crossref

PubMed

17.

F. T.

Wall

,

J. Am. Chem. Soc.

63

,

1862

(

1941

).

https://doi.org/10.1021/ja01852a016

Google Scholar

Crossref

18.

F. T.

Wall

,

J. Am. Chem. Soc.

66

,

2050

(

1944

).

https://doi.org/10.1021/ja01240a014

Google Scholar

Crossref

19.

F. R.

Mayo

and

F. M.

Lewis

,

J. Am. Chem. Soc.

66

,

1594

(

1944

).

https://doi.org/10.1021/ja01237a052

Google Scholar

Crossref

20.

T. E.

Ouldridge

,

Nat. Comput.

17

,

3

(

2018

).

https://doi.org/10.1007/s11047-017-9646-x

Google Scholar

Crossref

PubMed

21.

S.

Whitelam

,

R.

Schulman

, and

L.

Hedges

,

Phys. Rev. Lett.

109

,

265506

(

2012

).

https://doi.org/10.1103/physrevlett.109.265506

Google Scholar

Crossref

PubMed

22.

M.

Nguyen

and

S.

Vaikuntanathan

,

Proc. Natl. Acad. Sci. U. S. A.

113

,

14231

(

2016

).

https://doi.org/10.1073/pnas.1609983113

Google Scholar

Crossref

PubMed

23.

J. M.

Poulton

,

P. R.

ten Wolde

, and

T. E.

Ouldridge

,

Proc. Natl. Acad. Sci. U. S. A.

116

,

1946

(

2019

).

https://doi.org/10.1073/pnas.1808775116

Google Scholar

Crossref

PubMed

24.

J. M.

Poulton

and

T. E.

Ouldridge

,

New J. Phys.

23

,

063061

(

2021

).

https://doi.org/10.1088/1367-2630/ac0389

Google Scholar

Crossref

25.

J.

Juritz

,

J. M.

Poulton

, and

T. E.

Ouldridge

,

J. Chem. Phys.

156

,

074103

(

2022

).

https://doi.org/10.1063/5.0077865

Google Scholar

Crossref

PubMed

26.

P.

Sartori

and

S.

Pigolotti

,

Phys. Rev. Lett.

110

,

188101

(

2013

).

https://doi.org/10.1103/physrevlett.110.188101

Google Scholar

Crossref

PubMed

27.

P.

Sartori

and

S.

Pigolotti

,

Phys. Rev. X

5

,

041039

(

2015

).

https://doi.org/10.1103/physrevx.5.041039

Google Scholar

Crossref

28.

M.

Sahoo

,

N.

Arsha

,

P. R.

Baral

, and

S.

Klumpp

,

Phys. Rev. E

104

,

034417

(

2021

).

https://doi.org/10.1103/physreve.104.034417

Google Scholar

Crossref

PubMed

29.

Y.-S.

Song

,

Y.-G.

Shu

,

X.

Zhou

,

Z.-C.

Ou-Yang

, and

M.

Li

,

J. Phys.: Condens. Matter

29

,

025101

(

2016

).

https://doi.org/10.1088/0953-8984/29/2/025101

Google Scholar

Crossref

PubMed

30.

Y.

Song

and

C.

Hyeon

,

J. Phys. Chem. Lett.

11

,

3136

(

2020

).

https://doi.org/10.1021/acs.jpclett.0c00545

Google Scholar

Crossref

PubMed

31.

Q.-S.

Li

,

P.-D.

Zheng

,

Y.-G.

Shu

,

Z.-C.

Ou-Yang

, and

M.

Li

,

Phys. Rev. E

100

,

012131

(

2019

).

https://doi.org/10.1103/physreve.100.012131

Google Scholar

Crossref

PubMed

32.

S.

Pigolotti

and

P.

Sartori

,

J. Stat. Phys.

162

,

1167

(

2016

).

https://doi.org/10.1007/s10955-015-1399-2

Google Scholar

Crossref

33.

F.

Wong

,

A.

Amir

, and

J.

Gunawardena

,

Phys. Rev. E

98

,

012420

(

2018

).

https://doi.org/10.1103/physreve.98.012420

Google Scholar

Crossref

PubMed

34.

W. D.

Piñeros

and

T.

Tlusty

,

Phys. Rev. E

101

,

022415

(

2020

).

https://doi.org/10.1103/physreve.101.022415

Google Scholar

Crossref

PubMed

35.

C. H.

Bennett

,

BioSystems

11

,

85

(

1979

).

https://doi.org/10.1016/0303-2647(79)90003-0

Google Scholar

Crossref

PubMed

36.

P.

Gaspard

,

Phys. Rev. E

93

,

042420

(

2016

).

https://doi.org/10.1103/PhysRevE.93.042420

Google Scholar

Crossref

PubMed

37.

P.

Gaspard

,

Phys. Rev. E

93

,

042419

(

2016

).

https://doi.org/10.1103/physreve.93.042419

Google Scholar

Crossref

PubMed

38.

D.

Andrieux

and

P.

Gaspard

,

J. Chem. Phys.

130

,

014901

(

2009

).

https://doi.org/10.1063/1.3050099

Google Scholar

Crossref

PubMed

39.

R.

Rao

and

L.

Peliti

,

J. Stat. Mech.: Theory Exp.

2015

,

P06001

.

https://doi.org/10.1088/1742-5468/2015/06/p06001

Crossref

40.

K.

Banerjee

,

A. B.

Kolomeisky

, and

O. A.

Igoshin

,

Proc. Natl. Acad. Sci. U. S. A.

114

,

5183

(

2017

).

https://doi.org/10.1073/pnas.1614838114

Google Scholar

Crossref

PubMed

41.

D.

Chiuchiù

,

Y.

Tu

, and

S.

Pigolotti

,

Phys. Rev. Lett.

123

,

038101

(

2019

).

https://doi.org/10.1103/physrevlett.123.038101

Google Scholar

Crossref

PubMed

42.

T. E.

Ouldridge

and

P.

Rein ten Wolde

,

Phys. Rev. Lett.

118

,

158103

(

2017

).

https://doi.org/10.1103/physrevlett.118.158103

Google Scholar

Crossref

PubMed

43.

J. J.

Hopfield

,

Proc. Natl. Acad. Sci. U. S. A.

71

,

4135

(

1974

).

https://doi.org/10.1073/pnas.71.10.4135

Google Scholar

Crossref

PubMed

44.

J.

Ninio

,

Biochimie

57

,

587

(

1975

).

https://doi.org/10.1016/s0300-9084(75)80139-8

Google Scholar

Crossref

PubMed

45.

J. D.

Mallory

,

O. A.

Igoshin

, and

A. B.

Kolomeisky

,

J. Phys. Chem. B

124

,

9289

(

2020

).

https://doi.org/10.1021/acs.jpcb.0c06180

Google Scholar

Crossref

PubMed

46.

A.

Murugan

,

D. A.

Huse

, and

S.

Leibler

,

Phys. Rev. X

4

,

021016

(

2014

).

https://doi.org/10.1103/physrevx.4.021016

Google Scholar

Crossref

47.

A.

Murugan

,

D. A.

Huse

, and

S.

Leibler

,

Proc. Natl. Acad. Sci. U. S. A.

109

,

12034

(

2012

).

https://doi.org/10.1073/pnas.1119911109

Google Scholar

Crossref

PubMed

48.

Q.

Yu

,

A. B.

Kolomeisky

, and

O. A.

Igoshin

,

J. R. Soc., Interface

19

,

20210883

(

2022

).

https://doi.org/10.1098/rsif.2021.0883

Google Scholar

Crossref

PubMed

49.

M.

Sahoo

and

S.

Klumpp

,

J. Phys.: Condens. Matter

25

,

374104

(

2013

).

https://doi.org/10.1088/0953-8984/25/37/374104

Google Scholar

Crossref

PubMed

50.

M.

Ehrenberg

and

C.

Blomberg

,

Biophys. J.

31

,

333

(

1980

).

https://doi.org/10.1016/s0006-3495(80)85063-6

Google Scholar

Crossref

PubMed

51.

V.

Galstyan

and

R.

Phillips

,

J. Phys. Chem. B

123

,

10990

(

2019

).

https://doi.org/10.1021/acs.jpcb.9b08380

Google Scholar

Crossref

PubMed

52.

P.

Gaspard

and

D.

Andrieux

,

J. Chem. Phys.

141

,

044908

(

2014

).

https://doi.org/10.1063/1.4890821

Google Scholar

Crossref

PubMed

53.

T. L.

Hill

,

J. Theor. Biol.

10

,

442

(

1966

).

https://doi.org/10.1016/0022-5193(66)90137-8

Google Scholar

Crossref

PubMed

54.

T. L.

Hill

,

Proc. Natl. Acad. Sci. U. S. A.

85

,

2879

(

1988

).

https://doi.org/10.1073/pnas.85.9.2879

Google Scholar

Crossref

PubMed

55.

J. G.

Kemeny

and

J. L.

Snell

,

Finite Markov Chains

(

Springer

,

New York

,

1983

).

Google Scholar

56.

V.

Anantharam

and

P.

Tsoucas

,

Stat. Probab. Lett.

8

,

189

(

1989

).

https://doi.org/10.1016/0167-7152(89)90016-3

Google Scholar

Crossref

57.

A.

Wachtel

,

R.

Rao

, and

M.

Esposito

,

New J. Phys.

20

,

042002

(

2018

).

https://doi.org/10.1088/1367-2630/aab5c9

Google Scholar

Crossref

58.

H.-H.

Kohler

and

E.

Vollmerhaus

,

J. Math. Biol.

9

,

275

(

1980

).

https://doi.org/10.1007/bf00276029

Google Scholar

Crossref

PubMed

59.

F.

Cady

and

H.

Qian

,

Phys. Biol.

6

,

036011

(

2009

).

https://doi.org/10.1088/1478-3975/6/3/036011

Google Scholar

Crossref

PubMed

60.

P.

Gaspard

,

Phys. Rev. Lett.

117

,

238101

(

2016

).

https://doi.org/10.1103/physrevlett.117.238101

Google Scholar

Crossref

PubMed

61.

M.

Esposito

,

Phys. Rev. E

85

,

041125

(

2012

).

https://doi.org/10.1103/physreve.85.041125

Google Scholar

Crossref

62.

T. E.

Ouldridge

,

R. A.

Brittain

, and

P. R.

ten Wolde

, in

The Energetics of Computing in Life and Machines

, edited by

C.

Kempes

,

D. H.

Wolpert

,

P. F.

Stadler

, and

J. A.

Grochow

(

SFI Press

,

Santa Fe, NM

,

2019

), pp.

307

–

351

.

Google Scholar

Crossref

63.

T. M.

Cover

and

J. A.

Thomas

,

Elements of Information Theory

(

John Wiley & Sons, Inc.

,

Hoboken, NJ

,

2006

).

Google Scholar

64.

M.

Esposito

,

K.

Lindenberg

, and

C.

Van den Broeck

,

J. Stat. Mech.: Theory Exp.

2010

,

P01008

.

https://doi.org/10.1088/1742-5468/2010/01/p01008

Crossref

65.

D. T.

Gillespie

,

J. Comput. Phys.

22

,

403

(

1976

).

https://doi.org/10.1016/0021-9991(76)90041-3

Google Scholar

Crossref

66.

D.

Chiuchiù

,

J.

Ferrare

, and

S.

Pigolotti

,

Phys. Rev. E

100

,

062502

(

2019

).

https://doi.org/10.1103/physreve.100.062502

Google Scholar

Crossref

PubMed

67.

M.

Das

and

H.

Kantz

,

Phys. Rev. E

103

,

032110

(

2021

).

https://doi.org/10.1103/physreve.103.032110

Google Scholar

Crossref

PubMed

68.

V.

Galstyan

,

K.

Husain

,

F.

Xiao

,

A.

Murugan

, and

R.

Phillips

,

eLife

9

,

e60415

(

2020

).

https://doi.org/10.7554/elife.60415

Google Scholar

Crossref

PubMed

69.

P.

Gaspard

,

Phys. Rev. E

96

,

042403

(

2017

).

https://doi.org/10.1103/physreve.96.042403

Google Scholar

Crossref

PubMed

70.

D. A.

Harville

,

Matrix Algebra from a Statistician’s Perspective

(

Springer

,

New York

,

1998

).

Google Scholar

2023

Author(s)

A universal method for analyzing copolymer growth

I. INTRODUCTION

II. METHODS

A. Absorbing Markov chains

1. Expectations of an absorbing process are steady state averages of a “closed process”

2. Steady state averages of the closed process are calculated using the Markov chain tree theorem

3. Absorbing probabilities

4. Counting edge and cycle transitions

B. Copolymer methods

1. Philosophy of coarse-graining complex underlying copolymerization reaction networks

2. Identifying propensities in the coarse-grained model

3. Solving the coarse-grained model

4. Extracting properties of the fine-grained model from the solution of the coarse-grained model

5. Stalled growth

6. Limiting behavior

7. Simplification for factorizable propensities

III. EXAMPLE APPLICATIONS

A. Stalling behavior in a polymerization model with no neighbor–neighbor interactions

B. Balanced models of templated polymerization with autonomous separation

C. Hopfield’s kinetic proofreading in a model of templated copying with autonomous separation

IV. CONCLUSION

SUPPLEMENTARY MATERIAL

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: FACTORIZING SUMS OF SPANNING TREES

APPENDIX B: NORMALIZATION CONSTANT FOR EXAMPLE ABSORBING MARKOV PROCESS

APPENDIX C: EQUIVALENCE BETWEEN CHEMICAL WORK CALCULATED FROM EDGES AND CYCLES

APPENDIX D: CYCLES OF THE EXAMPLE ABSORBING PROCESS

APPENDIX E: NUMBER OF STEPS PER NET FORWARD STEP OF A RANDOM WALK

APPENDIX F: THE FREQUENCY AT STALL IS GIVEN BY THE DIAGONAL COFACTORS OF A MATRIX

APPENDIX G: THE FREQUENCIES IN THE IRREVERSIBLE LIMIT ARE GIVEN BY THE STEADY STATE OF A PROCESS OF THE COMPLETE GRAPH

APPENDIX H: SIMPLIFICATION OF RESULTS FOR FACTORIZABLE RATIOS OF PROPENSITIES

APPENDIX I: FREQUENCY FOR A BALANCED MODEL WITH TWO MONOMER TYPES IN THE SLOW POLYMERIZATION LIMIT

APPENDIX J: MODEL USED FOR BALANCED ON-RATE VS OFF-RATE DISCRIMINATION COMPARISONS

APPENDIX K: EQUATIONS FOR N-LOOP HOPFIELD MODEL

REFERENCES

Supplementary Material

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only