Assembly of heteropolymers via a network of reaction coordinates

TL;DR

This paper introduces a novel theoretical framework for modeling heteropolymer assembly using a continuous reaction coordinate, enabling estimation of error rates and speeds from the Gibbs free energy landscape, including complex error correction mechanisms.

Contribution

It presents a new approach applying reaction rate theory to heteropolymer assembly, capturing non-equilibrium dynamics and error correction in a unified continuous model.

Findings

The model estimates incorporation speed and error rates from free energy landscapes.

Application to simple and complex systems demonstrates versatility.

Incorporates error correction mechanisms like kinetic proofreading.

Abstract

In biochemistry, heteropolymers encoding biological information are assembled out of equilibrium by sequentially incorporating available monomers found in the environment. Current models of polymerization treat monomer incorporation as a sequence of discrete chemical reactions between intermediate meta-stable states. In this paper, we use ideas from reaction rate theory and describe non-equilibrium assembly of a heteropolymer via a continuous reaction coordinate. Our approach allows to estimate the copy error and incorporation speed from the Gibbs free energy landscape of the process. We apply our theory to several examples, from a simple reaction characterized by a free energy barrier to more complex cases incorporating error correction mechanisms such as kinetic proofreading.

Equations

$$\text{schemestart}h\text{arrow}(@c1--)<=>[+w][-w][-30]hw\text{arrow}(@c1--)<=>[+r][-r][30]hr\text{schemestop}$$

$$\Delta G^x(y)=G^{hx}(y)-G^{hx}(0),$$

$$G^{hx}(y)=$$

$$=$$

$$\frac{\mathrm{d}y}{\mathrm{d}t}=-\mu \frac{\mathrm{d}}{\mathrm{d}y}{G}^{hx}(y)+\sqrt{2D}\xi (t)$$

$$(y=0)\underset{{\pi }_{ϵ,0}^x}{\overset{{\pi }_{0,ϵ}^x}{\rightleftharpoons }}(y=ϵ)\underset{{\pi }_{1-ϵ,ϵ}^x}{\overset{{\pi }_{ϵ,1-ϵ}^x}{\rightleftharpoons }}(y=1-ϵ)\underset{{\pi }_{1,1-ϵ}^x}{\overset{{\pi }_{1-ϵ,1}^x}{\rightleftharpoons }}(y=1)$$

$${\pi }_{ϵ,0}^x={\pi }_{1-ϵ,1}^x=\frac{1}{3}+\mathcal{O}(ϵ).$$

$${\pi }_{0,ϵ}^x$$

$${\pi }_{ϵ,1-ϵ}^x$$

$${\pi }_{1-ϵ,ϵ}^x$$

$${\pi }_{1,1-ϵ}^x$$

$${\psi }^x(y)=$$

$$=$$

$$p_{\to }^x=$$

$$p_{\leftarrow }^x=$$

$$\frac{p_{\to }^x}{p_{\leftarrow }^x}=\exp \left[-\frac{\Delta G^x(1)}{k_{\mathrm{B}}T}\right].$$

$$\eta =\underset{N\to \infty }{\lim }\frac{N^w}{N}.$$

$$\frac{\eta }{1-\eta }=\underset{N\to \infty }{\lim }\frac{N^w}{N^r}=\underset{n\to \infty }{\lim }\frac{n_{\to }^w-n_{\leftarrow }^w}{n_{\to }^r-n_{\leftarrow }^r}$$

$$n_{\to }^r\sim$$

$$n_{\leftarrow }^r\sim$$

$$n_{\to }^w\sim$$

$$n_{\leftarrow }^w\sim$$

$$\frac{\eta }{1-\eta }=\frac{p_{\to }^w-\eta p_{\leftarrow }^w}{p_{\to }^r-(1-\eta )p_{\leftarrow }^r}$$

$$\frac{\eta }{1-\eta }=$$

$$\eta =\frac{1}{1+\exp \left[\frac{\Delta G^w(1)-\Delta G^r(1)}{k_{\mathrm{B}}T}\right]}\text{mbox}if\exp \left[\frac{\Delta G^w(1)-\Delta G^r(1)}{k_{\mathrm{B}}T}\right]\frac{{\int }_0^1\exp \left[\frac{\Delta G^r(y)}{k_{\mathrm{B}}T}\right]\mathrm{d}y}{{\int }_0^1\exp \left[\frac{\Delta G^w(y)}{k_{\mathrm{B}}T}\right]\mathrm{d}y}=1,$$

