1. Physics of Living Systems
Download icon

Stochastic yield catastrophes and robustness in self-assembly

  1. Florian M Gartner
  2. Isabella R Graf
  3. Patrick Wilke
  4. Philipp M Geiger
  5. Erwin Frey  Is a corresponding author
  1. Arnold Sommerfeld Center for Theoretical Physics (ASC) and Center for NanoScience (CeNS), Department of Physics, Ludwig-Maximilians-Universität München, Germany
Research Article
  • Cited 0
  • Views 1,266
  • Annotations
Cite this article as: eLife 2020;9:e51020 doi: 10.7554/eLife.51020

Abstract

A guiding principle in self-assembly is that, for high production yield, nucleation of structures must be significantly slower than their growth. However, details of the mechanism that impedes nucleation are broadly considered irrelevant. Here, we analyze self-assembly into finite-sized target structures employing mathematical modeling. We investigate two key scenarios to delay nucleation: (i) by introducing a slow activation step for the assembling constituents and, (ii) by decreasing the dimerization rate. These scenarios have widely different characteristics. While the dimerization scenario exhibits robust behavior, the activation scenario is highly sensitive to demographic fluctuations. These demographic fluctuations ultimately disfavor growth compared to nucleation and can suppress yield completely. The occurrence of this stochastic yield catastrophe does not depend on model details but is generic as soon as number fluctuations between constituents are taken into account. On a broader perspective, our results reveal that stochasticity is an important limiting factor for self-assembly and that the specific implementation of the nucleation process plays a significant role in determining the yield.

eLife digest

The self-assembly of a large biological molecule from small building blocks is like finishing a puzzle of magnetic pieces by shaking the box. Even though each piece of the puzzle is attracted to its correct neighbours, the limited control makes it very hard to finish the puzzle in a short amount of time.

The problem becomes even more difficult if several copies of the same puzzle are assembled in one box. If several puzzles start at the same time, the different parts might steal pieces from each other, making it impossible to successfully complete any of the puzzles. This is called a depletion trap. If the box is only shaken and there is no real control over individual pieces, these traps occur at random.

Overcoming these random depletion traps is an important challenge when assembling nanostructures and other artificial molecules designed by humans without wasting many, potentially expensive, components. Previous studies have shown that when multiple copies of the same structure are assembled simultaneously, slowing the rate of initiation increases the yield of correctly-made structures. This prevents new structures from stealing pieces from existing structures before they are fully completed.

Now, Gartner, Graf, Wilke et al. have used a mathematical model to show that changing the way initiation is delayed leads to different yields. This was especially true for small systems where fluctuations in the availability of the different pieces strongly enhanced the initiation of new structures. In these cases, the self-assembly process terminated undesirably with many incomplete structures.

Nanostructures have various applications ranging from drug delivery to robotics. These findings suggest that in order to efficiently assemble biological molecules, the concentrations of the different building blocks need to be tightly controlled. A question for further research is to investigate strategies that reduce fluctuations in the availability of the building blocks to develop more efficient assembly protocols.

Introduction

Efficient and accurate assembly of macromolecular structures is vital for living organisms. Not only must resource use be carefully controlled, but malfunctioning aggregates can also pose a substantial threat to the organism itself (Jucker and Walker, 2013; Drummond and Wilke, 2009). Furthermore, artificial self-assembly processes have important applications in a variety of research areas like nanotechnology, biology, and medicine (Zhang, 2003; Whitesides and Grzybowski, 2002; Whitesides et al., 1991). In these areas, we find a broad range of assembly schemes. For example, while a large number of viruses assemble capsids from identical protein subunits, some others, like the Escherichia virus T4, form highly complex and heterogeneous virions encompassing many different types of constituents (Zlotnick et al., 1999; Zlotnick, 2003; Hagan, 2014; Leiman et al., 2010). Furthermore, artificially built DNA structures can reach up to Gigadalton sizes and can, in principle, comprise an unlimited number of different subunits (Ke et al., 2012; Reinhardt and Frenkel, 2014; Gerling et al., 2015; Wagenbauer et al., 2017). Notwithstanding these differences, a generic self-assembly process always includes three key steps: First, subunits must be made available, for example by gene expression, or rendered competent for binding, for example by nucleotide exchange (Alberts and Johnson, 2015; Chen et al., 2008; Whitelam, 2015) (‘activation’). Second, the formation of a structure must be initiated by a nucleation event (‘nucleation’). Due to cooperative or allosteric effects in binding, there might be a significant nucleation barrier (Chen et al., 2008; Jacobs and Frenkel, 2015; Sear, 2007; Lazaro and Hagan, 2016; Hagan and Elrad, 2010). Third, following nucleation, structures grow via aggregation of substructures (‘growth’). To avoid kinetic traps that may occur due to irreversibility or very slow disassembly of substructures (Hagan et al., 2011; Grant et al., 2011), structure nucleation must be significantly slower than growth (Zlotnick et al., 1999; Ke et al., 2012; Reinhardt and Frenkel, 2014; Wei et al., 2012; Jacobs et al., 2015; Hagan and Elrad, 2010). Physically speaking, there are no irreversible reactions. However, in the biological context, self-assembly describes the (relatively fast) formation of long-lasting, stable structures. Therefore, at least part of the assembly reactions are often considered to be irreversible on the time scale of the assembly process. In this manuscript we investigate, for a given target structure, whether the nature of the specific mechanism employed in order to slow down nucleation influences the yield of assembled product. To address this question, we examine a generic model that incorporates the key elements of self-assembly outlined above.

Model definition

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,,S which assemble into rings of size L. The cases S=1 and 1<SL (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., 2016Michaels 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 α, 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 1231 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 Lnuc, which in terms of classical nucleation theory corresponds to the structure size at which the free energy barrier has its maximum. For l<Lnuc attachment of monomers to existing structures and decay of structures (reversible binding) into monomers take place at size-dependent reaction rates μl and δl, respectively (Figure 1). Here, we focus on identical rates μl=μ and δl=δ. A discussion of the general case is given in Appendix 4. Above the nucleation size, polymers grow by attachment of monomers with reaction rate νμ 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).

Schematic description of the model.

(a) Rings of size L are assembled from S different particle species. N monomers of each species are initially in an inactive state (blue) and are activated at the same per-capita rate α. Once active (green), species with periodically consecutive index can bind to each other. Structures grow by attachment of single monomers. Below a critical nucleation size (Lnuc), structures of size l (light yellow) grow and decay into monomers at size-dependent rates μl and δl, respectively. Above the critical size, polymers (dark yellow) are stable and grow at size-independent rate ν until the ring is complete (the absorbing state; red). (b) Illustration of depletion traps. If nucleation is slow compared to growth, initiated structures are likely to be completed. Otherwise, many stable nuclei will form that cannot be completed before resources run out.

We investigate two scenarios for the control of nucleation speed, first separately and then in combination. For the ‘activation scenario’ we set μ=ν (all binding rates are equal) and control the assembly process by varying the activation rate α. For the ‘dimerization scenario’ all particles are inherently active (α) and we control the assembly process by varying the dimerization rate μ (we focus on Lnuc=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).

Results

Deterministic behavior in the macroscopic limit

First, we consider the macroscopic limit, N1, and investigate how assembly yield depends on the activation rate α (activation scenario) and the dimerization rate μ (dimerization scenario) for Lnuc=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, αth or μ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.

Deterministic behavior in the macroscopic limit N1.

(a, b) Yield for different particle numbers N (symbols) and ring sizes L (colors) for Lnuc=2. Decreasing either (a) the activation rate (‘activation scenario’: μ=ν ) or (b) the dimerization rate (‘dimerization scenario’: α) achieves perfect yield. The stochastic simulation results (symbols) average over 16 realizations and agree exactly with the integration of the chemical rate equations (lines). The threshold values (Equation 1) are indicated by the vertical dashed lines. Plotting yield against the dimensionless quantity α/(νC) causes the curves for different C to collapse into a single master curve (inset in a). For both scenarios there is no dependency on the number of species S in the deterministic limit. (c, d) Illustration showing how depletion traps are avoided by either slow activation (c) or slow dimerization (d). If the activation or the dimerization rate is small (large) compared to the growth rate, assembly paths leading to complete rings are favored (disfavored). The color scheme is the same as in Figure 1. (e) Deterministically, the size distribution of polymers behaves like a wave, and is shown for large and small activation rate for L=60, Lnuc=2, N=10000 and μ=ν=1. The distributions are obtained from a numerical integration of the deterministic mean-field dynamics, Equation 6, and are plotted for early, intermediate and final simulation times. The respective percentage of inactive monomers and complete rings is indicated by the symbols in the scale bar on the left or right.

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 Lnuc=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 α<αth or μ<μth, where

(1) αth=PανμνC(L-L)3andμth=Pμν(L-L)2

(see Appendix 3) with proportionality constants Pα=[πΓ(2/3)/Γ(7/6)]3/35.77 and Pμ=π2/24.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 α and μ=ν, respectively): While αth decreases cubically with the ring size L, μth does so only quadratically. Furthermore, the threshold activation rate α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 N1. Finally, if α is finite and μ<ν, the interplay between the two slow-nucleation scenarios may lead to enhanced yield. This is reflected by the factor ν/μ in αth, and we will come back to this point later when we discuss the stochastic effects.

