We model the assembly of a fixed number of well-defined target structures from limited resources. Specifically, we consider a set of $S$ different species of constituents denoted by $1,\mathrm{\dots },S$ which assemble into rings of size $L$. The cases $S=1$ and $1<S\le L$ ($S=L$) are denoted as homogeneous and partially (fully) heterogeneous, respectively. The homogeneous model builds on previous work on virus capsid (Chen et al., 2008; Hagan et al., 2011), linear protein filament assembly (Michaels et al., 2016; Michaels et al., 2017; D'Orsogna et al., 2012) and aggregation and polymerization models (Krapivsky et al., 2010). The heterogeneous model in turn links to previous model systems used to study, for example, DNA-brick-based assembly of heterogeneous structures (Murugan et al., 2015; Hedges et al., 2014; D'Orsogna et al., 2013). We emphasize that, even though strikingly similar experimental realizations of our model exist (Gerling et al., 2015; Wagenbauer et al., 2017; Praetorius and Dietz, 2017), it is not intended to describe any particular system. The ring structure represents a general linear assembly process involving building blocks with equivalent binding properties and resulting in a target of finite size. The main assumption in the ring model is that the different constituents assemble linearly in a sequential order. In many biological self-assembling systems like bacterial flagellum assembly or biogenesis of the ribosome subunits the assumption of a linear binding sequence appears to be justified (Peña et al., 2017; Chevance and Hughes, 2008). In order to test the validity of our results beyond these constraints we also perform stochastic simulations of generalized self-assembling systems that do not obey a sequential binding order: i) by explicitly allowing for polymer-polymer bindings and ii) by considering the assembly of finite sized squares that grow independently in two dimensions (see Figures 6 and 7).

The assembly process starts with $N$ inactive monomers of each species. We use $C=N/V$ to denote the initial concentration of each monomer species, where $V$ is the reaction volume. Monomers are activated independently at the same per capita rate $\alpha$, and, once active, are available for binding. Binding takes place only between constituents of species with periodically consecutive indices, for example 1 and 2 or $S$ and 1 (leading to structures such as $\mathrm{\dots }1231\mathrm{\dots }$ for $S=3$); see Figure 1. To avoid ambiguity, we restrict ring sizes to integer multiples of the number of species $S$. Furthermore, we neglect the possibility of incorrect binding, for example species 1 binding to 3 or $S-1$. Polymers, that is incomplete ring structures, grow via consecutive attachment of monomers. For simplicity, polymer-polymer binding is disregarded at first, as it is typically assumed to be of minor importance (Zlotnick et al., 1999; Chen et al., 2008; Murugan et al., 2015; Haxton and Whitelam, 2013). To probe the robustness of the model, later we consider an extended model including polymer-polymer binding for which the results are qualitatively the same (see Figure 6 and the discussion). Furthermore, it has been observed that nucleation phenomena play a critical role for self-assembly processes (Ke et al., 2012; Wei et al., 2012; Reinhardt and Frenkel, 2014; Chen et al., 2008). So it is in general necessary to take into account a critical nucleation size, which marks the transition between slow particle nucleation and the faster subsequent structure growth (Michaels et al., 2016; Lazaro and Hagan, 2016; Morozov et al., 2009; Murugan et al., 2015). We denote this critical nucleation size by $L_{\mathrm{nuc}}$, which in terms of classical nucleation theory corresponds to the structure size at which the free energy barrier has its maximum. For $l<L_{\mathrm{nuc}}$ attachment of monomers to existing structures and decay of structures (reversible binding) into monomers take place at size-dependent reaction rates ${\mu }_l$ and ${\delta }_l$, respectively (Figure 1). Here, we focus on identical rates ${\mu }_l=\mu$ and ${\delta }_l=\delta$. A discussion of the general case is given in Appendix 4. Above the nucleation size, polymers grow by attachment of monomers with reaction rate $\nu \ge \mu$ per binding site. As we consider successfully nucleated structures to be stable on the observational time scales, monomer detachment from structures above the critical nucelation size is neglected (irreversible binding) (Murugan et al., 2015; Chen et al., 2008). Complete rings neither grow nor decay (absorbing state).

Figure 1

Download asset Open asset

We investigate two scenarios for the control of nucleation speed, first separately and then in combination. For the ‘activation scenario’ we set $\mu =\nu$ (all binding rates are equal) and control the assembly process by varying the activation rate $\alpha$. For the ‘dimerization scenario’ all particles are inherently active ($\alpha \to \mathrm{\infty }$) and we control the assembly process by varying the dimerization rate $\mu$ (we focus on $L_{\mathrm{nuc}}=2$). It has been demonstrated previously in Chen et al. (2008) and (Endres and Zlotnick, 2002; Hagan and Elrad, 2010; Morozov et al., 2009) that either a slow activation or a slow dimerization step are suitable in principle to retard nucleation and favour growth of the structures over the initiation of new ones. We quantify the quality of the assembly process in terms of the assembly yield, defined as the number of successfully assembled ring structures relative to the maximal possible number $NS/L$. Yield is measured when all resources have been used up and the system has reached its final state. We do not discuss the assembly time in this manuscript, however, in Appendix 5 we show typical trajectories for the time evolution of the yield in the activation and dimerization scenario. If the assembly product is stable (absorbing state), the yield can only increase with time. Consequently, the final yield constitutes the upper limit for the yield irrespective of additional time constraints. Therefore, the final yield is an informative and unambiguous observable to describe the efficiency of the assembly reaction.

We simulated our system both stochastically via Gillespie’s algorithm (Gillespie, 2007) and deterministically as a set of ordinary differential equations corresponding to chemical rate equations (see Appendix 1).