$$\frac{{\int }_0^1\exp \left[\frac{\Delta G^r(y)}{k_{\mathrm{B}}T}\right]\mathrm{d}y}{{\int }_0^1\exp \left[\frac{\Delta G^w(y)}{k_{\mathrm{B}}T}\right]\mathrm{d}y}\approx \exp \left[\frac{{\delta }^r-{\delta }^w}{k_{\mathrm{B}}T}\right]\sqrt{\frac{{\Sigma }^w}{{\Sigma }^r}}$$

$$v=$$

$$=$$

$$\mathcal{T}\sim n⟨\tau ⟩$$

$$v=\frac{c[p_{\to }^r-(1-\eta )p_{\leftarrow }^r+p_{\to }^w-\eta p_{\leftarrow }^w]}{⟨\tau ⟩}.$$

$$\Delta G^r(y)=-m_{r}y$$

$$\Delta G^w(y)=-m_{w}y.$$

$$\frac{\eta }{1-\eta }=\frac{m_w(1-e^{-m_r})[1-e^{-m_w}\eta ]}{m_r(1-e^{-m_w})[1+e^{-m_r}(1-\eta )]},$$

$$\eta =\frac{m_w}{m_r+m_w}\text{mbox}for{m}_r,m_w\gg 1.$$

$$\Delta G^x(y)=$$

Full text

I Introduction

DNA, RNA, and proteins are the building blocks of all living systems. These heteropolymers are assembled to match a template; only a very small number of mismatches with the template is tolerable for maintaining biological information and for correct functioning of cells. However, the binding energies of different monomers usually differ by only a few $k_{B}T$, where $k_{B}$ is the Boltzmann constant and $T$ the temperature. This means that, at physiological temperature, mismatches can not be completely suppressed Reynolds et al. (2010).

Our aim is to describe the chemical processes responsible for these errors. Specifically, we consider sequential assembly of heteropolymers where each incorporated monomer can be a right ($r$) or a wrong ($w$) match with a template. These two different outcomes can be represented as competing chemical reactions

[TABLE]

where $h$ is the heteropolymer produced so far, and $hr$/$hw$ are the same heteropolymer with an addition of a $r$/$w$ monomer at the tip, respectively. Each monomer incorporation is iteratively followed by a new one, so that the whole polymerization process is described by the tree-shaped network of chemical reactions Bennett (1979); Pigolotti and Sartori (2016) in Fig. 1a.

To achieve accurate and fast assembly, the reactions in Eq.(1) involve several intermediate steps, such as initial monomer discrimination Kunkel and Bebenek (2000), kinetic proofreading, Hopfield (1974); Ninio (1975); Kunkel and Bebenek (2000), and mismatch repair O’Donnell et al. (2013); Kunkel and Erie (2015). In general, each one of these error-correction mechanisms contribute simultaneously to polymerization accuracy, speed Rodnina and Wintermeyer (2001); Andrieux and Gaspard (2008); Johansson et al. (2008); Gaspard and Andrieux (2014); Banerjee et al. (2017); Sartori and Pigolotti (2015); Savir and Tlusty (2013); Sartori and Pigolotti (2013); Murugan et al. (2012), and energetic cost Cady and Qian (2009); Kramers (1940); Rao and Peliti (2015); Sartori and Pigolotti (2015); Wagoner and Dill (2019).

Two approaches can provide insight into the error-correction mechanisms underlying heteropolymer assembly. The first approach is to measure their kinetic rates under different experimental conditions Rodnina and Wintermeyer (2001). The second approach is to simulate heteropolymer assembly using molecular dynamics Bock et al. (2018). From the molecular dynamics, one can project the numerous degrees of freedom into a 1-dimensional collective variable called reaction coordinate Banushkina and Krivov (2016). The reaction coordinate simplifies a chemical process into a one-dimensional random motion Banushkina and Krivov (2016); Socci et al. (1996); Best and Hummer (2006). The parameters of this random motion depend on the underlying reactants dynamics Zwanzig (1961); Berezhkovskii and Szabo (2013); Lu and Vanden-Eijnden (2014) and on the projection technique Lu and Vanden-Eijnden (2014); Banushkina and Krivov (2016); Krivov and Karplus (2006).