In summary, for large particle numbers (N1), 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.

Stochastic effects in the case of reduced resources

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 ymax<1. Reducing the particle number N decreases this saturation value ymax until no finished structures are produced (ymax0). 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, ymax(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 Nyth 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>0th, on L can approximately be described by a power-law: N>0thLξ, with exponent ξ2.8 for L600. Consequently, for L=600 already more than 105 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 ymax 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 NSL, ymax 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 α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=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 S1, see Appendix 7. Moreover, increasing the nucleation size Lnuc, and with it the reversibility of binding, also increases ymax, 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.

Stochastic effects in the case of reduced resources.

(a, b) Yield of the fully heterogeneous system (S=L) for reduced number of particles (symbols) for L=60 and Lnuc=2 averaged over 1024 ensembles. In the activation scenario, at low activation rates the yield saturates at an imperfect value ymax, which decreases with the number of particles (a). This finding disagrees with the deterministic prediction (black line) of perfect yield for α0. In contrast, the dimerization scenario robustly exhibits the maximal yield of 1 for small N, in agreement with the deterministic prediction (black line) (b). (c) Diagram showing different regimes of ymax(N,L) in dependence of the particle number N and target size L (for the fully heterogeneous system S=L) as obtained from stochastic simulations in the limit α0. The minimal number of particles necessary to obtain a fixed yield increases in a strongly nonlinear way with the target size. The symbols along the line L=60 represent the saturation values of the yield curves in (a). (d) Dependence of ymax on the number of species S for fixed L=60 and fixed number of ring structures NS/L. Symbols indicate different values of the critical nucleation size Lnuc. The impact of stochastic effects strongly depends on the number of species under the constraint of a fixed total number of particles NS and fixed target size L. The homogeneous system is not subject to stochastic effects at all. Higher reversibility for larger Lnuc also mitigates stochastic effects.

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).

Cause and effect of stochasticity in the activation scenario.

(a) Illustration of the significance of stochastic effects when resources are sparse. Arrows indicate possible transitions and the probabilities in the depicted situation for sufficiently small activation rate α. For small α, the random order of activation alone determines the availability of monomers and therefore the order of binding. In the depicted situation, the complete structure is assembled only with probability 1/2. In all other cases, only fragments of the structure are assembled such that the final yield is decreased. (b) Polymer size distribution for the activation scenario (symbols) as obtained from stochastic simulations, in comparison with its deterministic prediction (lines) for S=L=100, N=1000 and Lnuc=2. Due to the enhanced number of nucleation events, the stochastic wave encompasses far more structures and moves more slowly. As a result, it does not quite reach the absorbing boundary.

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)).

Non-monotonic yield curves for a combination of slow dimerization and activation

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., μ<ν, α<). Figure 5 shows that, depending on the ratio μ/ν, the yield can become a non-monotonic function of α. In the regime where α is large, nucleation is dimerization-limited; therefore activation is irrelevant and the system behaves as in the dimerization scenario for α. Upon decreasing α 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 μ/ν dependence of αth, Equation 1), whilst the average number of active monomers is still high and fluctuations are negligible. Finally, a stochastic yield catastrophe occurs if α 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 μ/ν, 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).

Yield for a combination of slow dimerization and activation.

(a, b) Dependence of the yield of the fully heterogeneous system on the activation rate α for N=100 and different values of the dimerization rate (colors/symbols) for L=60 (a) and L=40 (b) (averaged over 1024 ensembles). For large activation rates the yield behaves deterministically (lines). In contrast, for small activation rates, stochastic effects (blue shading) lead to a decrease in yield. Depending on the parameters, the yield maximum is attained in either the deterministic, stochastic or intermediate regime. (c) Table summarizing the qualitative behavior of the yield (poor/intermediate/perfect) for a combination of dimerization and activation rates for both the deterministic and the stochastic limit. The columns correspond to low and high values of the dimerization rate, as indicated by the marker in the corresponding deterministic yield curve at the top of the column. Similarly, the rows correspond to low, intermediate and high activation rates. Arrows and colors indicate where and for which curve this behavior can be observed in (a) and (b). Deviations between the deterministic and stochastic limits are most prominent for low activation rates.

Robustness of the results to model modifications

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 L) bind at rate ν 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 N50 (Figure 6b).

Extended model including polymer-polymer binding.

(a) In the extended model, structures not only grow by monomer attachment but also by binding with another polymer (colored arrow). As before, binding only happens between periodically consecutive species with rate ν per binding site. So, the reaction rate for two polymers is identical to the one for monomer-polymer binding, ν. Furthermore, only polymers with combined length L can bind. All other processes and rules are the same as in the original model described in Figure 1. (b) The yield of the extended model as obtained from stochastic simulations is shown in dependence of the activation rate α for S=L=100, μ=ν=1, Lnuc=2 and different values of the number of particles per species, N (averaged over 1024 ensembles). The qualitative behavior is the same as for the original model. In particular, yield saturates (in the stochastic limit) at an imperfect value for slow activation rates. Note that for small particle numbers polymer-polymer binding results in an increase of the minimal yield (here for large activation rates). This is due to the fact that even in the case where a priori too many nucleation events happen, polymers can combine into final structures.

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.

Assembly of squares of size L×L from L different particle species.

(a) As in the ring models, there are N monomers of each species in the system. All particles are initially in an inactive state (blue) and are activated at the same per-capita rate α. Once active (green), species with neighboring position within the square (left/right, up/down) can bind to each other. Structures grow by attachment of single monomers until the square is complete (absorbing state). Depending on the number b of contacts between the monomer and the structure, the corresponding rate is bν. For simplicity, all polymers (yellow) are stable (Lnuc=2) and we do not consider polymer-polymer binding. (b) The yield of the square model as obtained from stochastic simulations is shown in dependence of the activation rate α for S=L=100, μ=ν=1 and different values of the number of particles per species, N (averaged over 256 ensembles). The qualitative behavior is the same as for the previous models: Whereas the yield is poor for large activation rates, it strongly increases below a threshold value and saturates (in the stochastic limit) at an imperfect value < 1 for small activation rates. The saturation value decreases with decreasing number of particles in the system.

Discussion

Our results show that different ways to slow down nucleation are indeed not equivalent, and that the explicit implementation is crucial for assembly efficiency. Susceptibility to stochastic effects is highly dependent on the specific scenario. Whereas systems for which dimerization limits nucleation are robust against stochastic effects, stochastic yield catastrophes can occur in heterogeneous systems when resource supply limits nucleation. The occurrence of stochastic yield catastrophes is not captured by the deterministic rate equations, for which the qualitative behavior of both scenarios is the same. Therefore, a stochastic description of the self-assembly process, which includes fluctuations in the availability of the different species, is required. The interplay between stochastic and deterministic dynamics can lead to a plethora of interesting behaviors. For example, the combination of slow activation and slow nucleation may result in a non-monotonic dependence of the yield on the activation rate. While deterministically, yield is always improved by decreasing the activation rate, stochastic fluctuations between species strongly suppress the yield for small activation rate by effectively enhancing the nucleation speed. This observation clearly demonstrates that a deterministically slow nucleation speed is not sufficient in order to obtain good yield in heterogeneous self-assembly. For example, a slow activation step does not necessarily result in few nucleation events although deterministically this behavior is expected. Thus, our results indicate that the slow nucleation principle has to be interpreted in terms of the stochastic framework and have important implications for yield optimization.

We showed that demographic noise can cause stochastic yield catastrophes in heterogeneous self-assembly. However, other types of noise, such as spatiotemporal fluctuations induced by diffusion, are also expected to trigger stochastic yield catastrophes. Hence, our results have broad implications for complex biological and artificial systems, which typically exhibit various sources of noise. We characterize conditions under which stochastic yield catastrophes occur, and demonstrate how they can be mitigated. These insights could usefully inform the design of experiments to circumvent yield catastrophes: In particular, while slow provision of constituents is a feasible strategy for experiments, it is highly susceptible to stochastic effects. On the other hand, irrespective of its robustness to stochastic effects, the experimental realization of the dimerization scenario relies on cooperative or allosteric effects in binding, and may therefore require more sophisticated design of the constituents (Sacanna et al., 2010; Zeravcic et al., 2017). Our theoretical analysis shows that stochasticity can be alleviated either by decreasing heterogeneity (presumably lowering realizable complexity) or by increasing reversibility (potentially requiring fine-tuning of bond strengths and reducing the stability of the assembly product). Alternative approaches to control stochasticity include the promotion of specific assembly paths (Murugan et al., 2015; Gartner, Graf and Frey, in preparation) and the control of fluctuations (Graf, Gartner and Frey, in preparation). One possibility to test these ideas and the ensuing control strategies could be via experiments based on DNA origami. Instead of building homogeneous ring structures as in Wagenbauer et al. (2017), one would have to design heterogeneous ring structures made from several different types of constituents with specified binding properties. By varying the opening angle of the ‘wedges’ (and thus the preferred number of building blocks in the ring) and/or the number of constituents, both the target structure size L as well as the heterogeneity of the target structure S could be controlled.