First, we consider the macroscopic limit, $N\gg 1$, and investigate how assembly yield depends on the activation rate $\alpha$ (activation scenario) and the dimerization rate $\mu$ (dimerization scenario) for $L_{\mathrm{nuc}}=2$. Here, the deterministic description coincides with the stochastic simulations (Figure 2a and b). For both high activation and high dimerization rates, yield is very poor. Upon decreasing either the activation rate (Figure 2a) or the dimerization rate (Figure 2b), however, we find a threshold value, ${\alpha }_{\mathrm{th}}$ or ${\mu }_{\mathrm{th}}$ , below which a rapid transition to the perfect yield of 1 is observed both in the deterministic and stochastic simulation. By exploiting the symmetries of the system with respect to relabeling of species, one can show that, in the deterministic limit, the behavior is independent of the number of species $S$ (for fixed $L$ and $N$, see Appendix 1). Consequently, all systems behave equivalently to the homogeneous system and yield becomes independent of $S$ in this limit. Note, however, that equivalent systems with differing $S$ have different total numbers of particles $SN$ and hence assemble different total numbers of rings.

Figure 2

Download asset Open asset

Decreasing the activation rate reduces the concentration of active monomers in the system. Hence growth of the polymers is favored over nucleation, because growth depends linearly on the concentration of active monomers while nucleation shows a quadratic dependence. Likewise, lower dimerization rates slow down nucleation relative to growth. Both mechanisms therefore restrict the number of nucleation events, and ensure that initiated structures can be completed before resources become depleted (see Figure 2c and d).

Mathematically, the deterministic time evolution of the polymer size distribution $c(l,t)$ is described by an advection-diffusion equation (Endres and Zlotnick, 2002; Yvinec et al., 2012) with advection and diffusion coefficients depending on the instantaneous concentration of active monomers (see Appendix 2). Solving this equation results in the wavefront of the size distribution advancing from small to large polymer sizes (Figure 2e). Yield production sets in as soon as the distance travelled by this wavefront reaches the maximal ring size $L$. Exploiting this condition, we find that in the deterministic system for $L_{\mathrm{nuc}}=2$, a non-zero yield is obtained if either the activation rate or the dimerization rate remains below a corresponding threshold value, that is if $\alpha <{\alpha }_{\mathrm{th}}$ or $\mu <{\mu }_{\mathrm{th}}$, where

(see Appendix 3) with proportionality constants ${\displaystyle P_{\alpha }=[\sqrt{\pi }\mathrm{\Gamma }(2/3)/\mathrm{\Gamma }(7/6){]}^3/3\approx 5.77}$ and $P_{\mu }={\pi }^2/2\approx 4.93$. These relations generalize previous results (Morozov et al., 2009) to finite activation rates and for heterogeneous systems. A comparison between the threshold values given by Equation 1 and the simulated yield curves is shown in Figure 2a,b. The relations highlight important differences between the two scenarios (where $\alpha \to \mathrm{\infty }$ and $\mu =\nu$, respectively): While ${\alpha }_{\mathrm{th}}$ decreases cubically with the ring size $L$, ${\mu }_{\mathrm{th}}$ does so only quadratically. Furthermore, the threshold activation rate ${\alpha }_{\mathrm{th}}$ increases with the initial monomer concentration $C$. Consequently, for fixed activation rate, the yield can be optimized by increasing $C$. In contrast, the threshold dimerization rate is independent of $C$ and the yield curves coincide for $N\gg 1$. Finally, if $\alpha$ is finite and $\mu <\nu$, the interplay between the two slow-nucleation scenarios may lead to enhanced yield. This is reflected by the factor $\nu /\mu$ in ${\alpha }_{\mathrm{th}}$, and we will come back to this point later when we discuss the stochastic effects.

In summary, for large particle numbers ($N\gg 1$), perfect yield can be achieved in two different ways, independently of the heterogeneity of the system - by decreasing either the activation rate (activation scenario) or the dimerization rate (dimerization scenario) below its respective threshold value.

Next, we consider the limit where the particle number becomes relevant for the physics of the system. In the activation scenario, we find a markedly different phenomenology if resources are sparse. Figure 3a shows the dependence of the average yield on the activation rate for different, low particle numbers in the completely heterogeneous case ($S=L$). Here, we restrict our discussion to the average yield. The error of the mean is negligible due to the large number of simulations used to calculate the average yield. Still, due to the randomness in binding and activation, the yield can differ between simulations. A figure with the average yield and its standard deviation is shown in Appendix 6. For very low and very high average yield, the standard deviation has to be small due to the boundedness of the yield. For intermediate values of the average, the standard deviation is highest but still small compared to the average yield. Thus, the average yield is meaningful for the essential understanding of the assembly process. Whereas the deterministic theory predicts perfect yield for small activation rates, in the stochastic simulation yield saturates at an imperfect value $y_{\max }<1$. Reducing the particle number $N$ decreases this saturation value $y_{\max }$ until no finished structures are produced ($y_{\max }\to 0$). The magnitude of this effect strongly depends on the size of the target structure $L$ if the system is heterogeneous. Figure 3c shows a diagram characterizing different regimes for the saturation value of the yield, $y_{\max }(N,L)$, in dependence of the particle number $N$ and the size of the target structure $L$ for fully heterogeneous systems $(S=L)$. We find that the threshold particle number $N_y^{th}$ necessary to obtain a fixed yield $y$ increases nonlinearly with the target size $L$. For the depicted range of $L$, the dependence of the threshold for nonzero yield, $N_{>0}^{th}$, on $L$ can approximately be described by a power-law: $N_{>0}^{th}\sim L^{\xi }$, with exponent $\xi \approx 2.8$ for $L\le 600$. Consequently, for $L=600$ already more than 10⁵ rings must be assembled in order to obtain a yield larger than zero. In Appendix 8 we included two additional plots that show the dependence of $y_{\max }$ on $N$ for fixed $L$ and the dependence on $L$ for fixed $N$, respectively. The suppression of the yield is caused by fluctuations (see explanation below) and is not captured by a deterministic description. Because these stochastic effects can decrease the yield from a perfect value in a deterministic description to zero (see Figure 3a), we term this effect ‘stochastic yield catastrophe’. For fixed target size $L$ and fixed maximum number of target structures $\frac{NS}{L}$, $y_{\max }$ increases with decreasing number of species, see Figure 3d. In the fully homogeneous case, $S=1$, a perfect yield of 1 is always achieved for $\alpha \to 0$. The decrease of the maximal yield with the number of species $S$ thus suggests that, in order to obtain high yield, it is beneficial to design structures with as few different species as possible. In large part this effect is due to the constraint $SN=\text{const}$, whereby the more homogeneous systems (small $S$) require larger numbers of particles per species $N$ and, correspondingly, exhibit less stochasticity. If $N$ is fixed instead of $SN$, the yield still initially decreases with increasing number of species $S$ but then quickly reaches a stationary plateau and gets independent of $S$ for $S\gg 1$, see Appendix 7. Moreover, increasing the nucleation size $L_{\mathrm{nuc}}$, and with it the reversibility of binding, also increases $y_{\max }$, see Figure 3(d). This indicates that, beside heterogeneity of the target structure, irreversibility of binding on the relevant time scale makes the system susceptible to stochastic effects.