While successful in describing protein folding Socci et al. (1996); Best and Hummer (2006); Klimov and Thirumalai (1997) and in modeling reaction rates Kramers (1940); Klimov and Thirumalai (1997), approaches based on reaction coordinates found little use in studies of polymerization speed and accuracy. In principle, both reactions in Eq. (1) can be described by means of a reaction coordinate (Fig. 1b). However, to study the complete polymerization process we need to join the reaction coordinates characterizing each branch in Fig. 1a. Mathematically, this amounts to impose appropriate boundary conditions at the nodes of the reaction network.

In this paper, we develop a model of heteropolymer assembly based on reaction coordinates, and use it to compute the accuracy and speed of polymerization in different conditions. The paper is organized as follows. In Section II, we introduce our model. From the reaction coordinate, we derive effective incorporation and removal probabilities of right and wrong monomers. In Section III, we compute the accuracy and speed of a general heteropolymer assembly. In Section IV we consider examples characterized by different Gibbs free energy landscapes. In Section V, we generalize our results to a case where the reaction leading to monomer incorporation is complemented by kinetic proofreading. Section VI is devoted to conclusions and perspectives.

II Model

We define our model of heteropolymer assembly with reaction coordinates through the following steps. We first introduce the reaction coordinate and the free energy landscape in each chemical reaction of the polymerization network. We then study the dynamics of the reaction coordinate dynamics and its boundary conditions at the nodes of the network. Finally, we compute the probabilities to incorporate/remove one monomer along each reaction coordinate.

II.1 Reaction coordinate and Gibbs free energy of the heteropolymer

We introduce the continuous reaction coordinate $y$ along each edge of the polymerization network, Fig. 1a. Without loss of generality, we choose the units of the reaction coordinate so that $y\in[0,1]$, where $y=0$ and $y=1$ correspond to $h$ and $hx$ respectively, i.e. to the states before and after monomer incorporation, see Figure 1.b.

Each point along this reaction coordinate is characterized by a Gibbs free energy $G^{hx}(y)$ (from now on simply ”free energy”). Such free energy depends on the previously incorporated sequence of monomers ($h$), on the candidate monomer to be incorporated ($x$) and on the stage of the incorporation process, i.e. the value of $y$. Implicitly, $G^{hx}(y)$ also depends on the reactant and product concentrations.

We introduce the free energy increments from the beginning of each incorporation reaction

[TABLE]

see Fig. 2. The free energy increments depend on the candidate monomer $x$ but not on the whole history of incorporated monomers $h$. With this notation, the (absolute) binding free energy of monomer $x$ is equal to $-\Delta G^{x}(1)$.

The free energy must be a continuous function of $y$, and must also vary continuously when crossing the nodes of the network in Fig. 1. This means that we can decompose the free energy at an arbitrary stage of the polymerization process as

[TABLE]

II.2 Stochastic dynamics of the reaction coordinate and boundary conditions

Because of thermal fluctuations, the reaction coordinate $y$ evolves according to a Langevin equation

[TABLE]

where $\mu$ is a mobility, $D$ is a diffusion coefficient, and $\xi(t)$ is white noise with $\langle\xi(t)\rangle=0$ and $\langle\xi(t)\xi(t^{\prime})\rangle=\delta(t-t^{\prime})$ Gardiner (2009). We assume that $D$ satisfies the Einstein relation $D=k_{B}T\mu$ with temperature $T$ and Boltzmann constant $k_{B}$. We also assume that $D$, $\mu$, and $T$ are constant. When the reaction coordinate reaches the boundaries, either $y=0$ or $y=1$, a new incorporation/removal reaction is commenced.

Equation (4) needs to be complemented by rules to specify which reaction initiates at the nodes of the reaction network. To this aim, we consider two intermediate values of the reaction coordinate: $y=\epsilon$ and $y=1-\epsilon$ with $\epsilon\ll 1$. Using these values we coarse-grain the evolution of the reaction coordinate $y$ as

[TABLE]