Moreover, the ideas presented in this manuscript are relevant for the understanding of intracellular self-assembly. In cells, provision of building blocks is typically a gradual process, as synthesis is either inherently slow or an explicit activation step, such as phosphorylation, is required. In addition, the constituents of the complex structures assembled in cells are usually present in small numbers and subject to diffusion. Hence, stochastic yield catastrophes would be expected to have devastating consequences for self-assembly, unless the relevant cellular processes use elaborate control mechanisms to circumvent stochastic effects. Further exploration of these control mechanisms should enhance the understanding of self-assembly processes in cells and help improve synthesis of complex nanostructures.

Materials and methods

All our simulation data was generated with either C++ or MATLAB. The source code is available at the eLife website.

Here we show the derivation of Equation 1 in the main text, giving the threshold values for the rate constants below which finite yield is obtained. The details can be found in Appendices 1–3.

Master equation and chemical rate equations

Request a detailed protocol

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 and species s in the system at time t by ns(t) with 2<L and 1sS. The species of a polymer is defined by the species of the respective monomer at its left end. Furthermore, n0s and n1s denote the number of inactive and active monomers of species s, respectively, and nL the number of complete rings. We signify the reaction rate for binding of a monomer to a polymer of size by ν. α denotes the activation rate and δ the decay rate of a polymer of size . By we indicate (ensemble) averages. The system governing the evolution of the first moments (the averages) of the {ns} is then given by:

(2a) ddtn0s=αn0s,
(2b) ddtn1s=αn0s=1L1ν(n1sns+1+n1sns)+=2Lnuc1k=s+1k=sδnk,
(2c) ddtn2s=ν1n1sn1s+1ν2n2sn1s+2ν2n2sn1s1δ2n2s1{2<Lnuc},
(2d) ddtns=ν1n1sn1+s1+ν1n1s+1n1sνnsn1s+νnsn1s1δns1{<Lnuc },
(2e) ddtnLs=νL1nL1sn1L+s1+νL1nL1s+1n1s.

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 α. 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 α, binding of active monomers to the left or to the right end of an existing structure of size at rate ν, and decay of below-critical polymers of size into monomers at rate δ (disassembly). Equations 2c and 2d describe the dynamics of dimers and larger polymers of size 3<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 𝟏{x<Lnuc} equals 1 if the condition x<Lnuc is satisfied and 0 otherwise. Note that a polymer of size 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).

Graphical illustration of Equations 2 and 6.

(a) Visualization of the gain and loss terms in the dynamics of the active monomers in Equation 2b. Gain of active monomers is due to activation of inactive monomers as well as decay of unstable polymers. Loss of active monomers is due to dimerization and attachment of monomers to larger polymers. (b) Visualization of the transitions between clusters of different sizes (without distinction of species). The first and second box represent the inactive and active monomers in the system, the subsequent boxes each represent the ensemble of polymers of a certain size. The arrows between the boxes show possible reactions and transitions with the reaction rates indicated accordingly. Each arrow starting from or leading to a box is associated with a corresponding loss or gain term on the right hand side of Equation 2 and Equation 6.

We simulated the Master equation underlying Equation 2 stochastically using Gillespie’s algorithm. For the following deterministic analysis, we neglect correlations between particle numbers {ns}, 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)

(3) nisnjk=nisnjks,k

Furthermore, for the expectation values it must hold

(4) ns=n1s

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

(5) c:=nsV,

we denote the concentration of any monomer or polymer species of size , 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)

(6a) ddtc0=-αc0,
(6b) ddtc1=αc0-2c1=1L-1νc+=2Lnuc-1lδc,
(6c) ddtc2=ν1c122ν2c1c2δ2c21{2<Lnuc},
(6d) ddtc=2ν1c1c12νc1cδc1{<L nuc},for 3<L,
(6e) ddtcL=2νL-1c1cL-1.

The rate constants ν 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 ν in Equation 6 can be interpreted in the usual units [litermol 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

(7) yield(t)=ScL(t)VSNL-1=cL(t)VLN.

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.

Effective description by an advection-diffusion equation

Request a detailed protocol

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 Lnuc= 2 and let ν1=μ and ν2=ν. Then, for the polymers with >2 we have

(8) tc=2νc1[c-1-c].

As a next step, we approximate the index {2,3,,L} indicating the length of the polymer as a continuous variable x[2,L] and define c(x=):=c. By A:=c1 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

(9) c(-1)=c()-xc()+12x2c(),

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

(10) tc(x)=-2νAxc(x)+νAx2c(x).

Equation 10 can be written in the form of a continuity equation tc(x)=-xJ(x) with flux J= 2νAc-νAxc. The flux at the left boundary x= 2 equals the influx of polymers due to dimerization of free monomers J(2,t)=μA2. This enforces a Robin boundary condition at x= 2

(11) 2νAc(2,t)-νAxc(2,t)=μA2.

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

(12) tA=αCe-αt-2μA2-2νA2Lc(x,t)𝑑x.

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= 2 from the boundary condition Equation 11, moves to the right with time-dependent velocity 2νA(t) due to the advection term, and broadens with a time-dependent diffusion coefficient ν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 dadv and the diffusively travelled distance ddiff exceeds the system size L-2

(13) dadv+ddiffL-2.

According to Equation 10, dadv=2ν0A(t)𝑑t and ddiff=2ν0A(t)𝑑t, giving as condition for the onset of finite yield

(14) 2ν0A(t)𝑑t=!14(1+4(L-2)-1)2L-L,

where the last approximation is valid for large L.

In order to obtain 0A(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):=2c(x,t)𝑑x (as long as yield is zero the upper boundary is irrelevant and we can consider L=). Equation 12 then reads

(15) ddtA=αCe-αt-2μA2-2νAB,

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

(16) ddtB=2tc(x,t)dx=2xJ(x,t)dx=J(,t)=0+J(2,t)=μA(t)2.

Measuring A and B in units of the initial monomer concentration C and time in units of (νC)-1 the equations are rewritten in dimensionless units as

(17a) ddtA=ωeωt2ηA22AB,
(17b) ddtB=ηA2,

where ω=ανC and η=μν. 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 ω and η

(18) 0Aω,η(t)𝑑t:=g(ω,η).

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

Dimerization scenario

Request a detailed protocol

The activation rate in the dimerization scenario is α, and instead of the term ωe-ωt in dA/dt, we set the initial condition A(0)=1 (and B(0)=0). Furthermore, η=μ/ν1 and we can neglect the term proportional to η in dA/dt. As a result,

dAdB=-2BηA.

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

(19) 2g(ω,η)=2π221η,

where for the evaluation of the integral we used the substitution ηA2dt=dB.

Activation scenario

Request a detailed protocol

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 ω1, we can again neglect the term proportional to η in dA/dt. This time, however, we have to keep the term ωe-ωt. As a next step, we assume that dA/dt is much smaller than the remaining terms on the right-hand side, ωe-ω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)=ωe-ωt/(2B(t)). Using it to solve dB/dt=ηA2 for B, and then to determine A, the totally travelled distance of the wave is deduced as

(20) 2g(ω,η)=232/3πΓ(2/3)6Γ(7/6)(ωη)-1/3.

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

(21) 2a(ηω)13LLα<αth:=PανμνC(LL)3
(22) or2bη12LLμ<μth:=Pμν(LL)2,

where a and b are numerical factors, and Pα= 8a35.77 and Pμ= 4b24.93. This verifies Equation 1 in the main text.

Appendix 1

Chemical reaction equations and the equivalence of models with different numbers of species

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.

Homogeneous structures