Figure 3

Download asset Open asset

The stochastic yield catastrophe is mainly attributable to fluctuations in the number of active monomers. In the deterministic (mean-field) equation the different particle species evolve in balanced stoichiometric concentrations. However, if activation is much slower than binding, the number of active monomers present at any given time is small, and the mean-field assumption of equal concentrations is violated due to fluctuations (for $S>1$). Activated monomers then might not fit any of the existing larger structures and would instead initiate new structures. Figure 4a illustrates this effect and shows how fluctuations in the availability of active particles lead to an enhanced nucleation and, correspondingly, to a decrease in yield. Due to the effective enhancement of the nucleation rate, the resulting polymer size distribution has a higher amplitude than that predicted deterministically (Figure 4b) and the system is prone to depletion traps. A similar broadening of the size distribution has been reported in the context of stochastic coagulation-fragmentation of identical particles (D'Orsogna et al., 2015).

Figure 4

Download asset Open asset

In the dimerization scenario, in contrast, there is no stochastic activation step. All particles are available for binding from the outset. Consequently, stochastic effects do not play an essential role in the dimerization scenario and perfect yield can be reached robustly for all system sizes, regardless of the number of species $S$ (Figure 3(b)).

So far, the two implementations of the ‘slow nucleation principle’ have been investigated separately. Surprisingly, we observe counter-intuitive behavior in a mixed scenario in which both dimerization and activation occur slowly (i.e., $\mu <\nu$, $\alpha <\mathrm{\infty }$). Figure 5 shows that, depending on the ratio $\mu /\nu$, the yield can become a non-monotonic function of $\alpha$. In the regime where $\alpha$ is large, nucleation is dimerization-limited; therefore activation is irrelevant and the system behaves as in the dimerization scenario for $\alpha \to \mathrm{\infty }$. Upon decreasing $\alpha$ we then encounter a second regime, where activation and dimerization jointly limit nucleation. The yield increases due to synergism between slow dimerization and activation (see $\mu /\nu$ dependence of ${\alpha }_{\mathrm{th}}$, Equation 1), whilst the average number of active monomers is still high and fluctuations are negligible. Finally, a stochastic yield catastrophe occurs if $\alpha$ is further reduced and activation becomes the limiting step. This decline is caused by an increase in nucleation events due to relative fluctuations in the availability of the different species (‘fluctuations between species’). This contrasts the deterministic description where nucleation is always slower for smaller activation rate. Depending on the ratio $\mu /\nu$, the ring size $L$ and the particle number $N$, maximal yield is obtained either in the dimerization-limited (red curves, Figure 5), activation-limited (blue curve, Figure 5b) or intermediate regime (green and orange curves, Figure 5).

Figure 5

Download asset Open asset

In our model, the reason for the stochastic yield catastrophe is that - due to fluctuations between species - the effective nucleation rate is strongly enhanced. Hence, if binding to a larger structure is temporarily impossible, activated monomers tend to initiate new structures, causing an excess of structures that ultimately cannot be completed. Natural questions that arise are whether (i) relaxing the constraint that polymers cannot bind other polymers or (ii) abandoning the assumption of a linear assembly path, will resolve the stochastic yield catastrophe. To answer these questions, we performed stochastic simulations for extensions of our model system showing that the stochastic yield catastrophe indeed persists. We start by considering the ring model from the previous section but take polymer-polymer binding into account in addition to growth via monomer attachment (Figure 6). In detail, we assume that two structures of arbitrary size (and with combined length $\le L$) bind at rate $\nu$ if they fit together, that is if the left (right) end of the first structure is periodically continued by the right (left) end of the second one. Realistically, the rate of binding between two structures is expected to decrease with the motility and thus the sizes of the structures. In order to assess the effect of polymer-polymer binding, we focus on the worst case where the rate for binding is independent of the size of both structures. If a stochastic yield catastrophe occurs for this choice of parameters, we expect it to be even more pronounced in all the ‘intermediate cases’. Figure 6 shows the dependence of the yield on the activation rate in the polymer-polymer model. As before, yield increases below a critical activation rate and then saturates at an imperfect value for small activation rates. Decreasing the number of particles per species, decreases this saturation value. Compared to the original model, the stochastic yield catastrophe is mitigated but still significant: For structures of size $S=L=100$, yield saturates at around 0.87 for $N=100$ particles per species and at around 0.33 for $N=10$ particles per species. We thus conclude that polymer-polymer binding indeed alleviates the stochastic yield catastrophe but does not resolve it. Since binding only happens between consecutive species, structures with overlapping parts intrinsically can not bind together and depletion traps continue to occur. Taken together, also in the extended model, fluctuations in the availability of the different species lead to an excess of intermediate-sized structures that get kinetically trapped due to structural mismatches. Note that in the extreme case of $N=1$, incomplete polymers can always combine into one final ring structure so that in this case the yield is always 1. Analogously, for high activation rates yield is improved for $N=10$ compared to $N\ge 50$ (Figure 6b).

Figure 6

Download asset Open asset

Kinetic trapping due to structural mismatches can occur in every (partially) irreversible heterogeneous assembly process with finite-sized target structure and limited resources. From our results, we thus expect a stochastic yield catastrophe to be common to such systems. In order to further test this hypothesis, we simulated another variant of our model where finite sized squares assemble via monomer attachment from a pool of initially inactive particles, see Figure 7. In contrast to the original model, the assembled structures are non-periodic and exhibit a non-linear assembly path where structures can grow independently in two dimensions. While the ring model assumes a sequential order of binding of the monomers, the square allows for a variety of distinct assembly paths that all lead to the same final structure. Note that, because of the absence of periodicity, the square model is only well defined for the completely heterogeneous case. Figure 7 depicts the dependence of the yield on the activation rate for a square of size $S=100$. Also in this case, we find that the yield saturates at an imperfect value for small activation rates. Hence, we showed that the stochastic yield catastrophe is not resolved neither by accounting for polymer-polymer combination nor by considering more general assembly processes with multiple parallel assembly paths. This observation supports the general validity of our findings and indicates that stochastic yield catastrophes are a general phenomenon of (partially) irreversible and heterogeneous self-assembling systems that occur if particle number fluctuations are non-negligible.

Figure 7

Download asset Open asset

We start with the general Master equation and derive the chemical rate equations (deterministic/mean-field equations) for the heterogeneous self-assembly process. We renounce to show the full Master equation here but instead state the system that describes the evolution of the first moments. To this end, we denote the random variable that describes the number of polymers of size $\mathrm{\ell }$ and species $s$ in the system at time $t$ by $n_{\mathrm{\ell }}^s(t)$ with $2\le \mathrm{\ell }<L$ and $1\le s\le S$. The species of a polymer is defined by the species of the respective monomer at its left end. Furthermore, ${\displaystyle n_0^s}$ and $n_1^s$ denote the number of inactive and active monomers of species $s$, respectively, and $n_L$ the number of complete rings. We signify the reaction rate for binding of a monomer to a polymer of size $\mathrm{\ell }$ by ${\displaystyle {\nu }_{\ell }}$. ${\displaystyle \alpha }$ denotes the activation rate and ${\delta }_{\mathrm{\ell }}$ the decay rate of a polymer of size $\mathrm{\ell }$. By $⟨\mathrm{\dots }⟩$ we indicate (ensemble) averages. The system governing the evolution of the first moments (the averages) of the $\{n_{\mathrm{\ell }}^s\}$ is then given by:

The different terms of this equation are illustrated graphically in Figure 8. The first equation describes loss of inactive particles due to activation at rate $\alpha$. Equation 2b gives the temporal change of the number of active monomers that is governed by the following processes: activation of inactive monomers at rate $\alpha$, binding of active monomers to the left or to the right end of an existing structure of size $\mathrm{\ell }$ at rate ${\nu }_{\mathrm{\ell }}$, and decay of below-critical polymers of size $\mathrm{\ell }$ into monomers at rate ${\delta }_{\mathrm{\ell }}$ (disassembly). Equations 2c and 2d describe the dynamics of dimers and larger polymers of size $3\le \mathrm{\ell }<L$, respectively. The terms account for reactions of polymers with active monomers (polymerization) as well as decay in the case of below-critical polymers (disassembly). The indicator function ${\mathrm{𝟏}}_{\{x<L_{\text{nuc}}\}}$ equals 1 if the condition $x<L_{\text{nuc}}$ is satisfied and 0 otherwise. Note that a polymer of size $\mathrm{\ell }\ge 3$ can grow by attaching a monomer to its left or to its right end whereas the formation of a dimer of a specific species is only possible via one reaction pathway (dimerization reaction). Finally, polymers of length $L$ – the complete ring structures – form an absorbing state and, therefore, include only the respective gain terms (cf Equation 2e).

Figure 8

Download asset Open asset

We simulated the Master equation underlying Equation 2 stochastically using Gillespie’s algorithm. For the following deterministic analysis, we neglect correlations between particle numbers $\{n_{\mathrm{\ell }}^s\}$, which is valid assumption for large particle numbers. Then the two-point correlator can be approximated as the product of the corresponding mean values (mean-field approximation)

Furthermore, for the expectation values it must hold

because all species have equivalent properties (there is no distinct species) and hence the system is invariant under relabelling of the upper index. By

we denote the concentration of any monomer or polymer species of size $\mathrm{\ell }$, where $V$ is the reaction volume. Due to the symmetry formulated in Equation 4, the heterogeneous assembly process decouples into a set of $S$ identical and independent homogeneous assembly processes in the deterministic limit. The corresponding homogeneous system then is described by the following set of equations that is obtained by applying (Equation 3, Equation 4) and (Equation 5) to (Equation 2)

The rate constants ${\nu }_{\mathrm{\ell }}$ in Equations 6 and 2 differ by a factor of $V$. For convenience, we use however the same symbol in both cases. The rate constants ${\nu }_{\mathrm{\ell }}$ in Equation 6 can be interpreted in the usual units ${\displaystyle [\frac{\text{liter}}{\text{mol sec}}]}$. Due to the symmetry, the yield, which is given by the quotient of the number of completely assembled rings and the maximum number of complete rings, becomes independent of the number of species $S$

Hence, it is enough to study the dynamics of the homogeneous system, Equation 6, to identify the condition under which non zero yield is obtained.

The dynamical properties of the evolution of the polymer-size distribution become evident if the set of ODEs, Equation 6, is rewritten as a partial differential equation. This approach was previously described in the context of virus capsid assembly (Zlotnick et al., 1999; Morozov et al., 2009). For simplicity, we restrict ourselves to the case $L_{\text{nuc}}=\mathrm{\hspace{0.17em}2}$ and let ${\nu }_1=\mu$ and ${\nu }_{\mathrm{\ell }\ge 2}=\nu$. Then, for the polymers with $\mathrm{\ell }>2$ we have

As a next step, we approximate the index $\mathrm{\ell }\in \{2,3,\mathrm{\dots },L\}$ indicating the length of the polymer as a continuous variable $x\in [2,L]$ and define $c(x=\mathrm{\ell }):=c_{\mathrm{\ell }}$. By $A:=c_1$ we denote the concentration of active monomers in the following to emphasize their special role. Formally expanding the right-hand side of Equation 8 in a Taylor series up to second order

one arrives at the advection-diffusion equation with both advection and diffusion coefficients depending on the concentration of active monomers $A(t)$

Equation 10 can be written in the form of a continuity equation ${\partial }_{t}c(x)=-{\partial }_{x}J(x)$ with flux $J=\mathrm{\hspace{0.17em}2}\nu Ac-\nu A{\partial }_{x}c$. The flux at the left boundary $x=\mathrm{\hspace{0.17em}2}$ equals the influx of polymers due to dimerization of free monomers $J(2,t)=\mu A^2$. This enforces a Robin boundary condition at $x=\mathrm{\hspace{0.17em}2}$

At $x=L$ we set an absorbing boundary $c(L,t)=\mathrm{\hspace{0.17em}0}$ so that completed structures are removed from the system. The time evolution of the concentration of active monomers is given by

The terms on the right-hand side account for activation of inactive particles, dimerization, and binding of active particles to polymers (polymerization).

Qualitatively, Equation 10 describes a profile that emerges at $x=\mathrm{\hspace{0.17em}2}$ from the boundary condition Equation 11, moves to the right with time-dependent velocity $2\nu A(t)$ due to the advection term, and broadens with a time-dependent diffusion coefficient $\nu A(t)$. In Appendices 2–3 we show how the full solution of Equations 10 and 11 can be found assuming knowledge of $A(t)$. Here, we focus only on the derivation of the threshold activation rate and threshold dimerization rate that mark the onset of non-zero yield. Yield production starts as soon as the density wave reaches the absorbing boundary at $x=L$. Therefore, finite yield is obtained if the sum of the advectively travelled distance $d_{\text{adv}}$ and the diffusively travelled distance $d_{\text{diff}}$ exceeds the system size $L-2$

According to Equation 10, $d_{\text{adv}}=2\nu \underset{0}{\overset{\mathrm{\infty }}{\int }}A(t)𝑑t$ and $d_{\text{diff}}=\sqrt{2\nu \underset{0}{\overset{\mathrm{\infty }}{\int }}A(t)𝑑t}$, giving as condition for the onset of finite yield

where the last approximation is valid for large $L$.

In order to obtain ${\int }_0^{\mathrm{\infty }}A(t)𝑑t$ we derive an effective two-component system that governs the evolution of $A(t)$. To this end, we denote the total number of polymers in Equation 12 by $B(t):={\int }_2^{\mathrm{\infty }}c(x,t)𝑑x$ (as long as yield is zero the upper boundary is irrelevant and we can consider $L=\mathrm{\infty }$). Equation 12 then reads

and the dynamics of $B$ is determined from the boundary condition, Equation 11

Measuring $A$ and $B$ in units of the initial monomer concentration $C$ and time in units of ${(\nu C)}^{-1}$ the equations are rewritten in dimensionless units as

where $\omega =\frac{\alpha }{\nu C}$ and $\eta =\frac{\mu }{\nu }$. Equation 17 describes a closed two-component system for the concentration of active monomers $A$ and the total concentration of polymers $B$. It describes the dynamics exactly as long as yield is zero. In order to evaluate the condition (14) we need to determine the integral over $A(t)$ as a function of $\omega$ and $\eta$

To that end, we proceed by looking at both scenarios separately. The numerical analysis, confirming our analytic results, is given in Appendix 3.

The activation rate in the dimerization scenario is $\alpha \to \mathrm{\infty }$, and instead of the term $\omega e^{-\omega t}$ in $\mathrm{d}A/\mathrm{d}t$, we set the initial condition $A(0)=1$ (and $B(0)=0$). Furthermore, $\eta =\mu /\nu \ll 1$ and we can neglect the term proportional to $\eta$ in $\mathrm{d}A/\mathrm{d}t$. As a result,

Solving this equation for $A$ as a function of $B$ using the initial condition $A(B=0)=1$, the totally travelled distance of the wave is determined to be

where for the evaluation of the integral we used the substitution $\eta A^2\mathrm{d}t=\mathrm{d}B$.

In the activation scenario, yield sets in only if the activation rate and thus the effective nucleation rate is slow. As a result, in addition to $\omega \ll 1$, we can again neglect the term proportional to $\eta$ in $\mathrm{d}A/\mathrm{d}t$. This time, however, we have to keep the term $\omega e^{-\omega t}$. As a next step, we assume that $\mathrm{d}A/\mathrm{d}t$ is much smaller than the remaining terms on the right-hand side, $\omega e^{-\omega t}$ and $-2AB$. This assumption might seem crude at first sight but is justified a posteriori by the solution of the equation (see Appendix 3). Hence, we get the algebraic equation $A(t)=\omega e^{-\omega t}/(2B(t))$. Using it to solve $\mathrm{d}B/\mathrm{d}t=\eta A^2$ for $B$, and then to determine $A$, the totally travelled distance of the wave is deduced as

Taken together, we therefore obtain two conditions out of which one must be fulfilled in order to obtain finite yield

where $a$ and $b$ are numerical factors, and $P_{\alpha }=\mathrm{\hspace{0.17em}8}{a}^3\approx 5.77$ and $P_{\mu }=\mathrm{\hspace{0.17em}4}{b}^2\approx 4.93$. This verifies Equation 1 in the main text.

In this section we derive the chemical rate equations (deterministic equations) for the self-assembly process as described in the main text. Furthermore, we show that for general $S$ in the deterministic limit the model is equivalent to a set of $S$ independent assembly processes with only one species.

First, we consider the homogeneous model ($S=\mathrm{\hspace{0.17em}1}$). By $c_{\mathrm{\ell }}(t)$ we denote the concentration of complexes of length $\mathrm{\ell }$ ($\mathrm{\ell }\ge 2$) at time $t$, $c_1(t)$ is the concentration of active monomers and $c_0(t)$ the concentration of inactive monomers at time $t$. In the following we will usually skip the time argument for better readability. We denote the reaction rate for binding of a monomer to a polymer of size $\mathrm{\ell }$ by ${\nu }_{\mathrm{\ell }}$. The model from the main text is recovered by setting ${\nu }_{\mathrm{\ell }}:={\mu }_{\mathrm{\ell }}$ if $\mathrm{\ell }<L_{\text{nuc}}$, and ${\nu }_{\mathrm{\ell }}:=\nu$ otherwise. The ensuing set of ordinary differential equations then reads:

The indicator function ${\mathrm{𝟏}}_{\{x<L_{\text{nuc}}\}}$ equals 1 if the condition $x<L_{\text{nuc}}$ is satisfied and 0 otherwise. The first equation describes loss of inactive particles due to activation at rate $\alpha$. It is uncoupled from the remainder of the equations and is solved by $c_0(t)=C{e}^{-\alpha t}$, with $C$ denoting the initial concentration of inactive monomers. The temporal change of the active monomers is governed by the following processes (Equation A1b): activation of inactive monomers at rate $\alpha$, binding of active monomers to existing structures at rate ${\nu }_{\mathrm{\ell }}$ (polymerization), and decay of below-critical polymers into monomers at rate ${\delta }_{\mathrm{\ell }}$ (disassembly). All binding rates appear with a factor of 2 because a monomer can attach to a polymer on its left or on its right end.

Note that there is a subtlety with the dimerization term ${\displaystyle 2{\nu }_1{c}_1^2}$ in Equation A1b: the dimerization term as well bears a factor of 2 because two identical monomers $A$ and $B$ can form a dimer in two possible ways, either as $AB$ or $BA$. Additionally, there is a stoichiometric factor of 2 for the monomers in this reaction. However, one factor of 2 is cancelled again because, assuming there are $n$ monomers, the number of ordered pairs of monomers that describe possible reaction partners is $\frac{1}{2}n(n-1)\approx n^2/2$ (if $n$ is large) rather than $n^2$ (the number of reaction partners when two different species react). This leaves us with a single factor of 2 like for all the other binding reactions.

Equations A1c and A1d describe the dynamics of dimers and larger polymers of size $3\le \mathrm{\ell }<L$, respectively. The terms account for reactions of polymers with active monomers (polymerization) as well as decay in the case of below-critical polymers (disassembly). The dimerization term in the equation for ${\partial }_t{c}_2$ lacks the factor of 2 because the stoichiometric factor is missing for the dimers as compared with the dimerization term for the monomers in the line above. Finally, polymers of length $L$ – the complete ring structures – form an absorbing state and therefore only include a reactive gain term (Equation A1e).

Next we consider systems with more than one particle species ($S>1$). The heterogeneous system can be described by dynamical equations equivalent to the homogeneous system. We show this starting from a full description that distinguishes both monomers and polymers into a set of different species $1,\mathrm{\dots },S$. The species of a polymer is defined by the species of the respective monomer at its left end. As polymers assemble in consecutive order of species, a polymer is uniquely determined by its length and species (i.e. species of leftmost monomer). In that sense, $c_{\mathrm{\ell }}^s$ with $0\le \mathrm{\ell }<L$ and $1\le s\le S$ denotes the concentration of a polymer of length $\mathrm{\ell }$ and species $s$ ($c_0^s$ and $c_1^s$ again denote inactive and active monomers of species $s$, respectively). For example, $c_4^5$ denotes the concentration of polymers [5678] if $S\ge 8$, or of polymers [5612] if $S=6$. Upper indices are always assumed to be taken modulo $S$ whenever they lie outside the range $[1,S]$. Therefore, the dynamics of the concentrations $c_{\mathrm{\ell }}^s$ with $3\le \mathrm{\ell }<L$ is given by

The terms on the right-hand side account for the influx due to binding of the respective polymers of length $\mathrm{\ell }-1$ with a monomer either on the right or on the left (first and second term), and for the outflux due to reactions of a polymer of length $\mathrm{\ell }$ and species $s$ with a monomer on the right or on the left (third and fourth term), as well as for decay into monomers for $\mathrm{\ell }<L_{\text{nuc}}$ (last term). For the dynamics of the dimers, however, there is only one gain term arising from dimerization:

Equivalently, for the active monomers we find:

Now we exploit the symmetry of the system with respect to the species index, that is, the upper index in $\{c_{\mathrm{\ell }}^s\}$: Since all species in the system are equivalent, the dynamic equations are invariant under relabelling of the upper indices. Consequently, it must hold that:

In other words, the upper index is irrelevant and can also be discarded. The variable $c_{\mathrm{\ell }}$ then denotes the concentration of any one polymer species of length $\mathrm{\ell }$. Taking advantage of this symmetry for the equations of the heterogeneous system, (Equation A2, Equation A3 and Equation A4), and collecting equal terms leads to a set of equations fully identical to those for the homogeneous system (Equation A1). We show the equivalence to the homogeneous model exemplarily for the dynamics of the polymers with size $\mathrm{\ell }\ge 3$ in Equation A2. Applying $c_{\mathrm{\ell }}^s(t)=c_{\mathrm{\ell }}(t)$ to Equation A2 yields for the dynamics of the concentration of an arbitrary polymer species of size $\mathrm{\ell }$:

which is identical to the respective dynamic Equation A1d for the homogeneous model. The other equations for the heterogeneous system reduce to those for the homogeneous system in an analogous manner.

Summarizing, we have shown that the (deterministic) heterogeneous assembly process decouples into a set of $S$ identical and independent homogeneous processes. In particular, yield, which is given by the quotient of the number of completely assembled rings and the maximal possible number of complete rings, becomes independent of $S$:

The dynamical properties of the evolution of the polymer size distribution become evident if the set of ODEs, Equation 1, is rewritten as a partial differential equation. This approach was previously described in the context of virus capsid assembly (Morozov et al., 2009; Zlotnick et al., 1999; Endres and Zlotnick, 2002) but we will restate the essential steps here for the convenience of the reader. To this end we interpret the length index of the polymer $\mathrm{\ell }\in \{2,3,\mathrm{\dots },L\}$ as a continuous variable that we rename $x\in [2,L]$. With such a continuous description in view we write $c(x=\mathrm{\ell }):=c_{\mathrm{\ell }}$ to denote the concentration of polymers of size $\mathrm{\ell }$.

Since the active monomers play a special role, we denote their concentration in the following by $A$. For simplicity we restrict our discussion to the case $L_{\text{nuc}}=\mathrm{\hspace{0.17em}2}$ and let ${\nu }_1=\mu$ and ${\nu }_{\mathrm{\ell }\ge 2}=\nu$. Generalizations to $L_{\text{nuc}}>2$ can be done in a similar way. Then, for the polymers with $\mathrm{\ell }\ge 3$ we have:

Formally, expanding the right-hand side in a Taylor series up to second order

we arrive at an advection-diffusion equation with both advection and diffusion coefficients depending on the concentration of active monomers $A(t)$,

Equation A9 can be written in the form of a continuity equation ${\partial }_{t}c(x)=-{\partial }_{x}J(x)$ with flux $J=\mathrm{\hspace{0.17em}2}\nu Ac-\nu A{\partial }_{x}c$. The flux at the left boundary, $x=\mathrm{\hspace{0.17em}2}$, equals the influx of polymers due to dimerization of free monomers, $J(2,t)=\mu A^2$. This enforces a Robin boundary condition at $x=\mathrm{\hspace{0.17em}2}$,

At $x=L$, we have an absorbing boundary $c(L,t)=\mathrm{\hspace{0.17em}0}$ so that completed structures are removed from the system. Furthermore, the time evolution of the concentration of active particles is given by

The terms on the right-hand side account for activation of inactive particles, dimerization, and binding of active particles to polymers (polymerization).

Qualitatively, Equation A9 describes a profile that emerges at $x=\mathrm{\hspace{0.17em}2}$ from the boundary condition, Equation A10, moves to the right with time dependent velocity $2\nu A(t)$ due to the advection term, and broadens with a time-dependent diffusion coefficient $\nu A(t)$. The concentration of active particles $A$ determines both the influx of dimers at $x=\mathrm{\hspace{0.17em}2}$, as well as the speed and diffusion of the wave profile.

Next, we derive an expression that solves Equation A9, assuming that we know $A(t)$. We start by solving Equation A9 at the left boundary $c(2,t)$, and then translate the resulting expression to obtain a solution for $c(x,t)$. To obtain $c(2,t)$ in dependence of $a(t)$ we can solve ${\displaystyle \frac{d}{dt}c(2,t)=\mu A^2-2\nu Ac(2,t)}$ (see Equation A1c) by ’variation of the constants’ as

With help of this expression we find $c(x,t)$: Given $c(2,t)$, the advective part of Equation A9,

is solved by

Here, $\tau (x,t)$ denotes the time when a particle now at position $x$ and time $t$ was at $x=2$. In other words, a particle at time $t$ and position $x$ has entered the system at $x=2$ at time $\tau (x,t)$. This ansatz solves the PDE (Equation A13) if and only if $\tau (x,t)$ satisfies

with $\tilde{A}$ being an arbitrary integral of $A$ such that ${\partial }_t\tilde{A}(t)=A(t)$ and ${\tilde{A}}^{-1}$ denoting its inverse. More easily, we find this form of $\tau$ by requiring that the integral over the velocity from time $\tau$ to $t$ equals the travelled distance $x-2$:

To include the diffusive contribution in Equation A13, we use the diffusion kernel,

with the time dependent diffusion constant $D(t)=\nu A(t)$. The kernel $k(x,y,t)$ accounts for the mass that has been diffusively transported from $y$ over a distance of $x$. Because the mass has entered the system at $x=2$ at time $\tau (y,t)$, it diffused for the time $t-\tau (y,t)$. The complete expression for $c(x,t)$ is then obtained as the convolution of $c_{\text{advec}}(x,t)$ (Equation A14), that is obtained from Equation A12 and Equation A15, and the diffusion kernel $k(x,y,t)$ (Equation A17):

Interpreting the terms in the equations and the general form of the solution, we are able to understand the qualitative behavior of the system. If both the activation and the dimerization rate are large, the system produces zero yield: both advection and diffusion are driven by the concentration of active monomers $A$. If activation is fast, the concentration of active monomers $A$ will become large initially since activation is faster than the reaction dynamics. Consequently, provided $\mu \sim \nu$, dimerization dominates over binding because it depends quadratically on $A$, see Equation A11. The reservoir of free particles then depletes quickly and cannot sustain the motion of the wave for long enough to reach the absorbing boundary, resulting in a very low yield. Only if either the activation rate is low enough or if $\mu \ll \nu$, the motion of the wave can be sustained until it reaches the absorbing boundary.

Based on the analysis from the previous section, we will now determine the threshold activation rate and threshold dimerization rate which mark the onset of non-zero yield. Yield production starts as soon as the density wave reaches the absorbing boundary at $x=L$. Therefore, finite yield is obtained if and only if the sum of the advectively travelled distance $d_{\text{adv}}$ and the diffusively travelled distance $d_{\text{diff}}$ exceeds the system size $L-2$:

The condition for the onset of non-zero yield is obtained by assuming equality in this relation. The advectively travelled distance is obtained from Equation A16 by setting the borders of the integral over the velocity to $\tau =0$ and $t=\mathrm{\infty }$:

The diffusively travelled distance is approximately given by the standard deviation of the Gaussian diffusion kernel, Equation A17, again with $\tau =0$ and $t=\mathrm{\infty }$,

Taken together, we obtain a condition for the onset of finite yield:

Substituting $y=\sqrt{2\nu \int A}$ and requiring that $y$ is positive, we solve the quadratic equation and find that Equation A22 is equivalent to

where the last approximation is valid for large $L$.

We determine the threshold values for the activation rate $\alpha$ and the dimerization rate $\mu$ by finding solutions of the dynamical equation for the active particles $A(t)$, Equation A11, such that the condition, Equation A23, is fulfilled. Thus, we start by deriving the dependence of ${\int }_0^{\mathrm{\infty }}A(t)𝑑t$ on $\alpha$ and $\mu$.

The concentration $c(x,t)$ appears in Equation A11 only in terms of an integral ${\int }_2^{L}c(x,t)𝑑x$, counting the total number of polymers in the system. As long as yield is zero there is no outflux of polymers at the absorbing boundary $x=L$ and the total number of polymers in the system only increases due to the influx at the left boundary $x=\mathrm{\hspace{0.17em}2}$. As long as yield is zero we can therefore equivalently consider the limit $L\to \mathrm{\infty }$. We denote the total number of polymers in Equation A11 by $B(t):=\int c(x,t)𝑑x$ for which the dynamics is determined from the boundary condition, Equation A10:

Hence, as long as yield is zero, the total number of polymers increases with the rate of the dimerization events. The system then simplifies to a set of two coupled ordinary differential equations for $A$ and $B$:

The dynamics of $A$ and $B$ is equivalent to a two-state activator-inhibitor system, where $A$ dimerizes into $B$ at rate $\mu$, and $B$ degrades (inhibits) $A$ at rate $2\nu$. Note that Equation A25 describes the exact dynamics of the active monomers $A$ and total number of polymers $B$ in the deterministic system as long as yield is zero. The system has therefore been greatly reduced from originally $SN$ coupled ODEs to now only two coupled ODEs.

For the further analysis it is useful to non-dimensionalize Equation A25 by measuring $A$ and $B$ in units of the initial concentration of inactive monomers $C$ and time in units of ${(\nu C)}^{-1}$:

with the remaining dimensionless parameters $\omega =\frac{\alpha }{\nu C}$ and $\eta =\frac{\mu }{\nu }$. We are interested in the integral over $A(t)$ as a function of $\omega$ and $\eta$,

which relates to the totally travelled distance of the wave. Note that, in case of zero yield, $2g(\omega ,\eta )$ is the total advectively travelled distance of the wave (cf. Equation A20) and the square of the diffusively travelled distance (cf. Equation A21).

The dimerization scenario is characterized by fast activation $\alpha \gg C\nu$ and slow dimerization $\mu \ll \nu$. For the dimensionless parameters these assumptions translate to $\eta \ll 1$ and $\eta \ll \omega$. Because for small $\eta \ll 1$ nucleation is much slower than growth we neglect the dimerization term in Equation A26a against the growth term. Furthermore, because $\eta \ll \omega$ activation happens on a fast time scale compared with nucleation and we may therefore integrate out the fast time scale assuming that all particles are activated instantaneously at the beginning. The system Equation A26 then reduces to

with the initial condition $A(0)=1$ and $B(0)=0$. We divide the first equation by the second one (formally applying the chain rule and the inverse function theorem) to obtain a single equation for the dynamics of $A(B)$:

where $A(B=0)=1$. This first order ODE can be solved by separation of variables and subsequent integration, yielding

Because the number of active monomers $A(t)$ must vanish for $t\to \mathrm{\infty }$, the final value of $B$ is

Thereby, we calculate the function $g(\eta )$ via variable substitution $dt=\frac{dB}{\eta A^2}$:

So, the dependence of the travelled distance of the wave on $\eta$ obeys a power law with exponent ${\displaystyle -\frac{1}{2}}$, confirming the previous result (Morozov et al., 2009). For the coefficient we find $\frac{\pi }{2\sqrt{2}}\approx 1.1107$.

Additionally, we can determine the time dependent solutions $A(t)$ and $B(t)$. Using the solution for $A(B)$ from Equation A30 in Equation A28b we obtain $B(t)$ as

We use this expression for $B(t)$ in Equation A28a to obtain $A(t)$. The resulting ODEs can again be solved by separation of variables as

In the activation scenario, $\alpha \ll C\nu$, such that $\omega \ll 1$ and $\omega \ll \eta$. As we know already that decreasing $\omega$ will slow down nucleation relative to growth we can again neglect the dimerization term in Equation A26a. In contrast to the dimerization scenario, however, we have to keep the activation term. Transforming time via $\tau :=1-e^{-\omega t}$ such that $\tau \in [0,1]$ and writing $a(\tau )=a(1-e^{-\omega t}):=A(t)$ and $b(\tau )=b(1-e^{-\omega t}):=B(t)$ the system in Equation A26 becomes:

with the initial condition $a(0)=b(0)=0$. The function $g(\omega ,\eta )$ transforms as

In the following we derive the asymptotic solution for $a(\tau )$ in the limit of small $\omega$ in order to evaluate the integral in Equation A36. In the limit $\tau \to 1$ ($\iff t\to \mathrm{\infty }$) both $a(\tau )$ and $\frac{d}{d\tau }a(\tau )$ will become small whereas $b(\tau )$ increases monotonically. The reaction term in Equation A35a is furthermore weighted by a factor $\frac{1}{\omega }$ which will become large if $\omega \ll 1$. We therefore postulate that for sufficiently large $\tau$ the derivative $\frac{d}{d\tau }a(\tau )$ is much smaller than the two terms on the right-hand side of Equation A35a and hence negligible. This assumption has to be justified a posteriori with the obtained solution. Neglecting the derivative term $\frac{d}{d\tau }a$ in (Equation A35a) reduces the equation to an algebraic equation and we find