where the quantities $\pi_{\tilde{y},y}^{x}$ are first-passage probability from $y$ to $\tilde{y}$. For example, $\pi_{1-\epsilon,\epsilon}^{x}$ is the probability that the reaction coordinate reaches $y=1-\epsilon$ from $y=\epsilon$ without having reached $y=0$ before.

The representation in Eq.(5) separates the dynamics in proximity of the nodes of Fig. 1.b, from the dynamics in the interval $[\epsilon,1-\epsilon]$. Thanks to this separation, we use detailed balance, probability conservation close to the nodes, and the continuity of $G^{hx}(y)$ to compute the first-passage probabilities $\pi_{\epsilon,0}^{x}$ and $\pi_{1-\epsilon,1}^{x}$ (see Appendix A). This procedure results in

[TABLE]

We compute the first-passage probabilities in the interior applying standard techniques Gardiner (2009); Iyer-Biswas and Zilman (2016) to the Fokker-Planck equation associated to Eq. (4). We obtain

[TABLE]

with

[TABLE]

where the last equality follows from the relation $D=k_{B}T\mu$.

II.3 Effective probabilities of monomer incorporation/rejection

From the probabilities $\pi_{\tilde{y},y}^{x}$, we now compute the effective probabilities $p_{\rightarrow}^{x}$ and $p_{\leftarrow}^{x}$ to incorporate and reject monomer $x$ along each edge of the reaction network in Figure 1. To this end, we assume that the coarse grained dynamics in Eq.(5) is at steady state. We then use adiabatic elimination Pigolotti and Vulpiani (2008) to obtain (see Appendix B)

[TABLE]

Equations (9) relate the free energy landscapes $G^{x}(y)$ and the incorporation/removal probabilities of the polymerization process. They are consistent with the detailed balance condition

[TABLE]

which connects the ratios of forward and backward probabilities to the binding free energy $-\Delta G^{x}(1)$, see Fig. 2.

III Results

We now address the accuracy and speed of a polymerization process in the reaction coordinate framework. We consider a copy polymer made up of a number $N^{r}$ of right monomers and $N^{w}$ of wrong monomers with $N=N^{r}+N^{w}$. For large $N$, we define the error rate

[TABLE]

To compute $\eta$ from the incorporation and removal probabilities $p_{\rightarrow}^{x}$ and $p_{\leftarrow}^{x}$, we first recast Eq.(11) into the implicit equation

[TABLE]

where we have introduced the numbers $n_{\rightarrow}^{r}$, $n_{\leftarrow}^{r}$, $n_{\rightarrow}^{w}$ and $n_{\leftarrow}^{w}$ of $r$ and $w$ incorporation and removal reactions which have occurred in the process, and $n$ is the total number of observed chemical reactions. For large $n$ we have

[TABLE]

where $c=[p_{\rightarrow}^{r}+(1-\eta)p_{\leftarrow}^{r}+p_{\rightarrow}^{w}+\eta p_{\leftarrow}^{w}]^{-1}$ is a normalization constant so that $n=n_{\rightarrow}^{r}+n_{\leftarrow}^{r}+n_{\rightarrow}^{w}+n_{\leftarrow}^{w}=n$. Substituting Eqs. (13) into Eq. (12) gives

[TABLE]

Equation (14) is a general ”self-consistency” relation for the error rate that holds also for discrete models of polymerization Sartori and Pigolotti (2015); Pigolotti and Sartori (2016); Bennett (1979). In our case, we substitute Eqs. (9) in Eq. (14) and take the limit $\epsilon\to 0$, obtaining

[TABLE]

Solving Eq. 15 for $\eta$ yields an explicit expression of the error rate from the energy potentials.

Equation 15 identifies different regimes of error correction. To identify a first regime, we observe that

[TABLE]

In the regime where Eq. (16) holds, the error depends only on the binding free energy difference $\Delta G^{w}(1)-\Delta G^{r}(1)$. This regime is called energetic discrimination regime in the literature Sartori and Pigolotti (2013); Pigolotti and Sartori (2016). Systems near equilibrium operates in this regime because the Boltzmann factors of the binding free energies determine, via detailed balance, the probabilities to incorporate different monomers.