First, we consider the homogeneous model (S= 1). By c(t) we denote the concentration of complexes of length (2) at time t, c1(t) is the concentration of active monomers and c0(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 by ν. The model from the main text is recovered by setting ν:=μ if <Lnuc, and ν:=ν otherwise. The ensuing set of ordinary differential equations then reads:

(A1a) ddtc0=-αc0,
(A1b) ddtc1=αc0-2c1=1L-1νc+=2Lnuc-1lδc,
(A1c) ddtc2=ν1c122ν2c1c2δ2c21{2<Lnuc},
(A1d) ddtc=2ν1c1c12νc1cδc1{<L nuc},for 3<L,
(A1e) ddtcL=2νL-1c1cL-1.

The indicator function 𝟏{x<Lnuc} equals 1 if the condition x<Lnuc is satisfied and 0 otherwise. The first equation describes loss of inactive particles due to activation at rate α. It is uncoupled from the remainder of the equations and is solved by c0(t)=Ce-α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 α, binding of active monomers to existing structures at rate ν (polymerization), and decay of below-critical polymers into monomers at rate δ (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 2ν1c12 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 12n(n-1)n2/2 (if n is large) rather than n2 (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<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 tc2 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).

Heterogeneous structures

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,,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, cs with 0<L and 1sS denotes the concentration of a polymer of length and species s (c0s and c1s again denote inactive and active monomers of species s, respectively). For example, c45 denotes the concentration of polymers [5678] if S8, 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 cs with 3<L is given by

(A2) ddtcs=ν1c1sc1+s1+ν1c1s+1c1sνcsc1s+νcsc1s1δcs1{<Lnuc}.

The terms on the right-hand side account for the influx due to binding of the respective polymers of length -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 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 <Lnuc (last term). For the dynamics of the dimers, however, there is only one gain term arising from dimerization:

(A3) ddtc2s=ν1c1sc1s+1ν2c2sc1s+2ν2c2sc1s1δ2c2s1{2<Lnuc}.

Equivalently, for the active monomers we find:

(A4) ddtc1s=αCe-αt-c1s=1L-1ν(cs+1+cs-)+=2Lnuc-1k=s+1-k=sδck.

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

(A5) cs(t)=ck(t),foranys,kSatanytimet.

In other words, the upper index is irrelevant and can also be discarded. The variable c then denotes the concentration of any one polymer species of length . 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 3 in Equation A2. Applying cs(t)=c(t) to Equation A2 yields for the dynamics of the concentration of an arbitrary polymer species of size :

ddtc=ν1c1c1+ν1c1c1νcc1νcc1δc1{<Lnuc}.=2ν1c1c12νcc1δc1{<Lnuc},

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:

(A6) yield(t)=ScL(t)SNL-1=cL(t)LN.

Appendix 2

Effective description of the evolution of the polymer size distribution as an advection-diffusion equation

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 {2,3,,L} as a continuous variable that we rename x[2,L]. With such a continuous description in view we write c(x=):=c to denote the concentration of polymers of size .

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 Lnuc= 2 and let ν1=μ and ν2=ν. Generalizations to Lnuc>2 can be done in a similar way. Then, for the polymers with 3 we have:

(A7) tc()=2νA[c(-1)-c()].

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

(A8) c(-1)=c()-xc()+12x2c(),

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

(A9) tc(x)=-2νAxc(x)+νAx2c(x).

Equation A9 can be written in the form of a continuity equation tc(x)=-xJ(x) with flux J= 2νAc-νAxc. The flux at the left boundary, x= 2, equals the influx of polymers due to dimerization of free monomers, J(2,t)=μA2. This enforces a Robin boundary condition at x= 2,

(A10) 2νAc(2,t)-νAxc(2,t)=μA2.

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

(A11) tA=αCe-αt-2μA2-2νA2Lc(x,t)𝑑x.

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= 2 from the boundary condition, Equation A10, moves to the right with time dependent velocity 2νA(t) due to the advection term, and broadens with a time-dependent diffusion coefficient νA(t). The concentration of active particles A determines both the influx of dimers at x= 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 ddtc(2,t)=μA22νAc(2,t) (see Equation A1c) by ’variation of the constants’ as

(A12) c(2,t)=0tμA(t~)2exp[2t~tνA(t)dt]dt~.

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

(A13) tc~(x)=-2νAxc~(x).

is solved by

(A14) cadvec(x,t)=c(2,τ(x,t)).

Here, τ(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 τ(x,t). This ansatz solves the PDE (Equation A13) if and only if τ(x,t) satisfies

(A15) τ(x,t)=A~-1(A~(t)-x-22ν)

with A~ being an arbitrary integral of A such that tA~(t)=A(t) and A~-1 denoting its inverse. More easily, we find this form of τ by requiring that the integral over the velocity from time τ to t equals the travelled distance x-2:

(A16) τt2νA(t)𝑑t=x-2.

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

(A17) k(x,y,t)=(4πτ(y,t)tD(t))-1/2exp(-x24τ(y,t)tD(t)),

with the time dependent diffusion constant D(t)=ν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 τ(y,t), it diffused for the time t-τ(y,t). The complete expression for c(x,t) is then obtained as the convolution of cadvec(x,t) (Equation A14), that is obtained from Equation A12 and Equation A15, and the diffusion kernel k(x,y,t) (Equation A17):

(A18) c(x,t)=cadvec(s,t)k(x-s,s,t)𝑑s=c(2,τ(s,t))k(x-s,s,t)𝑑s.

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 μν, 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 μν, the motion of the wave can be sustained until it reaches the absorbing boundary.

Appendix 3

Threshold values for the activation and dimerization rate

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 dadv and the diffusively travelled distance ddiff exceeds the system size L-2:

(A19) dadv+ddiffL-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 τ=0 and t=:

(A20) dadv=02νA(t)𝑑t.

The diffusively travelled distance is approximately given by the standard deviation of the Gaussian diffusion kernel, Equation A17, again with τ=0 and t=,

(A21) ddiff=2ν0A(t)𝑑t.

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

(A22) 2ν0A(t)𝑑t+2ν0A(t)𝑑t=L-2.

Substituting y=2νA and requiring that y is positive, we solve the quadratic equation and find that Equation A22 is equivalent to

(A23) 2ν0A(t)𝑑t=y2=14(1+4(L-2)-1)2L-L,

where the last approximation is valid for large L.

We determine the threshold values for the activation rate α and the dimerization rate μ 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 0A(t)𝑑t on α and μ.

The concentration c(x,t) appears in Equation A11 only in terms of an integral 2Lc(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= 2. As long as yield is zero we can therefore equivalently consider the limit L. We denote the total number of polymers in Equation A11 by B(t):=c(x,t)𝑑x for which the dynamics is determined from the boundary condition, Equation A10:

(A24) ddtB=2tc(x,t)𝑑x=2-xJ(x,t)dx=-J(,t)=0+J(2,t)=μA(t)2.

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:

(A25a) ddtA=αCeαt2μA22νAB
(A25b) ddtB=μA2.

The dynamics of A and B is equivalent to a two-state activator-inhibitor system, where A dimerizes into B at rate μ, and B degrades (inhibits) A at rate 2ν. 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 (νC)-1:

(A26a) ddtA=ωeωt2ηA22AB,
(A26b) ddtB=ηA2,

with the remaining dimensionless parameters ω=ανC and η=μν. We are interested in the integral over A(t) as a function of ω and η,

(A27) 0Aω,η(t)𝑑t:=g(ω,η),

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

Analysis of the dimerization scenario

The dimerization scenario is characterized by fast activation αCν and slow dimerization μν. For the dimensionless parameters these assumptions translate to η1 and ηω. Because for small η1 nucleation is much slower than growth we neglect the dimerization term in Equation A26a against the growth term. Furthermore, because ηω 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

(A28a) ddtA=2AB,
(A28b) ddtB=ηA2,

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):

(A29) dAdB=-2ηBA,

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

(A30) A(B)=1-2ηB2.

Because the number of active monomers A(t) must vanish for t, the final value of B is

(A31) B:=B(t=)=η2.

Thereby, we calculate the function g(η) via variable substitution dt=dBηA2:

(A32) g(η)=0A(t)𝑑t=0BA(B)dBηA(B)2=1η0BdB1-2ηB2=π22η-12.

So, the dependence of the travelled distance of the wave on η obeys a power law with exponent 12, confirming the previous result (Morozov et al., 2009). For the coefficient we find π221.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

(A33) B(t)=η2tanh(2ηt).

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

(A34) A(t)=1cosh(2ηt).

Analysis of the activation scenario

In the activation scenario, αCν, such that ω1 and ωη. As we know already that decreasing ω 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 τ:=1-e-ωt such that τ[0,1] and writing a(τ)=a(1-e-ωt):=A(t) and b(τ)=b(1-e-ωt):=B(t) the system in Equation A26 becomes:

(A35a) ddτa=12ω(1τ)ab,
(A35b) ddτb=ηω(1τ)a2,

with the initial condition a(0)=b(0)=0. The function g(ω,η) transforms as

(A36) g(ω,η)=0A(t)𝑑t=1ω01a(τ)1-τ𝑑τ.

In the following we derive the asymptotic solution for a(τ) in the limit of small ω in order to evaluate the integral in Equation A36. In the limit τ1 (t) both a(τ) and ddτa(τ) will become small whereas b(τ) increases monotonically. The reaction term in Equation A35a is furthermore weighted by a factor 1ω which will become large if ω1. We therefore postulate that for sufficiently large τ the derivative ddτa(τ) 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 ddτa in (Equation A35a) reduces the equation to an algebraic equation and we find

(A37) a=ω(1-τ)2b.

Using this result in Equation A35b we can solve for b by separation of variables and subsequent integration:

(A38) b(τ)=(ωη)13(34τ-38τ2)13.

From Equation A37 we immediately obtain a(τ):

(A39) a(τ)=ω23η131-τ(6τ-3τ2)13:=ω23η13h(τ),

where by h(τ) we denote the part of the solution that depends only on τ. Hence, we find that a and hence also ddτa scale like ω23, and will thus become small if ω1 and τ is large enough. Therefore the solution is consistent and justifies the approximation in which we neglected the derivative term in the limit of small ω and sufficiently large τ.

Note that consistency of the solution with the approximation is a sufficient criterion for the validity of the approximation: We can solve the system for A and B in Equation A35 iteratively by defining

ddτai1=12ω(1τ)aibi,
ddτbi=ηω(1τ)ai2.

Assuming that for i, ai and bi converge to the correct solutions a(τ) and b(τ) when starting with a0=0, we obtain a1 and b1 as given by Equation A39 and Equation A38 and can iteratively refine the approximation. The next iteration step then reads: ddτa1=1-2ω(1-τ)a2b2. As a1ω23 we know that the left-hand side will be small and a1 and b1 solve the system if the left-hand side equals 0. Writing a2=a1+a~2 and b2=b1+b~2 this gives:

(A40) ddτa1=12ω(1τ)(a1+a~2)(b1+b~2)2ω(1τ)(a1b~2+b1a~2).

From dimensional analysis it follows that the correction terms a~2 and b~2 must scale like a~2ω43 and b~2ω and are hence much smaller than the first order approximations a1 and b1. Higher order corrections will give even smaller contributions showing that if ddτa11, a1 is indeed a very good approximation.

In the limit τ0, however, the expression for a(τ) in Equation A39 diverges and consistency is violated. Hence, the obtained solution is valid only for sufficiently large τ.

We fix some small ϵ>0 such that the approximation can be assumed to be sufficiently good if ddta<ϵ. Furthermore, we define τϵ such that ddτa<ϵ for all τ>τϵ. Using Equation A39 we can write this as ddτh<ϵη13/ω23 for all τ>τϵ, where the left-hand side, ddτh, depends only on τ. Hence, by decreasing ω we can make τϵ arbitrarily small: limω0τϵ=0. In order to calculate g(ω,η) the integral in Equation A36 can be separated in a domain where the approximation a(τ) is accurate and a domain where the correct solution a~(τ) deviates strongly from a(τ):

(A41) g(ω,η)=1ω0τϵa~(τ)1-τ𝑑τ+1ωτϵ1a(τ)1-τ𝑑τ.

We see from Equation A35a that ddτa~=1 describes an upper bound to a~ showing that a~(τ)τ. Therefore we can bound the contribution of the first integral as 0τϵa~(τ)1-τ𝑑τ0τϵτ1-τϵ𝑑τ=12τϵ21-τϵ. Because this upper bound for the integral goes to 0 if ω and hence τϵ become small the first integral will become negligible against the second one. Asymptotically, we therefore only need to consider the second integral with the solution for a(τ) as given by Equation A39:

(A42) g(ω,η)=(ωη)1301(6τ3τ2)13dτ=(ωη)1303dz6z131z3==323π Γ(23)6 Γ(76)(ωη)130.8969(ωη)13,

where we used the substitution τ=11z/3 and Γ(x) is the (Euler) Gamma function. So, in the limit of small ω, g scales with ω and η with identical exponent 13. This contrasts the dimerization scenario where g as well as A and B depend only on η and are independent of ω (cf. Equation A32, A33 and A34).

Numerical analysis and the threshold values for the rate constants

In order to confirm the results of the last two paragraphs and to see how g(ω,η) behaves in the intermediate regime where ω and η are of the same order of magnitude we also investigate the function g(ω,η) numerically. For that purpose we numerically integrate the ODE-system for A(t) and B(t) in Equation A26 for different values of ω and η with a semi-implicit method. Subsequently, we integrate the solution A(t) using an adaptive recursive Simpson’s rule. Plotting g in dependence of ω for fixed η on a double-logarithmic scale reveals a rather simple bipartite form of g, see Appendix 3—figure 1a:

(A43) g(ω,η)={g1(η)ω-13ω1g2(η)ω1.
Appendix 3—figure 1
Fit of g(ω,η) on log-log scale.

The function g(ω,η)=0Aω,η(t)𝑑t describes (half) the travelled distance of the profile of the polymer size distribution in dependence of ω=ανC and η=μν. Marker points show solutions for g(ω,η) as obtained numerically from integration of Equation A26. Red lines are linear fits on log-log scale. In (a) we plot g(ω,η) for fixed η (here exemplarily for η= 0.01) over 25 orders of magnitude in ω and find a markedly bipartite behavior: For small ω the dependence on ω is perfectly matched by a power law with exponent 13 and η-dependent coefficient g1(η), whereas for large ω it is a constant g2(η). (b) Plotting g2(η)=g(ω=,η) in dependence of η reveals again strictly bipartite behavior. Here, however, only the branch for small η1 is relevant. With the coefficient g1(η) that can be determined in a similar way this leads to the final form of g(ω,η) as given by Equation A46.

The transition between these two regimes is rather sharp so that g is best described in a piecewise fashion

(A44) g(ω,η)=max(g1(η)ω-13,g2(η)).

Next, we plot the coefficients g1(η) and g2(η) against η. Here we find that g1(η)=aη-13 with a=const0.90 and g2(η) is again bipartite with a sharp kink in between (Appendix 3—figure 1b):

(A45) g2(η)=min(bη-12,bη-0.85),

where b1.11 and b1.37. The transition between both regimes is at η1.82. The second regime is not relevant for self-assembly since it refers to both large ω and large η, hence the travelled distance 2g is too small to give finite yield in this regime. Therefore, we discard the second regime and obtain as final result

(A46) g(ω,η)=max(a(ηω)-13,bη-12),

with a0.90 and b1.11. This confirms perfectly the exponents as well as the coefficients found in the last two paragraphs. It is, however, surprising that there is such a sharp transition between both regimes, which allows to define g(ω,η) in a piecewise fashion. This behavior must be the result of a series of lower oder terms in g(ω,η) which are unimportant in the limits ωη and ηω but cause the sharp transition when ω and η are of the same order of magnitude.

Finally, we return to our original task of finding the threshold values of the activation and dimerization rate for the onset of yield. Using our result for g(ω,η) in Equation A23 we find as necessary and sufficient condition to obtain finite yield in the deterministic system:

(A47) 2max(a(ηω)-13,bη-12)L-L.

Alternatively, we can state this result as two separate conditions out of which at least one must be fulfilled to obtain finite yield:

(A48) 2a(ηω)13LLα<αth:=PανμνC(LL)3
(A49) or2bη12LLμ<μth:=Pμν(LL)2

where Pα= 8a35.77 and Pμ= 4b24.93. This verifies Equation 1 in the main text.

Appendix 4

Impact of the implementation of sub-nucleation reactions

In the main text we focused our discussion on irreversible binding Lnuc= 2. In this section we investigate the effect of different implementations of the sub-nucleation reactions.

In general, perfect yield is trivially achieved if the complete ring is the only stable structure. However, yield can be maximal already for smaller nucleation sizes Lnuc depending on the explicit decay rate δ. In the deterministic limit without the dimerization and activation mechanisms (μ=ν, α ) a rapid transition from zero yield to perfect yield occurs in dependence of the critical nucleation size (see Appendix 4—figure 1). The threshold value in this case is approximately half the ring size and is weakly affected by the decay rate δ. In order to obtain finite yield for small nucleation sizes, an extremely high decay rate would be necessary. Hence, maximizing the yield solely by increasing the nucleation size is not very feasible.

Appendix 4—figure 1
Yield maximization due to increased nucleation size.

Without activation and dimerization mechanism (α,μ=ν) the yield can still be optimized by increasing the critical nucleation size Lnuc. However, a significant improvement is only achieved for critical sizes larger than half the ring size. Above, a rapid transition to perfect yield takes place. Below no effect is observed at all. Increasing δ shifts the onset of yield to slightly smaller critical nucleation sizes. Other parameters: L= 60, N= 10000.

In our model, the subcritical reaction rates μi may take different values. Here, we want to restrict our discussion to two scenarios. First, all rates have an identical value μi=μ and second, the rates increase linearly up to the super-nucleation reaction rate: μi=μ+(ν-μ)i-1Lnuc-1.

In the deterministic limit, both implementations show the same qualitative behavior as the dimerization mechanism with Lnuc= 2 in the main text (see Appendix 4—figure 2). The only relevant aspect for the final yield is the extend to which nucleation is slowed down in total. In the constant scenario all reaction steps contribute equally. As a results there is a strong dependence on the number of such reaction steps, that is on the critical nucleation size. If however, the reaction rates increase linearly with the size of the polymers, the dimerzation rate dominates. Only in the case μν finite yield is observed at all. In this limit the dimerization rate is much smaller than the subsequent growth rates. The explicit form of the different μi is not of major importance for the yield. The total slowdown of nucleation is the central feature. Structure decay does not play any role for intermediate nucleation sizes.

Appendix 4—figure 2
Yield for the dimerization mechanism (α) with different nucleation sizes (colors).

(a) If all sub-nucleation growth rates are identical (μi=μ) increasing the nucleation size increases the threshold value μth. The slow down of nucleation due to the individual sub-nucleation steps in total determines the yield. (b) If the sub-nucleation growth rates increase linearly (μi=μ+(ν-μ)i-1Lnuc-1) no dependence on the nucleation size is observed. The dimerization rate μ1=μ (which is the most limiting step) dominates entirely. Other parameters: L= 60, N= 10000, δ= 1.

The last question we want to address is how the combination of activation and dimerization mechanism and the corresponding non-monotonic behavior is affected by the nucleation size. Again, we compare constant sub-nucleation growth with a linearly increasing growth rate (see Appendix 4—figure 3). In the deterministic regime both implementations behave qualitatively similar as the dimerization mechanism discussed in the main text. However, in both cases the stochastic yield catastrophe is less pronounced. For the constant growth rates a saturation of the maximal yield is observed for sufficiently low μ. If the profile is linear this effect is weaker as compared to the constant case and a dependency on the explicit value of μ is still observed. The saturation value is not reached for these reactions rates.

Appendix 4—figure 3
Combined mechanisms for different nucleation sizes (symbols) and dimerization rates (color).

(a) If the sub-nucleation growth rates are identical (μi=μ) the stochastic yield catastrophe is weakened but still has a drastic impact. The qualitative behavior remains unchanged. (b) For a linearly increasing sub-nucleation growth rate (μi=μ+(ν-μ)i-1Lnuc-1) in the deterministic regime no changes are observed at all. The effect of the stochastic yield catastrophe is less pronounced. This improvement is mainly caused by structure decay which mitigates stochastic fluctuations. However, a slight dependency of the saturation value on the rate μ is observed. Other parameters: L= 60, S=L, N= 100, δ= 0.1.

Taking all our results for the sub-nucleation behavior together we draw the following conclusions: First, structure decay by itself it not very efficient in order to maximize yield. Second, the explicit choice of the sub-nucleation rates is of minor importance for the qualitative behavior. The system behaves similarly to the case Lnuc= 2. Third, larger nucleation sizes mitigate the stochastic yield catastrophe in general.

Appendix 5

Time evolution of the yield in the activation and dimerization scenario

In the main text we focus on the final yield, which represents the maximal yield that can be obtained in the assembly reaction for t. Here, we briefly discuss the temporal evolution of the yield in the two scenarios. Appendix 5—figure 1 shows the yield as a function of time for the dimerization scenario (blue) and the activation scenario (red) for the corresponding parameters indicated in the plot. Drawn lines show the evolution of the yield in the stochastic simulation whereas dashed lines represent its deterministic evolution obtained by integrating the corresponding mean-field rate equations (only shown for the activation scenario). In both scenarios, yield production sets in after a short lag time (Hagan and Elrad, 2010). The emergence of a lag time can be understood in terms of the interpretation of the assembly process as the progression of a travelling wave (see Sec. B). The travelling wave thereby describes the polymer size distribution and the time that is needed for the wave to reach the absorbing boundary equals the lag time for yield production observed in Appendix 5—figure 1. After the lag time, the yield increases very abruptly in the dimerization scenario and a bit more continually in the activation scenario. Since monomers are provided gradually in the activation scenario, the emerging wave is flatter and extends over a larger range (in polymer size space) as compared to the dimerization scenario. Consequently, yield production is more gradual in the activation scenario than in the dimerization scenario. For the same reason, the dimerization scenario is generally ‘faster’ or more time efficient than the activation scenario. For a detailed analysis of the time efficiency of these and other self-assembly scenarios we refer the reader to our manuscript in preparation (Gartner, Graf and Frey, in preparation).

In all depicted situations, the yield increases monotonically with time. This is, of course, generally true since the completed ring structures define an absorbing state in our system. The final yield, which is indicated in the right bar, therefore represents the upper limit for the yield that can be achieved in the assembly reaction. Appendix 5—figure 1 shows that the temporal yield curves initially are rather steep and quickly reach a value that lies within 10% of the final yield (‘quickly’ thereby refers to the respective time scale), before the curves flatten and increase more slowly. This underlines that the final yield is a meaningful observable that not only describes the upper limit for the yield but also approximates the typical yield of the assembly reaction under appropriate time constraints that are not too restrictive (on the time scale set by the respective lag time).

Appendix 5—figure 1
Time evolution of the yield in the activation and dimerization scenario.

The time dependence of the yield is depicted for a dimerization scenario (blue) with μ=5×10-4 and N=100 and for two activation scenarios (red) with α=0.1 and N=2×102 and N=104, respectively, for target structures of size L=20. Drawn lines show the time evolution of the stochastic systems while dashed lines describe the time evolution in the corresponding deterministic systems (where the final yield may be higher in the activation scenario). In all cases the yield increases monotonically with time. The final yield, that is indicated in the right bar, represents the upper limit of the yield at any time. Yield production in the activation scenario is generally more gradual than in the dimerization scenario. Therefore, the dimerization scenario is, in general, more time efficient than the activation scenario.

Appendix 6

Standard deviation of the yield

In the main text, the analysis focuses on the average yield. A priori it is, however, not apparent that this average quantity is informative, in particular due to the strong effect of stochasticity in the system. Here, we thus take a step forward to complement this picture by additionally considering a simple measure for the fluctuations of the yield, its standard deviation. Appendix 6—figure 1 is an extension of Figure 3a in the main text, showing the dependence of the average yield and its sample standard deviation on the activation rate. Since yield is always positive, the standard deviation of the yield has to be small if the average yield is close to 0 (N=500 in Appendix 6—figure 1). The same holds true for average yield close to 1 as the yield is bounded by one from above (N=5000 in Appendix 6—figure 1). For intermediate values of the average yield, the standard deviation is highest but still small compared to the average yield (N=1000 in Appendix 6—figure 1). The average yield is, thus, meaningful. Naturally the ratio of the standard deviation compared to the average yield also depends on the number of particles per species N and on the number of species S. Generally speaking, for higher N and S, this ratio decreases (see Appendix 7—figure 1 for the dependency on S).

Appendix 6—figure 1
Average yield and its sample standard deviation.

For average yield close to 0 or close to 1, the standard deviation has to be small due to the boundedness of the yield to the interval [0, 1]. For intermediate values, the standard deviation is highest. Its value is, however, still considerably smaller than the average yield. The parameters are L= 60, S=L, μ=ν=1 and different particle numbers N (colors/symbols). To obtain the average yield, the yield has been averaged over 1000 simulations. The standard deviation corresponds to the unbiased sample standard deviation.

Appendix 7

Influence of the heterogeneity of the target structure for fixed number of particles per species

Figure 3d in the main text shows how the maximal yield ymax depends on the number of species S if the ring size L and the number of possible ring structures NS/L is fixed. This comparison for fixed NS is motivated by the question which role the heterogeneity of a structure plays for assembly efficiency if a certain number of structures should be realized. Figure 3d illustrates that a higher number of species S (more heterogeneous structures) leads to a lower maximally possible yield, suggesting that it is beneficial to build structures with as few different species as possible. However, this situation does not correspond to the deterministically equivalent case of fixed number of particles per species N (note, though, that in the deterministic case the maximally possible yield is always 1, namely for α0). Instead, for higher number of species S, the number of particles per species N1/S decreases. How does the heterogeneity of the structures S alter the maximally possible yield if L and N (instead of L and NS) are fixed? Appendix 7—figure 1 shows how the maximal yield ymax and its standard deviation (obtained as average yield and sample standard deviation for α=10-8 when the yield has well saturated and the dynamics (except for the timescale) get independent of the exact value of the rate-limiting activation rate) depend on the number of species S. For homogeneous structures S=1 yield is always perfect since in this case there can be no fluctuations between species. As a result, the average yield is 1 and the standard deviation is 0. For increasing S, the average yield decreases until it levels off for S1. This behavior indicates that indeed the decreasing number of particles per species N for larger S is essential for the decrease of the maximal yield with S in Figure 3d. As mentioned above, the standard deviation is largest for small S>1 and decreases with S.

Appendix 7—figure 1
Influence of the heterogeneity of the target structure on the yield for fixed number of particles per species N.

The maximal yield and its standard deviation (obtained as average yield and sample standard deviation for α=10-8) are plotted against the number of species S making up the structure of size L=60. The number of particles per species N=1000 is fixed. Yield drops from a perfect value of 1 for S=1 to a smaller value and levels off for S1. The standard deviation is largest for small S (except for S=1 where the yield is always perfect) and decreases with increasing number of species.

Appendix 8

Dependence of the maximal yield ymax in the activation scenario on N and L

Figure 3c in the main text characterizes the dependence of the maximal yield ymax in the activation scenario as a ‘phase diagram’ distinguishing different regimes of ymax in dependence of the particle number N and target size L. Supplementing this figure in the main text, Appendix 8—figure 1 shows the maximum yield that is obtained in the activation scenario in the limit α0 for fixed L in dependence of N (Appendix 8—figure 1a) as well as for fixed N in dependence of L (Appendix 8—figure 1b). For larger particle number N, the maximal yield exhibits a transition from 0 to 1 over roughly three orders of magnitude. Increasing L shifts the transition to larger N. The threshold particle number where the transition starts is characterised by Nth>0(L) (see main text). Approximately, for L600, we find Nth>0(L)L2.8 (cf. main text, Figure 3c). Similarly, decreasing the target size L for fixed N, the maximal yield exhibits a transition from 0 to 1 over roughly one order of magnitude in L. The corresponding threshold value Lth>0 as a function of N is obtained as the inverse function of Nth>0(L). Hence, at least for N105, approximately it holds Lth>0(N)N0.36. Since ymax is largely independent of the number of species S for fixed N and L (see Appendix 7), the maximal yield in the activation scenario (for Lnuc=2) can be fully characterized as a function ymax(N,L) of N and L. Hence, ymax can roughly be expressed in terms of the threshold particle number Nth>0(L) as

(A50) ymax(N,L){1 if N>103Nth>0(L)<1 if Nth>0(L)<N<103Nth>0(L)=0 if N<Nth>0(L)
Appendix 8—figure 1
Dependence of the maximal yield y𝐦𝐚𝐱 in the activation scenario on N and L.

For each data point, ymax was determined as the average yield of 100 independent stochastic simulations of the activation scenario with α=10-12. (a) Variation of the particle number N for different target sizes L. The maximal yield increases from 0 to 1 over roughly three order of magnitude in N. The onset of the transition depends on L. (b) Variation of the target size L for different particle numbers N. Increasing the target size L with N being fixed causes the maximal yield to drop to 0. The transition from 1 to 0 spans roughly one order of magnitude in L and its position is determined by N.

As can be seen from Figure 3c in the main text, the transition line between zero and nonzero yield slightly flattens with increasing L. Hence, the power law Nth>0(L)L2.8 (and similarly for Lth>0) only holds approximately and for a restricted range in L and N. The asymptotic behavior of Nth>0 in the limit L remains elusive.

References

  1. 1
    Molecular Biology of the Cell
    1. B Alberts
    2. A Johnson
    (2015)
    New York: Garland Science.
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
    A Kinetic View of Statistical Physics
    1. PL Krapivsky
    2. S Redner
    3. E Ben-Naim
    (2010)
    Cambridge University Press.
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42

Decision letter

  1. Frank Jülicher
    Reviewing Editor; Max Planck Institute for the Physics of Complex Systems, Germany
  2. Naama Barkai
    Senior Editor; Weizmann Institute of Science, Israel
  3. Pablo Sartori
    Reviewer; The Rockefeller University, United States

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

The authors study the role of fluctuations for the yield of self-assembly of heterogenous molecular complexes in specific arrangements. The approach captures the combination of subunit activation, nucleation of complex formation and complex assembly, starting from a set of different subunits in solution. This study shows in particular that when nucleation is limited by activation of subunits, yield can drop catastrophically due to stochasticity of activation and nucleation. This occurs in regimes where yield would be high when fluctuations are neglected. This work provides deep new insights in the role of stochasticity for the reliable self-assembly of molecular complexes and a framework for the study of molecular assembly in cells.

Decision letter after peer review:

Thank you for submitting your article "Stochastic yield catastrophes and robustness in self-assembly" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Naama Barkai as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Pablo Sartori (Reviewer #2).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

The authors study the role of fluctuations for the yield of self-assembly of heterogenous molecular structures. An elegant and simple model is introduced which captures fluctuations of both activation of monomers and stochastic nucleation. The paper first recapitulates known results that slow nucleation is key to high yield. It is shown that in the deterministic or mean field limit slow nucleation by activation or by cooperative assembly (dimerization) play a similar role. However, in the presence of fluctuations, both scenarios can show very different behaviors. In particular the yield can drop in the activation scenario due to fluctuations which the authors term stochastic yield catastrophe. Interesting is the general case where activation and dimerization occur together. The authors show that in this case there can be an optimal activation rate at which yield is maximized. The work provides fundamental new insights in the self-assembly of macromolecular structures in cells.

However, the referees also raised several points that need to be addressed in a revision. In particular the depth and relevance of some of the results is not fully clear. The authors should address the specific points given below.

Essential revisions:

The relevance and strength of the results remains partly unclear. Several points need clarification:

1) The authors use the term "Stochastic yield catastrophe" to describe the drop of yield for increased fluctuations. In contrast to the catastrophic drop of yield in the deterministic case, it is unclear whether for these stochastic effects the term stochastic yield catastrophe is meaningful and appropriate as the drop is gradual and seems soft rather than representing a transition.

2) In Figure 3D the effect of increasing the number of species for a fixed size L is studied. It is shown that more species results in lower maximal yield. This could simply be an effect of the constraint NS = constant, which results in a decreasing number of components per species, N. The authors should check whether, when keeping N = constant, changing S still has an effect on the yield.

3) Still in Figure 3, the characterization of the "stochastic yield catastrophe" seems incomplete in view of its relevance to the authors. It seems that there is a certain dependence ymax(N,L,S) (forgetting about Lnuc for now). Immediate questions that are not addressed are: how does ymax scale with N, for fixed L and S? How does it then scale with L? And how with the heterogeneity, L/S? One step in this direction is Figure 3D, in which the authors show ymax(S) for several values of L, but variable N (why then the label N = 1000 in the panel?). More such plots should be presented to clarify the key properties of the system.

4) There are additional questions which relate to the relevance of the results in realistic situations. The model is elegant but very much simplified. For example, growth is fully reversible which may not be the case in many systems of macromolecular assembly. It is unclear how general and robust the results of the simplified model are if the simplifications are relaxed.

5) The time needed for assembly is not discussed. In a biological context it is essential that assembly occurs sufficiently fast. High yield in the long term alone is not enough for assembly to happen. I understand that only looking at yield at long times is an elegant simplification of the problem, but it remains unclear how useful these insights are for real situations.

https://doi.org/10.7554/eLife.51020.sa1

Author response

Essential revisions:

The relevance and strength of the results remains partly unclear. Several points need clarification:

1) The authors use the term "Stochastic yield catastrophe" to describe the drop of yield for increased fluctuations. In contrast to the catastrophic drop of yield in the deterministic case, it is unclear whether for these stochastic effects the term stochastic yield catastrophe is meaningful and appropriate as the drop is gradual and seems soft rather than representing a transition.