To identify a second error-correction regime in Eq. (15), we consider the case where $\Delta G^{r}(y)$ and $\Delta G^{w}(y)$ are characterized by energy barriers with heights $\delta^{r}$ and $\delta^{w}$ respectively (see Figure 1.b and Kramers Kramers (1940)). When such barriers are large, we can approximate the integrals in Eq. (15) by using the Laplace method Bender and Orszag (1978)

[TABLE]

where $\Sigma^{r}$ and $\Sigma^{w}$ are the curvatures of $\Delta G^{r}(y)$ and $\Delta G^{w}(y)$ at their maxima, respectively. Equation 17 implies that activation barriers suppress the polymerization error via the second term in round brackets in Eq. (15). The regime where this suppression occurs is the kinetic discrimination regime Sartori and Pigolotti (2013); Pigolotti and Sartori (2016). The first factor on the right-hand side of Eq. (17) represents the contribution of a difference in activation energy barrier $\delta^{r}-\delta^{w}$ to accuracy. This effect is also present in models based on discrete-step reactions Bennett (1979); Cady and Qian (2009); Sartori and Pigolotti (2013); Pigolotti and Sartori (2016). The factor $\sqrt{\Sigma^{w}/\Sigma^{r}}$ is a correction to activation energies based on the width of the activation barriers. This factor permits kinetic discrimination at equal barrier heights, provided that the barrier for right monomers is significantly more narrow than for wrong monomers.

We estimate the average polymerization speed using a similar argument to that leading to Eq. (15). For large number $N$ of incorporated monomers, the average speed is equal to $N$ divided the total time $\mathcal{T}$ needed to assemble the polymer

[TABLE]

where we expressed $N$ in terms of the number of incorporation/removal reactions. For large $N$ we can approximate the polymerization time as

[TABLE]

where $\langle\tau\rangle$ is the average time it takes to either incorporate or remove a monomer. Substituting Eqs. (13) and (19) into Eq. (18) gives the estimate for the polymerization speed

[TABLE]

The numerator of Eq. (18) is the probability of an incorporation minus the probability of a removal, while the denominator provides the timescale of these events. In practice, calculating $\langle\tau\rangle$ is not straightforward since one has to take into account contributions from incorporation attempts that are not finalized. In Appendix C, we provide a more formal derivation of Eq. (20), together with an explicit expression for $\langle\tau\rangle$.

IV Examples

To address the validity and practical implications of Eqs. (15) and (20) we consider two examples of potentials $\Delta G^{r}(y)$ and $\Delta G^{w}(y)$. In both cases, we work in dimensionless units by fixing $D=1$, $k_{\mathrm{B}}T=1$, and $\mu=1$.

IV.1 Linear potential

As first example we consider linear free energy landscapes

[TABLE]

Despite their simplicity, the potentials in Eq. (21) are useful to understand the physics of the process. Upon increasing the slopes $m_{r}$ and $m_{w}$, polymerization becomes increasingly irreversible. Substituting the potentials Eq. (21) into the expression for the error, Eq. (15) and performing the integrals we obtain

[TABLE]

which implies

[TABLE]

The exact solution of Eq. (22) shows that the error is approximately a function of $m_{w}/m_{r}$ when $m_{r}$, $m_{w}$ are large, as predicted by Eq. (23), Fig. 3a. We compared the predictions from Eqs. (22) and (20) with numerical simulations of the incorporation process from Eq.(4). Our theory yields reliable predictions for a broad range of parameters, Fig. 3c and 3d.

IV.2 Potential with an activation barrier

As a second example we consider the potential

[TABLE]

where $a_{x}$, $b_{x}$ and $c_{x}$ are monomer-dependent parameters that control the shape of the free energy potentials. Key features of the potential of Eq. (24) are the binding energy $-\Delta G^{x}(1)$, the height of the activation barrier $\delta^{x}$ and its width $\sigma^{x}=4c_{x}$, Fig. 4.

We study this model for different cases, corresponding to different parameter choices. In the first case we fix $-\Delta G^{r}(1)=-\Delta G^{w}(1)$ upon choosing $b_{r}=b_{w}=b$ and $c_{r}=c_{w}=1/20$. This enforces a kinetic discrimination regime Pigolotti and Sartori (2016) where the binding energy $-\Delta G^{r}(1)$ quantifies the degree of irreversibility. For highly irreversible processes, the error $\eta$ should mainly depend on the activation energy difference $\delta^{r}-\delta^{r}$, see Eq. (17). We also expect that the reaction speed should increase for more irreversible processes. Equations (15) and (20) confirm such qualitative picture, see Figure 5a and b. Also in this case, numerical simulations are in excellent quantitative agreement with our theory, Fig. 5c and 5d.

As a second case, we fix $a_{r}=a_{w}=5$ and $b_{r}=b_{w}=1$. In this way we have that $-\Delta G^{r}(1)\approx-\Delta G^{w}(1)$ and $\delta^{r}\approx\delta^{w}$. Energetics alone would not permit monomer discrimination in this case Pigolotti and Sartori (2016). However, Eq. (17) predicts that the difference in the barrier widths $\sigma_{r}$ and $\sigma_{w}$ should allow to discriminate $r$ and $w$ monomers (see Figure 6.a). We confirmed the existence of such kinetic discrimination regime with numerical simulations, Fig. 6.c.

V Kinetic proofreading

In this Section we sketch a generalization of our framework to include kinetic proofreading Hopfield (1974); Ninio (1975). We assume that the reaction $h\rightleftharpoons hx$ can be decomposed into three sub-reactions

[TABLE]

where each sub-reaction occurs with probabilities $p_{\rightarrow}^{i,x}$s and $p_{\leftarrow}^{i,x}$s, and $hx^{*}$ is an intermediate meta-stable complex. The extra pathway $hx^{*}\rightleftharpoons h$ represents kinetic proofreading. Such reaction can improve accuracy when driven towards the reactants $h$, so that wrong monomers undergo an additional checkpoint. Hopfield (1974); Pigolotti and Sartori (2016).

Every sub-reactions in Eq. (25) is described by its own reaction coordinate $y$ which evolves according to a Langevin equation

[TABLE]

where $G^{i,hx}(y)$ is the free energy landscapes along the $i$-th sub-reaction. Also in this case we take $y\in[0,1]$ for all sub-reactions, with $y=1$ always in the direction of incorporation of monomer $x$. Similarly to Eq. (3), we decompose the free energies for the sub-reactions as

[TABLE]

where we specified that monomer $\tilde{x}$ was incorporated before attempting to incorporate monomer $x$. Here, $G^{3,h\tilde{x}x}(y)$ depend on the direction of the sub-reaction because the heteropolymer total energy now depends also on the sequence of sub-reactions.

We now compute the probabilities $p_{\rightarrow}^{i,x}$s and $p_{\leftarrow}^{i,x}$ with $i\in\{1,2,3\}$ from Eq. (26) with the same procedure which leads to Eq.(9). This yields

[TABLE]

with

[TABLE]

Equations (29) state that the all sub-reactions from reactants $h$ and $hx^{*}$ respectively can start with equal probabilities.

To obtain an equation for $\eta$, we need to compute the effective incorporation and removal probabilities $p_{\rightarrow}^{x}$ and $p_{\leftarrow}^{x}$ in Eq. (14) from Eqs. (28) and (29). To this end we assume that the reactions in Eq.(25) are at steady state. We then use adiabatic elimination Pigolotti and Vulpiani (2008) to obtain (see Appendix D)

[TABLE]

Substituting Eqs. (28) and (30) in Eq. (14) finally provides an expression for $\eta$ in terms of the free energy landscapes $G^{i,x}(y)$.

VI Conclusions

In this paper, we described assembly of heteropolymers by means of continuous reaction coordinates. In the simplest cases, our results are consistent with those derived for reactions occurring in discrete steps Andrieux and Gaspard (2008); Gaspard and Andrieux (2014); Banerjee et al. (2017); Sartori and Pigolotti (2015, 2013); Murugan et al. (2012); Pigolotti and Sartori (2016); Bennett (1979). Moreover, our formalism reveals discrimination mechanisms that are not easily described with discrete reactions. One example is the possibility to discriminate according to barrier widths, as described by Eq. (17) and confirmed in simulations, Fig. 6c.

For simplicity, in this paper we developed our formalism by means of a reaction coordinate characterized by a Markovian dynamic. In general, only specific projection techniques yield reaction coordinates with negligible non-Markovian contributions Banushkina and Krivov (2016); Lu and Vanden-Eijnden (2014); Krivov and Karplus (2006); Krivov (2018), and the resulting Langevin equation might not be in the form of Eq. (4). Our framework can be adapted to such situations as well as to non-Markovian reaction coordinates, describing for example enzymes undergoing slow conformational changes.