The question alerted us to the fact that we have to be more explicit about the reason why we use the term’ catastrophe’. The stochastic yield catastrophe is not the change in yield with respect to the parameter α. The term is intended to refer to situations in which we deterministically expect a perfect yield but due to stochastic effects we may end up with zero yield instead. This drop in yield of up to 100% compared to the deterministically expected value is what we refer to as a ’stochastic catastrophe’. We have revised the definition of the term in the subsection “Stochastic effects in the case of reduced resources” of the main text to make this point clearer. Furthermore, we extended our discussion of the nonlinearities in the stochastic effects with respect to different model parameters (especially the ring size L) to give a more complete picture of their significance as well. This directly links to answer 3 in which those aspects are discussed in detail.

– We have adapted the introduction of the term “Stochastic yield catastrophe” in subsection “Stochastic effects in the case of reduced resources”.

2) In Figure 3D the effect of increasing the number of species for a fixed size L is studied. It is shown that more species results in lower maximal yield. This could simply be an effect of the constraint NS = constant, which results in a decreasing number of components per species, N. The authors should check whether, when keeping N = constant, changing S still has an effect on the yield.

Indeed, if N (instead of NS) is fixed in Figure 3D, the maximal yield becomes independent of S for large S >> 1. However, this saturation value of the yield differs from the yield for the homogeneous system S = 1, where the maximal yield is always 1 as the system behaves deterministically. So, the yield (sharply) drops from a perfect value of 1 (for S = 1) to a non-perfect value for S > 1 and then levels off. Intriguingly, for some cases, one even observes non-monotonic behavior of the yield with respect to the number of species: Yield decreases from 1 (for S = 1) to a lower value and then increases again before it saturates for S >> 1. This rather non-intuitive behavior needs a more thorough discussion which we provide in detail in a follow-up manuscript (in preparation). In short, we believe that similar to the deterministic description also the stochastic limit shows equivalent behavior for different number of species S, as long as N and L are fixed and S >> 1.

In this manuscript, we decided to focus on the case where NS = const since this ensures that all compared systems can build the same number of rings, NS/L. In this context, Figure 3D suggests that in order to achieve high yield, it is beneficial to build structures that are as homogeneous as possible. So, we believe that the relevance of our findings does not change. Nonetheless, in order to pick up the important point made by the referees, we added a remark to the manuscript and a figure that shows the flattening of the maximal yield for S >> 1.

– We have added a remark to subsection “Stochastic effects in the case of reduced resources”.

– We have added subsection “Influence of the heterogeneity of the target structure for fixed number of particles per species” and a figure to the Appendix.

3) Still in Figure 3, the characterization of the "stochastic yield catastrophe" seems incomplete in view of its relevance to the authors. It seems that there is a certain dependence ymax(N,L,S) (forgetting about Lnuc for now). Immediate questions that are not addressed are: how does ymax scale with N, for fixed L and S? How does it then scale with L? And how with the heterogeneity, L/S? One step in this direction is Figure 3D, in which the authors show ymax(S) for several values of L, but variable N (why then the label N=1000 in the panel?). More such plots should be presented to clarify the key properties of the system.

We are happy to take this opportunity to clarify the dependence of ymax(N,L,S) on N,L and S. As already explained in point (2) and shown by the added plot, ymax quickly becomes independent of S when S > 1. Hence, what remains to be described is the dependence of ymax on N and L. Phenomenologically, by increasing the particle number N (for fixed L), we encounter a threshold value Nth where the yield starts to increase sigmoidally from 0 to 1 (roughly over three orders of magnitude in N), quite similar to the behavior of the deterministic yield with respect to α and µ in Figure 1. Analogously, for fixed particle number N, by decreasing L, ymax exhibits a rapid transition from 0 to 1 at some threshold Lth. We included two additional plots in the Appendix that show the dependences of ymax on N and L, respectively. In the main text, in Figure 3C, we characterize the dependence of ymax on N and L in combined form as a “phase diagram” that shows the four regimes ymax = 0, 0 < ymax < 0.5, 0.5 < ymax < 0.99 and 0.99 < ymax in NL phase space (we now include the regime 0.99 < ymax to show that the transitions from 0 to almost perfect yield extend over a finite range in N or L, respectively).