The framework described here is microscopically reversible. This allows to characterize non-equilibrium work and heat exchanges during the polymerization process from the diffusive dynamics of the reaction coordinate, similarly to recent studies of the ATP synthase Lucero et al. (2019); Kasper and Sivak (2019) and small-scale technological devices Neri et al. (2015); López-Suárez et al. (2016). This analysis would permit to characterize thermodynamic limits of information processing of these processes Bérut et al. (2012); Jun et al. (2014); Chiuchiú (2015); Sartori and Pigolotti (2015).

Acknowledgements.

This work was supported by JSPS KAKENHI Grant Number JP18K03473 (to DC and SP).

Appendix A First passage time probabilities at the nodes

Because of detailed balance, the probabilities $\pi_{\epsilon,0}^{x}$ and $\pi_{1-\epsilon,1}^{x}$ are related to the free energy difference when passing from one edge of the reaction network to another, i.e.

[TABLE]

where we specified that the monomer $\tilde{x}\in\{r,w\}$ was incorporated before monomer $x\in\{r,w\}$. After the incorporation of $\tilde{x}$, the enzyme can catalyze three reactions: removal of $\tilde{x}$ or incorporation of either $r$ or $w$. The probabilities of these three events must be normalized

[TABLE]

Combining Eqs.(A)-(32) gives

[TABLE]

Substituting Eq.(3) into Eq.(33), taking the limit of small $\epsilon$, using the continuity of $G^{h\tilde{x}x}(y)$ and then renaming $\tilde{x}$ with $x$ finally gives Eq.(6).

Appendix B Effective incorporation and removal probabilities

The dynamics of the reaction coordinate $y$ in the coarse grained description of Eq.(5) obeys a Markov chain

[TABLE]

where the first-passage probabilities appear as transition probabilities, and the quantities $P^{x}_{0}(\nu)$, $P_{\epsilon}(\nu)$, $P_{1-\epsilon}(\nu)$, and $P_{1}(\nu)$ are the probabilities that the reaction coordinate reaches the point $y=0$, $y=\epsilon$, $y=1-\epsilon$ and $y=1$ after $\nu$ consecutive transitions respectively. The external fluxes in Eqs. (34a) and (34d) are the probability fluxes from the remaining reactions which originate from the nodes $y=0$ and $y=1$ in the network of Figure 1.a.

To simplify Eq. (34) we perform adiabatic elimination Pigolotti and Vulpiani (2008) of the intermediate states $y=\epsilon$ and $y=1-\epsilon$: we impose the steady state regime $P_{\epsilon}^{x}(\nu+1)=P_{\epsilon}^{x}(\nu)$ and $P_{\epsilon}^{x}(\nu+1)=P_{\epsilon}^{x}(\nu)$ in Eqs.(34b) and (34c) respectively, we solve Eqs.(34b)-(34c) for $P_{\epsilon}^{x}(\nu)$ and $P_{1-\epsilon}^{x}(\nu)$, and we finally substitute the result back into Eqs.(34a), (34d). This yield the effective Markov chain

[TABLE]

where we have defined the effective probabilities $p_{\rightarrow}^{x}$ and $p_{\leftarrow}^{x}$ to incorporate a monomer ($h\rightarrow hx$) and remove a monomer ($h\leftarrow hx$) respectively as

[TABLE]

Substituting Eq.(6)-(8) into Eq. (36) and then expanding for small $\epsilon$ finally gives Eq. (9).

Appendix C Derivation of the polymerization speed via reaction coordinates.

To derive the polymerization speed, we consider a mean field formulation of the polymerization process in Figure 1.a where the enzyme can remove any monomer in the copy heteropolymer. Removal of $r$ and $w$ monomers occurs with probabilities $1-\eta$ and $\eta$ respectively. This assumption simplifies the reaction tree of Figure 1.a into the closed network of Fig. 7.a, where the incorporation and removal probabilities $p_{\rightarrow}^{x}$ and $p_{\leftarrow}^{x}$ are defined as in Eq. (9).

We now introduce the reaction coordinate in this mean field description, Fig. 7.b. For later convenience, we also consider the values of the reaction coordinate $y=0$ $y=\epsilon$, $y=1-\epsilon$ and $y=1$ together with the probabilities $\pi_{\tilde{y},y}^{x}$s defined in Eqs.(6) and (II.2).

Using the scheme in Figure 7.b, we define the probability $P_{0,1}(\zeta)$ that $y=0$ or $y=1$ after $\zeta$ consecutive transitions, and the probabilities $P_{\epsilon}^{r}(\zeta)$, $P_{1-\epsilon}^{r}(\zeta)$, $P_{\epsilon}^{w}(\zeta)$, $P_{1-\epsilon}^{w}(\zeta)$ that $y=\epsilon$ or $y=1-\epsilon$ for the $r$ and $w$ monomer after $\zeta$ consecutive transitions. These probabilities evolves according to the Markov chain

[TABLE]

where

[TABLE]

and

[TABLE]

where $c^{\prime}=[\pi_{\epsilon,0}^{r}+(1-\eta)\pi_{1-\epsilon,1}^{r}+\pi_{\epsilon,0}^{w}+\eta\pi_{1-\epsilon,1}^{w}]$ is a normalization constant. We now define the matrices

[TABLE]

which contain the contribution of each transition to the heteropolymer length $N$ and the polymerization time $\mathcal{T}$. The time increments $\langle\mathrm{d}\tau\rangle_{y,\tilde{y}}^{x}$ in $\mathtt{J}^{\mathcal{T}}$ are the first passage times from $\tilde{y}$ to $y$ Gardiner (2009). In particular we have that

[TABLE]

with

[TABLE]

The remaining first passage times $\langle\mathrm{d}\tau\rangle_{\epsilon,0}^{r}$, $\langle\mathrm{d}\tau\rangle_{1-\epsilon,1}^{r}$, $\langle\mathrm{d}\tau\rangle_{\epsilon,0}^{w}$ and $\langle\mathrm{d}\tau\rangle_{1-\epsilon,1}^{w}$, are assumed equal to zero for simplicity. Physically, this assumption is justified when binding and unbinding of monomers is much faster than processing a monomer into a finalized incorporation.

Using Eq.(40) we define the tilted matrix $\mathtt{B}$ with components

[TABLE]

and dummy variables $q_{N}$, and $q_{\tau}$. For large values of $\zeta$, the largest eigenvalue of $\mathtt{B}$ coincides with the scaled cumulant generating function of $N$ and $\mathcal{T}$, see Touchette (2009). The implicit function theorem then implies

[TABLE]

where $\det\left[\mathtt{B}-\lambda\mathtt{I}\right]$ is the characteristic polynomial of $\mathtt{B}$. To compute $v$ we finally use that

[TABLE]

which is equivalent to Eq.(18). Substituting Eqs.(6), (II.2), (39), (40) and (43) into Eqs.(45) and then taking the leading order for small $\epsilon$ yields Eq. (20), where

[TABLE]

and $c$ is defined as in Eq. (13).

Appendix D Effetive incorporation and removal probabilities for the Kinetic proofreading example

To compute the incorporation and removal probabilities for the Kinetic proofreading case we mimic the procedure that leads to Eq.(9). We consider the probabilities $P_{h}(\xi)$, $P_{hx^{*}}(\xi)$ and $P_{hx}(\xi)$ to obtain the reactants $h$, $hx^{*}$, and $hx$ after $\xi$ sub-reactions of Eq. (25). These probabilities evolve according to the Markov chain

[TABLE]

where the external fluxes are the probability fluxes of the other sub-reactions entering the nodes $y=0$ and $y=1$. At steady state, we simplify Eq. (47) with adiabatic eliminationPigolotti and Vulpiani (2008): we impose $P_{hx^{*}}(\xi+1)=P_{hx^{*}}(\xi)$ into Eq. (47b), solve it for $P_{hx^{*}}(\xi)$ and substitute the solution in Eqs. (47a) and (47c). This yields, after some rearrangements,

[TABLE]

with effective incorporation/removal probabilities $p_{\rightarrow}^{x}$ and $p_{\leftarrow}^{x}$ defined as in Eq.(30).