In order to assess the relevance of stochastic effects for a particular system it is, therefore, crucial to know Nth(L), the threshold particle number as a function of L, whose inverse gives Lth(N). The function Nth(L) is exactly described by the line that separates the two regimes ymax = 0 and ymax > 0 in the “phase diagram” in Figure 3C. On the double logarithmic scale, we find that there is an approximate power-law dependence of ymax on L with an exponent of around 2.8. It is not totally clear to us why there is this dependence with such a high exponent. Simple scaling arguments that we tried could only explain an exponent of maximally 1.5. We investigate the underlying reasons for this dependence in a technical follow-up project (manuscript in preparation). The strongly nonlinear dependence of Nth on L partly explains why we argue that the expression “stochastic yield catastrophe” is justified, as the relevance of this effect strongly increases with the size of the target structure. We extended the caption of Figure 3C and the description in the main text in order to clarify the relevance of the stochastic yield catastrophe and the dependence of ymax on the system parameters.

The label “N = 1000” in Figure 3D was, of course, incorrect and misleading. We thank the attentive referee for this remark and have corrected the label.

– We have included the regime ymax > 0.99 in Figure 3C and extended the description in the Figure caption and the main text.

– We have added subsection “Dependence of the maximal yield ymax in the activation scenario on N and L” and two additional figures to the Appendix.

– We have corrected the false label in Figure 3D.

4) There are additional questions which relate to the relevance of the results in realistic situations. The model is elegant but very much simplified. For example, growth is fully reversible which may not be the case in many systems of macromolecular assembly. It is unclear how general and robust the results of the simplified model are if the simplifications are relaxed.

In our model two different modes of binding exist. Below a polymer size Lnuc binding is fully reversible whereas above Lnuc irreversible growths takes place. This irreversibility is meant to represent a separation of time scales between assembly and disassembly of a stable structure. Introducing a very small dissociation rate would not affect our results significantly on the time scale of observation (see answer to question 5) but makes the definition of the yield ambiguous. If binding is reversible for all unfinished structures (Lnuc = L − 1), trivially, only a yield of 1 is possible at the end of the assembly process. However, the process then may take an unphysical amount of time. The relevance of the size of Lnuc for the stochastic yield catastrophe is studied in Figure 3D. The robustness of all other effects with respect to the size of Lnuc (which corresponds to the degree of reversibility) is discussed in Appendix 4. Furthermore, the two additional models illustrated in Figure 6A and Figure 7A were introduced to test for robustness against other model assumptions like monomer binding and linear growth. In all cases, we still find stochastic yield reduction, which leads us to the conclusion that the observed effects do not rely on model specifics. We added remarks to the main text to clarify the relation between Lnuc and the existence of reversible and irreversible binding. Finally, we want to remark that the assumption of a finite Lnuc < L and irreversible growth above this threshold have been applied successfully in experimental studies, as for example of the Brome Mosaic Virus Capsids (Chevance and Hughes, 2008).

– We made the relation between reversible and irreversible binding and Lnuc more explicit in the main text.

5) The time needed for assembly is not discussed. In a biological context it is essential that assembly occurs sufficiently fast. High yield in the long term alone is not enough for assembly to happen. I understand that only looking at yield at long times is an elegant simplification of the problem, but it remains unclear how useful these insights are for real situations.

The importance of time in the biological context is a very pertinent observation and we fully agree with the referee that this question arises naturally: time efficiency plays a significant role in biological and also artificial self-assembly. Due to its significance and due to the following three reasons we believe the question of assembly time and time efficiency is beyond the scope of the current study. First, a thorough discussion of time in the context of self-assembly turned out to be extensive and hence would have considerably increased the length of the manuscript. Second, a simultaneous discussion of both factors, stochasticity and time efficiency, would complicate the analysis (if not treated separately) and obscure the significance and relevance of each one of these factors. Third, the final yield still is a most significant observable that characterizes the efficiency of the assembly process in an explicit way: Since the yield monotonically increases with time, the final yield represents an upper bound for the yield at any finite time. In other words, the final yield is the maximum yield that can be achieved in the assembly process irrespective of time constraints. As such, the final yield is an important and informative observable that can be defined and analyzed in an unambiguous way. For those reasons, we decided to focus here entirely on the analysis of the final yield and how it is affected by stochastic effects and to treat the time aspect separately.

However, in order to give the reader an impression of what the time evolution of the yield typically looks like, we now included a corresponding plot and a new section in the Appendix. The plot exhibits typical time trajectories of the yield in both scenarios and shows, in particular, that the yield increases monotonically with time. It also shows that typically, after some delay time, the yield increases rather quickly to a value that lies within 10% of the final yield and then continues to grow more slowly. This underlines that the final yield is a meaningful observable that describes the typical outcome of the assembly reaction well, under appropriate time constraints that are not too restrictive.

– We have added a paragraph to clarify the significance of the observable of the final yield and to state more explicitly that the time efficiency is not considered here.

– We have added subsection “Time evolution of the yield in the dimerization and activation scenario” and a figure to the Appendix.

https://doi.org/10.7554/eLife.51020.sa2

Article and author information

Author details

  1. Florian M Gartner

    Arnold Sommerfeld Center for Theoretical Physics (ASC) and Center for NanoScience (CeNS), Department of Physics, Ludwig-Maximilians-Universität München, München, Germany
    Contribution
    Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Project administration
    Contributed equally with
    Isabella R Graf and Patrick Wilke
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9801-4288
  2. Isabella R Graf

    Arnold Sommerfeld Center for Theoretical Physics (ASC) and Center for NanoScience (CeNS), Department of Physics, Ludwig-Maximilians-Universität München, München, Germany
    Contribution
    Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Project administration
    Contributed equally with
    Florian M Gartner and Patrick Wilke
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9169-9109
  3. Patrick Wilke

    Arnold Sommerfeld Center for Theoretical Physics (ASC) and Center for NanoScience (CeNS), Department of Physics, Ludwig-Maximilians-Universität München, München, Germany
    Contribution
    Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Project administration
    Contributed equally with
    Florian M Gartner and Isabella R Graf
    Competing interests
    No competing interests declared
  4. Philipp M Geiger

    Arnold Sommerfeld Center for Theoretical Physics (ASC) and Center for NanoScience (CeNS), Department of Physics, Ludwig-Maximilians-Universität München, München, Germany
    Contribution
    Conceptualization, Validation, Investigation, Visualization, Project administration
    Competing interests
    No competing interests declared
  5. Erwin Frey

    Arnold Sommerfeld Center for Theoretical Physics (ASC) and Center for NanoScience (CeNS), Department of Physics, Ludwig-Maximilians-Universität München, München, Germany
    Contribution
    Conceptualization, Resources, Supervision, Funding acquisition, Validation, Methodology, Project administration
    For correspondence
    frey@lmu.de
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8792-3358

Funding

Deutsche Forschungsgemeinschaft (GRK2062)

  • Patrick Wilke

Deutsche Forschungsgemeinschaft (QBM)

  • Florian M Gartner
  • Isabella R Graf

Aspen Center for Physics (PHY-160761)

  • Erwin Frey

Deutsche Forschungsgemeinschaft (EXC-2094 - 390783311)

  • Erwin Frey

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Nigel Goldenfeld for a stimulating discussion, and Raphaela Geßele and Laeschkir Hassan for helpful feedback on the manuscript. This research was supported by the German Excellence Initiative via the program ‘NanoSystems Initiative Munich’(NIM) and was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094–390783311. FMG and IRG are supported by a DFG fellowship through the Graduate School of Quantitative Biosciences Munich (QBM). We also gratefully acknowledge financial support by the DFG Research Training Group GRK2062 (Molecular Principles of Synthetic Biology). Finally, EF thanks the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611, for their hospitality and inspiring discussions with colleagues.

Senior Editor

  1. Naama Barkai, Weizmann Institute of Science, Israel

Reviewing Editor

  1. Frank Jülicher, Max Planck Institute for the Physics of Complex Systems, Germany

Reviewer

  1. Pablo Sartori, The Rockefeller University, United States

Publication history

  1. Received: August 12, 2019
  2. Accepted: February 4, 2020
  3. Accepted Manuscript published: February 5, 2020 (version 1)
  4. Version of Record published: March 23, 2020 (version 2)

Copyright

© 2020, Gartner et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,266
    Page views
  • 213
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Physics of Living Systems
    William T Ireland et al.
    Research Article Updated

    Advances in DNA sequencing have revolutionized our ability to read genomes. However, even in the most well-studied of organisms, the bacterium Escherichia coli, for ≈65% of promoters we remain ignorant of their regulation. Until we crack this regulatory Rosetta Stone, efforts to read and write genomes will remain haphazard. We introduce a new method, Reg-Seq, that links massively parallel reporter assays with mass spectrometry to produce a base pair resolution dissection of more than a E. coli promoters in 12 growth conditions. We demonstrate that the method recapitulates known regulatory information. Then, we examine regulatory architectures for more than 80 promoters which previously had no known regulatory information. In many cases, we also identify which transcription factors mediate their regulation. This method clears a path for highly multiplexed investigations of the regulatory genome of model organisms, with the potential of moving to an array of microbes of ecological and medical relevance.

    1. Physics of Living Systems
    Sean Fancher, Andrew Mugler
    Research Article

    Morphogen profiles allow cells to determine their position within a developing organism, but not all morphogen profiles form by the same mechanism. Here we derive fundamental limits to the precision of morphogen concentration sensing for two canonical mechanisms: the diffusion of morphogen through extracellular space and the direct transport of morphogen from source cell to target cell, e.g., via cytonemes. We find that direct transport establishes a morphogen profile without adding noise in the process. Despite this advantage, we find that for sufficiently large values of profile length, the diffusion mechanism is many times more precise due to a higher refresh rate of morphogen molecules. We predict a profile lengthscale below which direct transport is more precise, and above which diffusion is more precise. This prediction is supported by data from a wide variety of morphogens in developing Drosophila and zebrafish.